Characterization of Microbiota Composition and Presence of Selected Antibiotic Resistance Genes in Carriage Water of Ornamental Fish

International trade with ornamental fish is gradually recognized as an important source of a wide range of different antibiotic resistant bacteria. In this study we therefore characterized the prevalence of selected antibiotic resistance genes in the microbiota found in the carriage water of ornamental fish originating from 3 different continents. Real-time PCR quantification showed that the sul1 gene was present in 11 out of 100 bacteria. tet(A) was present in 6 out of 100 bacteria and strA, tet(G), sul2 and aadA were present in 1–2 copies per 100 bacteria. Class I integrons were quite common in carriage water microbiota, however, pyrosequencing showed that only 12 different antibiotic gene cassettes were present in class I integrons. The microbiota characterized by pyrosequencing of the V3/V4 variable region of 16S rRNA genes consisted of Proteobacteria (48%), Bacteroidetes (29.5%), Firmicutes (17.8%), Actinobacteria (2.1%) and Fusobacteria (1.6%). Correlation analysis between antibiotic resistance gene prevalence and microbiota composition verified by bacterial culture showed that major reservoirs of sul1 sul2, tet(A), tet(B) tet(G), cat, cml, bla, strA, aacA, aph and aadA could be found among Alpha-, Beta- and Gammaproteobacteria with representatives of Enterobacteriaceae, Pseudomonadaceae, Rhizobiaceae and Comamonadaceae being those most positively associated with the tested antibiotic resistance genes.


Introduction
International trade with ornamental fish is gradually being recognized as an important source of a wide range of different antibiotic resistant bacteria [1,2] since antibiotics in ornamental fish breeding are used both for therapy and as prophylactic treatment before transport [3]. This results in frequent and efficient positive selection of antibiotic resistant bacterial clones which, during periods of antibiotic withdrawal, may serve as reservoir of antibiotic resistance genes. Once an antibiotic resistant clone has been selected by antibiotic use, conditions in ornamental fish rearing, i.e. ambient temperature and aquatic environment, allow for horizontal transfer of antibiotic resistance genes to other population members by genetic elements of varying mobility, e.g. conjugative plasmids, transposons or integrons [4,5]. The mobility of elements transferring antibiotic resistance is not limited by species or genus and genetic elements with the same antibiotic resistance genes can be found across a broad range of different bacterial species [6,7]. Although absolute volumes of carriage water are quite small, frequent use of antibiotics, international trade, close contact to humans and bacteria adapted to an aquatic environment may result in the long term survival of such bacteria in new ecological niches where they may serve as a source of new (combinations of) antibiotic resistance genes for 'domestic' microbiota.
One of the methods used traditionally for assessing antibiotic resistance in complex bacterial populations is bacterial culture on nutrient agars with and without antibiotics. However, this provides only information on the antibiotic resistance in bacterial species which can grow under the given conditions and what may introduce certain bias. This is why culture-independent techniques such as quantitative real-time PCR or next generation sequencing are therefore increasingly used for characterizing the prevalence of a particular gene in a given bacterial community [8][9][10]. In addition, these techniques also allow for taxonomic characterization of the bacterial community itself [11,12] although they cannot differentiate between live and dead cells. In this study we therefore tested samples of ornamental fish carriage water originating from 3 different continents immediately after their import to the Czech Republic with the aim to assess the prevalence of selected antibiotic resistance genes and characterize the microbiota composition in these samples. Since initial data indicated an increased abundance of integrons, we also determined the integron gene cassettes present in the microbiota of ornamental fish carriage water in detail. Finally we predicted the association of particular taxa with antibiotic resistance and confirmed the predictions by culture and antibiotic resistance characterization of isolates from selected bacterial families.

Results
Quantitative real-time PCR sul1 and tet(A) were the most prevalent genes out of those tested during initial sample screening. The median prevalence of sul1 in microbiota of ornamental fish carriage water was 0.11 which corresponds to the presence of the sul1 gene in 11 out of 100 bacteria (Fig. 1). The median prevalence of tet(A) in the microbiota of ornamental fish carriage water was 0.06 corresponding with the presence of tet(A) in approx. 6 out of 100 bacteria. The prevalence of strA, tet(G), sul2 and aadA was between 1-2 copies per 100 bacteria. The least prevalent were cat and tet(B) genes being present in less than 1 out of 10,000 bacteria. However, although cat and tet(B) were detected at a low prevalence on average, there were samples in which these genes were present at the same prevalence as the rest of the target genes ( Fig. 1).

Pyrosequencing of integron gene cassettes
Since the sul1 and aadA genes were present in ornamental fish water microbiota at a high prevalence (Fig. 1) and as both these genes are present in class I integrons [5,7], we hypothesized that class I integrons must have been common in genomes of microbiota of ornamental fish carriage water. To address this hypothesis, we performed PCR with primers specific for 59CS and 39CS ends of class I integrons and this PCR resulted in positive amplification in 41 out of 48 samples (Fig. 2).
To identify gene cassettes present in class I integrons, the resulting PCR products were pooled and pyrosequenced. Pyrosequencing resulted into 223,877 reads longer than 200 bp. However, 188,201 reads were excluded by GS De Novo Assembler as containing repeated, chimeric or other corrupted sequences. Isotigs were therefore assembled from 18,806 reads and the remaining 16,870 reads represented singletons which, when submitted to GenBank, did not match any entry.
Altogether 33 different genes were detected in integron structures. Out of these, qacE, dfrA22 and aadA13 were the most frequent whereas ereA, arr-8 or qnrVC1 were detected only sporadically (Table 1). Only two integron gene cassettes with predicted functions, nit2 and estX, were not associated with antibiotic resistance. nit2 is a member of the nitrilase family which acts as a tumor suppressor in human cells [13] and estX encodes esterase X of unknown function found in integrons [14]. When the genes were grouped based on the type and mechanism of resistance, there were only 12 different antibiotic gene cassettes present in class I integrons (qacE, dfr, aadA, aadB, aacA, aphA, cat, cmlA, blaP1, ereA, arr8 and qnr).
In the next step we verified the pyrosequencing results performed on the pooled sample by real-time PCR quantification of aacA, ereA, dfrA22, cmlA1, blaP1, aphA, aadA, qnrVC1 and arr-8 using individual samples. The real-time PCRs were performed with two different template DNAs, either with 5CS-3CS PCR products or with carriage water DNA. Real-time PCR using the former template allowed us to verify the results of pyrosequencing as the same template was subjected to the analysis. The results of real-time PCR using the latter template allowed us to predict which antibiotic resistance genes were dominantly associated with integrons and which were associated with genetic elements different from integrons. When average values were calculated from individual samples after real-time PCR using 5CS-3CS PCR products as templates, similar results to those obtained by pyrosequencing were recorded (Fig. 3). Performing real-time PCR with carriage water DNA as a template showed that aacA and ereA genes must have been commonly encoded by genetic elements different from integrons as their relative representation in total carriage water DNA was higher than in amplified 5CS-3CS PCR products. On the other hand, the lower representation of dfrA genes in carriage water DNA when compared with amplified 5CS-3CS PCR products indicated that this gene must have been common in integrons and was not commonly encoded by other genetic elements (Fig. 3).

Microbiota composition determined by 16S rRNA pyrosequencing
Because of the high antibiotic resistance gene prevalence in the microbiota of ornamental fish carriage water, next we were interested in the microbiota composition in carriage water and predicting the bacterial taxa forming a reservoir of antibiotic resistance in such populations.
Knowing the microbiota composition and antibiotic resistance gene prevalence for each of the 48 samples, in the next step we correlated these two parameters with each other. Into this analysis we included only the families which were present in more than a half of all samples, i.e. we excluded minority population members which provided false results when included. There were 3 main family groups according to their relationship to the prevalence of tested antibiotic resistance genes. The first group comprised 11 bacterial families which exhibited a negative correlation with the majority of antibiotic resistance genes quantified in this study (Fig. 5). Out of these, Lachnospiraceae, Ruminococcaceae, Lactobacillaceae, Clostridiaceae I, Peptostreptococcaceae, Rickenellaceae and Bacteroidaceae are common to the intestinal tract of warmblooded animals and fish [11,[15][16][17] which may explain their origin in the carriage water.
The second cluster comprised of 29 bacterial families which exhibited a moderate positive correlation with the majority of antibiotic resistance genes tested in this study (Fig. 5). Representatives of Enterobacteriaceae and Moraxellaceae were clustered to this group.
The third group included 9 families which exhibited a high positive correlation with the majority of antibiotic resistance genes tested in this study. Representatives of Pseudomonadaceae, Rhizobiaceae, Commamonadaceae and Flavobacteriaceae were clustered to this group (Fig. 5).
In the last experiment we experimentally tested the predicted correlations in representatives of families Pseudomonadaceae (28 isolates), Rhizobiaceae (n = 4), Flavobacteriaceae (n = 8), Comamonadaceae (n = 22), Enterobacteriaceae (n = 17) and Moraxellaceae (n = 16). A new set of 42 carriage water samples had to be used for culture because the original samples were used only for DNA purification and not for bacterial culture. All the isolates were tested both for the presence of target genes and for the phenotypic expression of antibiotic resistance by disk diffusion assay.
Antibiotic resistance genes targeted in this study were detected the most frequently in isolates belonging to the families Pseudomonadaceae, Enterobacteriaceae and Rhizobiaceae, although for Rhizobiaceae the conclusion was based only on 4 isolates. Isolates of families Comamonadaceae and Moraxellaceae encoded the tested antibiotic resistance genes the least frequently and with a single exception, none of the tested antibiotic resistance genes were confirmed in Flavobacteriaceae (Fig. 6). Except for Flavobacteriaceae, bacterial culture followed by gene detection confirmed the predictions based on the correlation of microbiota composition and antibiotic gene prevalence.
Antibiotic resistance tested in the same isolates by disk diffusion assay was extremely high. More than half of the isolates in all families were resistant to ampicillin, tetracycline, and sulphonamides. Resistance to tetracycline was the most common because between 75 to 100% of isolates in all 6 families were resistant to its activity. Isolates in families Pseudomonadaceae and Flavobacteriaceae exhibited the highest levels of antibiotic (multi)resistance out of the 6 families compared. In all families, but in Flavobacteriaceae in particular, it could be seen that the resistant phenotypes were only partially explained by the presence of tested genes and the resistances must have been caused genes different from those targeted in this study (Fig. 6).

Discussion
The prevalence of strA, sul1, sul2, tetA, tetG, aadA, aacA and cml in microbiota of ornamental fish carriage water was very high and the sul1 gene was present in more than 10% of all bacteria. This prevalence was 2-4 logs higher than in the feces of farm animals, which we characterized recently [10].
Despite a wide geographical distribution of sample origins, pyrosequencing of class I integron gene cassettes identified only a limited pool of genes. DqacE, dfrA and aadA genes were among those recorded most frequently, similar to previous findings [5,7,18]. On the other hand, we never recorded tetracycline resistance genes as part of integrons despite the fact that tet(A) and tet(G) genes were quite common in carriage water microbiota. The comparison of individual gene prevalence in total DNA and amplified integron sequences also showed that certain genes were common to class I integrons but were infrequently encoded outside of integrons. dfrA22 coding for trimethoprim resistance is such an example as it was common in integrons but underrepresented in total DNA. On the other hand, the aac gene responsible for aminoglycoside resistance was a very common gene in the microbiota of ornamental fish but this gene was infrequently associated with integrons. This corresponds to studies of Hopkins et al. and Galimand et al. who showed that aac genes are a common part of transposons and conjugative plasmids separate from class I integrons [19,20]. Interestingly, since class I integron gene cassettes were identified in this study without any bias and because these genes positively correlated with families belonging to Alpha-, Betaand Gammaproteobacteria, this genetic element seems to be restricted to these classes belonging to phylum Proteobacteria and only rarely or not at all spread to bacterial species of phyla Firmicutes and Bacteroidetes.
Since an environment consisting of closed transport bags with ornamental fish was analyzed, the carriage water must have been influenced by fish gut microbiota. However, the composition of ornamental fish carriage water samples was similar to that of fresh water aquatic environments [12,21], and mostly differed from fish gut microbiota [16,22,23], sea water microbiota, or even fecal material of warm blooded animals [11,15,17,24,25].
Correlation analysis between microbiota composition and prevalence of selected antibiotic resistance genes was performed to predict taxa likely serving as the most important antibiotic resistance gene reservoirs. When performing such an analysis, a certain number of limiting points must be taken into consideration. Firstly, even of the taxa harboring antibiotic resistance genes, not all isolates are resistant and absolute correlation in acquired antibiotic resistances cannot be expected. Secondly, different microbiota members may code for the same antibiotic resistance gene thus making the contribution of individual microbiota members rather difficult to determine. Thirdly, under-represented taxa may carry an antibiotic gene at a high frequency but this will unlikely cause a change in global antibiotic resistance gene prevalence and their positive correlation with a particular antibiotic resistance gene may remain unrecognized. Considering all these possibilities and also the fact that the predictions and bacterial culture verifications were performed in 2 different sets of samples, the fact that tet(A), tet(G) or sul1 genes were the most frequently identified in Pseudomonadaceae and Rhizobiaceae and less frequently in Comamonadaceae is in agreement with the correlation analysis. On the other hand, the widespread presence of strA among different taxa made predictions on its preferential association with a particular taxon less accurate. The correlation analysis also indicated that none of the families for which we collected isolates should be associated with dfr genes and in agreement, this gene was not detected in any single bacterial isolate. Finally, it would be mistaken to conclude that taxa which were identified in this study as negatively correlating with the tested genes represent an antibiotic sensitive part of a bacterial population. Though we did not prove this experimentally, results recorded for Flavobacteriaceae rather suggest that the majority of taxa forming carriage water microbiota are antibiotic resistant but due to the presences of genes different from those tested in this study.
In conclusion, the microbiota of carriage water from ornamental fish is commonly resistant to antibiotics. Moreover, since the same antibiotic resistance genes were present in bacteria belonging to different classes, ornamental fish breeding represents an ecological niche where horizontal gene transfer must be quite common. Predictions derived from culture-independent protocols allowed for the identification of the most frequent reservoirs of tested antibiotic resistance genes and together with culture-based confirmation indicated that horizontal gene transfer is limited to the class level and that the spread of the same genes among bacterial phyla is infrequent. Although absolute volumes of carriage water are quite small and acquired antibiotic resistance is small in scope when compared with intrinsic and natural resistances of certain bacterial species, ornamental fish breeding and trade represent an underestimated ecological niche with common presence of antibiotic resistance. The selection pressure then leads to selection of new combination of antibiotic resistance which can be transferred to domestic microbiota and/or to pathogens.

Ethics Statement
Owner of the company was aware that the samples were being collected for this study and gave their permission for doing so. Such permissions can be obtained from AC who, as a coauthor, was responsible for carriage water samples collection.

Sample collection and DNA isolation
Forty-eight samples of ornamental fish carriage water were collected immediately after import to the Czech Republic. Countries of origin included Brazil (n = 1), China (n = 2), Hong Kong (n = 3), Indonesia (n = 10), Israel (n = 3), Nigeria (n = 1), Peru (n = 5), Singapore (n = 9), Thailand (n = 2) and Vietnam (n = 12). The bacteria were pelleted from the carriage water by centrifugation at 30006g for 15 min at 4uC and the DNA was immediately extracted using QIAamp DNA Stool Mini kit following the manufacturer's protocol (Qiagen). The quality and concentration of the purified DNA was determined spectrophotometrically and DNA was stored at 220uC until use. Table 1 Quantitative real-time PCR sul1, sul2, strA, tet(A), tet(B), tet(G), cat and aadA genes were quantified by real-time PCR using primers listed in Table 2. Amplification of 16S rRNA genes was used as a reference to determine the total amount of bacterial DNA in each sample. PCR was performed in 3ml volumes in 384-well microplates using QuantiTect SYBR Green PCR Master Mix (Qiagen) and a Nanodrop pipeting station (Innovadyne) for dispensing PCR mixtures. PCR and signal detection was performed with a LightCycler II (Roche) with an initial denaturation at 95uC for 15 min followed by 40 cycles of denaturation at 95uC for 20 s, primer annealing at 61uC for 30 s and extension at 72uC for 30 s. Each sample was subjected to real-time PCR in triplicate and the mean values of the triplicates were used for subsequent analysis.
The Ct values of genes of interest were normalized to an average Ct value of 16S rRNA gene amplification (DCt) and the relative amount of each gene of interest in given bacterial population was finally calculated as 2 2DCt .
16S rRNA gene PCR, pyrosequencing and data analysis DNA samples from carriage water were individually amplified over the V3/V4 variable region of 16S rRNA [11]. PCR was performed using HotStarTaq Plus Master Mix Kit (Qiagen, Hilden, Germany) with initial denaturation at 95uC for 15 min followed by 30 cycles of denaturation at 94uC for 40 s, primer annealing at 55uC for 45 s and extension at 72uC for 1 min. The PCR products were separated on a 1.5% agarose gel and extracted from the gel with QIAquick Gel Extraction kit (Qiagen, Hilden,    . Antibiotic resistances in isolates belonging to families predicted as highly associated with antibiotic resistance in the microbiota of carriage water of ornamental fish. Data shows the percentage of isolates within a given family encoding the gene indicated on the X-axis out of all isolates tested. Black columns, results based on PCR detection of gene as indicated. Grey columns, percentage of antibiotic resistance determined in the same isolates by disk diffusion assay (ampicillin resistance is shown with bla gene, the same tetracycline resistance data are shown in tet(A), tet(B) and tet(G) genes, streptomycin resistance is shown in strA, aacA, aph and aadA genes, chloramphenicol resistance is shown in cml and cat genes, sulphonamide resistance is shown in sul1 and sul2 genes and nalidixic acid resistance is shown in qnr gene. Resistance to erythromycin (ere), trimethoprim (dfr) and rifampicin (arr) was not determined by disk diffusion assay and therefore are not shown. doi:10.1371/journal.pone.0103865.g006 Germany). Purified PCR products were subjected to library preparation using the GS Amplicon Library Preparation Kit (Roche, Basel, Switzerland), emulsion PCR and pyrosequencing with the GS Junior Titanium kit strictly following the manufacturer's instructions (Roche, Basel, Switzerland). The pyrosequencing was performed with the GS Junior 454 sequencer (Roche, Basel, Switzerland).
Data analysis was performed as described previously [11]. In brief, fasta and qual files generated as an output of the pyrosequencing were uploaded into Qiime software [26]. Quality trimming criteria included no mismatch in MID sequences and a maximum of 1 mismatch in primer sequences. The obtained sequences with a qual score higher than 20 were all shortened to a length of 350 bp and classified with RDP Seqmatch with an OTU discrimination level set to 97%. Diversity analyses on OTU clusters were performed both using all sequences available for each sample and using the same number of randomly selected sequences adjusted to the number of sequences available for the sample with the lowest coverage. Finally, UniFrac analysis [27] followed by principal coordinate analysis (PCoA) was used to characterize the diversity in the microbial populations tested.

Integron PCR, pyrosequencing and data analysis
Gene cassettes present in integrons were PCR amplified using 59CS and 39CS primers [28] and Taq PCR Master Mix Kit (Qiagen, Hilden, Germany). The PCR cycling conditions consisted of initial denaturation at 94uC for 3 min followed by 35 cycles of denaturation at 92uC for 40 s, primer annealing at 55uC for Table 2. List of primers used in this study.

Verification of predicted associations of particular bacterial families and antibiotic resistance genes
To verify the predicted associations of particular bacterial families and antibiotic resistance genes, an additional 42 carriage water samples (i.e. in addition to 48 samples mentioned above) were collected and used for culture on Columbia blood agar, McConkey agar, Anacker and Ordal agar without any antibiotics at 25uC for 48-72 hours. Five morphologically different colonies were picked from each plate for species identification using MALDI Biotyper (Bruker Daltronik, Bremen, Germany), as recommended by the manufacturer. Phenotypic resistance to ampicillin, chloramphenicol, nalidixic acid, streptomycin, sulfonamides and tetracycline was tested by the disk diffusion method. Since interpretative criteria have not been defined for fish bacteria, the evaluation of results was based only on the presence/absence of inhibition zones.
DNA from 95 successfully identified colonies was extracted by boiling for 15 min and used as DNA template in real-time PCR with primers listed in Table 2. Real time PCR was used to maintain the same conditions, e.g. primer pairs or PCR mix, as in all previous PCR. In this case, real time PCR was not used for any quantification. PCRs in which the Ct values of the antibiotic resistance gene differed from the Ct value of 16S rRNA gene by less than 4 cycles and the melting temperature after the denaturation step following PCR was within 0.5uC of the expected range were considered positive whilst all other results were considered negative. The rather high 4 cycle difference in Ct values was used considering the different copy numbers of rRNA genes in different species as well as the potential presence of antibiotic resistance genes in (multicopy) plasmids.

Statistics
Pyrosequencing results from 16S rRNA gene sequencing were analyzed by UniFrac analysis followed by principal coordinate analysis (PCoA) present in Qiime software. The Spearman's rank correlation of microbiota with the prevalence of particular antibiotic resistance genes was visualized as a cluster heatmap using Ward's method. The analysis was performed with MATLAB version 2013a (MathWorks, Natick, MA, USA).