Bacteria From the Multi-Contaminated Tinto River Estuary (SW, Spain) Show High Multi-Resistance to Antibiotics and Point to Paenibacillus spp. as Antibiotic-Resistance-Dissemination Players

Bacterial resistance to antibiotics is an ever-increasing phenomenon that, besides clinical settings, is generally assumed to be prevalent in environmental soils and waters. The analysis of bacteria resistant to each one of 11 antibiotics in waters and sediments of the Huelva’s estuary, a multi-contaminated environment, showed high levels of bacteria resistant mainly to Tm, among others. To further gain knowledge on the fate of multi-drug resistance (MDR) in environmental bacteria, 579 ampicillin-resistant bacteria were isolated tested for resistance to 10 antibiotics. 92.7% of the isolates were resistant to four or more antibiotic classes, indicating a high level of multi-resistance. 143 resistance profiles were found. The isolates with different MDR profiles and/or colony morphologies were phylogenetically ascribed based on 16S rDNA to phyla Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes, including 48 genera. Putative intrinsic resistance was detected in different phylogenetic groups including genera Altererythrobacter, Bacillus, Brevundimonas, Erythrobacter, Mesonia, Ochrobactrum, and Ponticaulis. Correlation of the presence of pairs of the non-intrinsic-resistances in phylogenetic groups based on the kappa index (κ) highlighted the co-habitation of some of the tested pairs at different phylogenetic levels. Maximum correlation (κ = 1.000) was found for pairs CzR/TcR in Betaproteobacteria, and CcR/TcR and EmR/SmR in Sphingobacteriia at the class level, while at the genus level, was found for CcR/TcR and NxR/TmR in Mesonia, CzR/TmR and EmR/KmR in Paenibacillus, and CcR/EmR and RpR/TcR in Pseudomonas. These results could suggest the existence of intra-class and intra-genus-transmissible genetic elements containing determinants for both members of each pair. Network analysis based on κ values higher than 0.4 indicated the sharing of paired resistances among several genera, many of them centered on the Paenibacillus node and raising the hypothesis of inter-genera transmission of resistances interconnected through members of this genus. This is the first time that a possible hotspot of resistance interchange in a particular environment may have been detected, opening up the possibility that one, or a few, bacterial members of the community could be important promoters of antibiotic resistance (AR) dissemination in this environment’s bacterial population. Further studies using the available isolates will likely give insights of the possible mechanisms and genetic elements involved.


INTRODUCTION
The antibiotic resistance (AR) of bacterial pathogens is currently a worldwide problem with severe consequences for the treatment of infectious diseases (World Health Organization [WHO], 2018). However, AR existed in prehistoric (Perry et al., 2016) and more recent times (Devault et al., 2017); and its origin is not directly related to the clinical use of antibiotics. It was probably triggered by the appearance of antibiotics in the environment when their producers needed to be protected from their effects (Jiang et al., 2017). The environmental role of these compounds is not yet clear and for a long time was considered only ecological. This function was demonstrated only in a few cases (Currie et al., 1999;Neeno-Eckwall et al., 2001;Haas and Défago, 2005). Antibiotic concentrations in nature may not be enough for this purpose (Fajardo et al., 2009), while low concentrations have pleiotropic effects that are not related to competition (Davies et al., 2006;Linares et al., 2006;Yim et al., 2007;Fajardo and Martinez, 2008;Sengupta et al., 2013;Andersson and Hughes, 2014;Goneau et al., 2015;Martínez, 2017;Xiong et al., 2017). However, low individual concentrations of a cocktail of antibiotics could have a combined effect [reviewed by Danner et al. (2019)]. Additionally, sublethal concentrations select resistant phenotypes from mutations or horizontal gene transfer (HGT) that sometimes use mobile genetic elements (MGEs) (Kohanski et al., 2010;Baharoglu and Mazel, 2011;Gutierrez et al., 2013;Haaber et al., 2017;Lekunberri et al., 2017). The AR in the environment is a natural process but is also promoted by anthropogenic pollutants like antibiotics, biocides, heavy metals, hydrocarbons, pesticides, and nanomaterials, among others (Baker-Austin et al., 2006;Dealtry et al., 2014;Pal et al., 2015Pal et al., , 2017Chen et al., 2017;Poole, 2017;Wang et al., 2018;Nguyen et al., 2019). Therefore, naturally antibiotic-resistant bacteria (ARB) could be a prime source for antibiotic resistance genes (ARGs) found in pathogens as certain findings may indicate (Rodríguez et al., 2004;Poirel et al., 2005Poirel et al., , 2012Yang J. et al., 2013;Jiang et al., 2017), but may later have evolved and be selected.
The interaction of microorganisms from natural, clinical, and farming sources could select multi-resistant bacteria (MRB), severely affecting environments, change biodiversity, and modify evolution paths in favor of resistants (Ding and He, 2010;Gupta, 2011;Proia et al., 2013;Laverman et al., 2015;Hiltunen et al., 2017;Grenni et al., 2018). The major challenge posed by ARB is their acquisition of MDR. Several mechanisms for the accumulation of resistances are known (Cambray et al., 2010;Michael et al., 2012;Brown-Jaque et al., 2015;Carraro et al., 2016Carraro et al., , 2017Ramsay and Firth, 2017). However, MDR spreads in the environment and the resistance patterns that are passed-on together or transferred to pathogenic bacteria, or how these bacteria can model their environments remain unknown (Chamosa et al., 2017). Studies to estimate the risk of HGT between environmental and pathogenic bacteria have identified hotspots of transfer including wastewater treatment plants (WWTPs), which are recognized as one of the main sources of antibiotics pollution in surface waters (Guo et al., 2017;Manaia et al., 2018). Few experimental studies on the specific transmissible elements and their dissemination capabilities have been conducted. The structure, function, and inter-relationships of the environmental microbiomes are complex phenomena only partially studied and their comprehension needs to use multidisciplinary approaches and methodologies, including different culture-dependent and independent approaches, to analyze in a "one health" approach how the AR affects environmental microbial biodiversity and limits the use of antibiotics (Stefani et al., 2015;Ellington et al., 2017;Hiltunen et al., 2017;Tripathi and Cytryn, 2017).
Recent studies have shown that many rivers worldwide contain relevant concentrations of antibiotics Wilkinson and Boxall, 2019) up to 50 µg/L in some African countries or up to 10 µg/L in European ones (Danner et al., 2019). The prescriptions of antibiotics in Spain has risen from 2012 to 2018 (European Centre for Disease Prevention and Control [ECDC], 2017) 1 and several studies on the concentration of antibiotics in Spain's surface waters have shown concentrations of 1.3 µg/L on average [reviewed in Danner et al. (2019)]. Running surface waters such as rivers and their estuaries in anthropogenically influenced areas are a means for antibiotic resistance dissemination through their courses and receiving waters. Some rivers receive effluents from WWTPs facilities, which are not able in most cases to completely eliminate pharmaceuticals such as antibiotics, ARBs, or ARGs, even with the implementation of disinfection protocols (Loos et al., 2012;Suzuki et al., 2017;Barancheshme and Munir, 2018;Ju et al., 2019). So, these bacteria and chemicals carried by the treated wastewater can affect the environment by selection of antibiotic resistant bacteria, changing biodiversity and increasing the levels of resistant and MRB (Sanseverino et al., 2018). Also soils and sediments can receive antibiotics and other pollutants that are not completely degraded and can impact their microbial activity and diversity (Cycon et al., 2019). In addition, other substances can be contaminating rivers and estuaries when they receive effluents from industries established within their borders.
The Huelva estuary (SW, Spain) receives waters from two acidic and heavy metal-rich rivers that mix with seawater from the Atlantic Ocean. The WWTP of Huelva, a city with about 145,000 inhabitants, processes urban wastewaters about 240.000 population equivalents 2 ; and emits its effluents into the estuary. Also, effluents from a number of industries around the estuary, including oil refineries, metallurgic industries, fertilizers plants, etc., (Figure 1) arrive in the estuary to create a multicontaminated environment (Davis et al., 2000;Sainz et al., 2004;Pérez-López et al., 2011;Amils, 2016;Papaslioti et al., 2018). Moreover, some studies have reported the presence of radioactivity in the estuary that is related to the fertilizer production that uses uranium-containing phosphogypsum, which accumulates in sediments (Martínez-Aguirre and García-León, 1997;Elbaz-Poulichet et al., 1999). As stated above, several studies have reported that the appearance of antibiotic resistance can be facilitated by the presence in the environment of heavy metals that are abundant in the Huelva estuary; hydrocarbons that come from oil refineries' management or accidental spills; and for antibiotics that are likely provided by the WWTP. Due to the prevalence of and magnitude of the metals pollution, previous studies have been limited to this aspect.
This study reports our conclusions as to the prevalence of culturable ARB in the surface water and sediments of the Huelva estuary with an aim to contribute to the current knowledge on the diversity and fate of MDR of culturable environmental bacteria as well as to collect biological materials for further study on the transmission of antibiotic resistance and the MGEs involved. 579 ampicillin-resistant bacteria were isolated and their MDR levels and resistance profile diversity were determined. The analysis of these data identified putative intrinsic resistances in different phylogenetic groups as well as correlations of the presence of pairs of acquired resistances mainly at the genus level. A network analysis of the best correlations suggests that members of the Paenibacillus genus could be important players in the dissemination of resistance in this environment.

Huelva's Estuary Sampling
Samples were collected in areas named H1 and H2 (Figure 1)

Culture and Isolation of ARB
In order to quantify total culturable and resistant bacteria, dilutions in the culture media (10 0 -10 −6 ) of the liquid samples and culture-media-suspended sediments were inoculated by extension in Petri dishes containing autoclaved solid media. Two culture medias were used: marine [55.1 g/l Difco Marine agar (BD, Sparks, MD, United States)] and nutritive [3 g/l meat extract (Merck, Darmstadt, Germany), 5 g/l NaCl (Merck), 10 g/l bactopeptone (Labs. Conda, Madrid, Spain), 15 g/l agar (Labs. Conda)]. After autoclaving (121 • C, 20 , 1 atm), media were left to cool down to 50 • C and when needed, supplemented with one out of 11 antibiotics at the following final concentrations: 100 µg/ml for streptomycin (Sm) (Duchefa Biochemie, Haarlem, Netherlands); 50 µg/ml for ampicillin (Ap) and vancomycin (Vm) (Duchefa), ceftazidime (Cz), kanamycin (Km), and nalidixic acid (Nx) (Sigma-Aldrich, St. Louis, MO, United States); 25 µg/ml for chloramphenicol (Cc) and erythromycin (Em) (Duchefa), and tetracycline (Tc) (Sigma); 16 µg/ml for trimethoprim (Tm) (Duchefa); and 8 µg/ml for rifampicin (Rp) (Duchefa). To all the media, cycloheximide (Duchefa) was added at the final concentration of 75 µg/ml. Three dishes per dilution were inoculated and cultured at 30 • C for 5 days before the colony counts were done. Bacteria growing on Ap-containing media (Ap R ) and showing different colony morphologies were selected for isolation. When many colonies with the same morphology appeared from a sample, one out of every three was collected. Isolation from the initial plating colonies was performed by successive passages on solid media containing the same selective agent until isolated colonies of the same morphology were obtained during three consecutive replatings. Isolates were stored at −80 • C in 50% (v/v) of glycerol (87%, Merck) and the corresponding liquid culture medium.
Multi-resistance of the isolates was tested on the solid medium used for their isolation containing individual antibiotics at the concentrations indicated above. Dishes inoculated with Ap served as growth control.

DNA Extraction
DNA was isolated from individual bacterial colonies according to Kai et al. (2006), or using the Ultraclean Microbial DNA Isolation Kit (MOBIO, Carlsbad, CA, United States) according to the manufacturer's instructions [except for a pretreatment with lysozyme 10 mg/ml (Sigma) for 1 h at 37 • C].
Isolates were initially assigned to the genus or species with the highest sequence similarity using NCBI BLAST (Altschul et al., 1990); and were further identified compared with type strains in the phylogenetic trees built with the sequences at the Ribosomal Database Project (RDP) (RRID: SCR_006633) (Cole et al., 2014). Clustal X 4 was used for alignments and building of neighbor-joining trees, with 1000 repetitions for bootstrapping.

Statistical Analysis
Data were analyzed using SPSS Statistics V21.0 software (IBM Corporation, New York, NY, United States) and Microsoft Excel 2013. The normality of the distribution of the variables was evaluated by using the Kolmogorov-Smirnov test.
To characterize MDR level the multiple-antibiotic resistance index (MAR) was calculated for each isolate and for the groups of them considered as the median (x) of the isolates' MARs (Krumperman, 1983).
The non-parametric tests of Mann-Whitney and Kruskal-Wallis were performed to compare the MARs median. Statistical correlations between the sampling areas or culture media with the observed resistance to each antibiotic were analyzed using contingency tables and contrasting Chi-square.
The kappa indexes (κ) for association of resistance pairs were calculated for groups with eight or more isolates. Networking among pairs with κ0.4 by class and genus were represented using gephi 5 .
For all analysis, results were considered statistically significant for p < 0.05.
Distributions of resistants among zones or phylogenetic groups are presented by using the Circos online application 6 .
Relative abundances of resistance to each antibiotic (as% of resistants/total culturable bacteria) (Figure 2) were, in general, higher for bacteria tested on the marine rather than on the nutritive medium. The highest was obtained for Tm-resistants, irrespective of the sampling zone or the culture medium. However, Tc-resistants were common among the bacteria grown on marine medium but scarce among those from nutritive one. Percentages of resistants to Cc, Em, or Rp, were low, but higher in the nutritive than in the marine medium.

Ampicillin-Resistant Bacteria
579 ampicillin-resistant bacteria were isolated corresponding to 286 and 293 colonies from H1 and H2 areas, respectively. From each zone, 70.63 and 74.06%, respectively, were isolates from the marine medium. Table 1 summarizes the distribution of isolates per sample.  Table 2 shows the order of prevalences. Distributions of resistants to each antibiotic were significantly different among the samples (H1L, H2L, H1S, and H2S) and with respect to the culture media for almost all antibiotics (p < 0.05).

MDR of Ampicillin-Resistant Isolates
A total of 143 different resistance profiles were identified (Supplementary Table S4), with 2.1% of resistants only to Ap, and 2.6% to all the antibiotics. The most frequent profiles included resistance to eight antibiotics (128 isolates and 22 profiles), while the most diverse included resistance to 6 or 7 antibiotics (26 profiles each). 92.7% of the isolates were MRB [resistance to 4 or more antibiotic classes (Narciso-da-Rocha and Manaia, 2016)]. Table 1 summarizes the profiles/isolates ratios in each group.

Phylogenetic Identification of Bacterial Isolates
Partial 16S rRNA sequences of the 345 isolates considered different by their MDR profiles and colony morphology were obtained, but phylogenetic assignation of some of them sometimes gave different results, mostly at the species level, when comparisons were made using either sequences at GenBank or phylogenetic trees (Supplementary Figures S1-S7); however, most affiliations were coincident at the genus level (Supplementary Table S5). From H1 and H2 areas, 168 and 141 isolates, respectively, were considered to be different based on sample origin, multi-resistance profiles and 16S rDNA sequences. These Ribbons connect antibiotic to which resistant isolates have been found and the groups of isolates in which they were found. Ribbons wideness represents percentages of bacteria resistant to the corresponding antibiotic among the isolates of the corresponding considered group.

Group of isolates
Order of prevalence isolates were assigned to phyla Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria and classes Alphaproteobacteria, Bacilli, Gammaproteobacteria, Actinobacteria, Flavobacteriia, Betaproteobacteria, Sphingobacteriia, and Cytophagia ( Figure 4A). Moreover, 48 different genera were identified, 52.1% were ubiquitous, 27.1% specifically of H1, and 20.8% of H2 ( Figure 4B). Alphaproteobacteria showed the greatest diversity in genera and Cytophagia the least. At the genus level, the highest was for Stenotrophomonas (x = 1.000) and the lowest was for Pseudomonas (x = 0.273). All Ochrobactrum, Stenotrophomonas, and Streptomyces were resistant to at least six antibiotics.

Levels of MDR Among Phylogenetically Ascribed Isolates
Sm R was the most frequently observed and Cc R the least observed. The prevalent resistances (Figure 5) for the phylum/class Actinobacteria were Km R and Tm R , and the less common were Em R and Vm R . For Bacteroidetes, Km R and Sm R were the most frequent, while Em R was the least. For Firmicutes, Tm R and Vm R were the most and least frequently found, respectively. For Proteobacteria, the most prevalent was Sm R and the least was Cc R .

Putative Intrinsic or Acquired Resistance of the Isolates
Resistance to a specific antibiotic was considered as putatively intrinsic in a particular phylogenetic group given the high majority (≥90%) of isolates that belong to that group were resistant (Figure 5). Thus, at the phylum level (Figure 5A), Bacteroidetes were intrinsically resistant to Km, Proteobacteria to Sm, and InR was neither observed in Firmicutes down to class nor in Actinobacteria down to genus.
At the class level (Figure 5B), Alphaproteobacteria showed InR only to Sm, and Betaproteobacteria to Em, Nx, and Sm, but for Gammaproteobacteria no InR was detected at this level or down to genus. Class/Order Flavobacteriia/Flavobacteriales showed InR to Km, Nx and Sm and class Sphingobacteriia to Cz and Km.
At the order level (Figure 5C), Caulobacteriales showed InR to Sm, Rhizobiales to Nx, Sm and Tm, and Sphingomonadales to Nx and Sm.
At the family level (Figure 5D), considering only those with more than one genus or species, Bacillaceae isolates appeared intrinsically resistant to Tm; Brucellaceae to Cz, Em, Sm, Tm, and Vm; Caulobacteriaceae to Cz, Nx and Sm; Erythrobacteriaceae to Nx and Sm; and Sphingomonadaceae to Sm.
Five genera of Alphaproteobacteria showed InR: Altererythrobacter to Sm and Tm; Brevundimonas to Cz, Nx, and Sm; Erythrobacter to Nx, Sm, and Tc; Ochrobactrum to Cz, Em, Sm, Tm, and Vm; Ponticaulis only to Sm; all Bacillus isolates were resistant to Tm; and all Mesonia isolates to Sm. No genus showed InR to Cc, Km, or Rp; and all isolates of the genera Altererythrobacter, Bacillus, and Mesonia were susceptible to Cc, Vm, and Em, respectively ( Figure 5E).

Correlations of Resistances Pairs in Each Phylogenetic Group
The kappa correlation index (κ) for pairs of antibiotic resistances was determined for the whole set of isolates and for the phylogenetic groups that included at least eight isolates.
Considering all the isolates, no significant correlation appeared, but specific phylogenetic groups showed several positive and negative correlations. Only some of the positive ones reached κ > 0.400, corresponding to moderate and up to perfect correlations (κ = 0.810-1.000) (Supplementary Table S6). Perfect correlations appeared for Cz R /Tc R in Betaproteobacteria and Cc R /Tc R and Em R /Sm R in Sphingobacteriia. At the genus level, these correlations were found in Mesonia for Cc R /Tc R and Nx R /Tm R , in Microbacterium for Cc R /Rp R , in Paenibacillus for Cc R /Em R , Cc R /Km R , Em R /Nx R , Km R /Nx R , Cz R /Tm R and Em R /Km R , in Pseudoalteromonas for Tm R /Vm R , and, finally, in Pseudomonas for Cc R /Em R and Rp R /Tc R . These correlations configure networks, which were represented for κ values higher than 0.4 (Figure 6). In these networks, some pairs of resistances with moderate or optimal correlations are shared by several genera while others appeared only in specific ones.

Environmental Conditions of the Sampling Zones
The samples used in this work contained heavy metals with a slightly higher level of the most abundant elements in the H1 area and with lower salinity and conductivity when compared to the H2 area. This result may be due to the effect of a higher influence from the Tinto River water, which constitutes the leading heavy metal source of the H1 station. The industrial plants around H2 also dump metals and metalloids into the FIGURE 5 | Distribution among phylogenetic groups of isolates resistant to the tested antibiotics. Graphs correspond to (A) phylum, (B) class, (C) order, (D) family, and (E) genus. Ribbons connect the antibiotic to which resistance was detected with the phylogenetic group to which resistant isolates were assigned. The wideness of each ribbon represents the percentage of isolates within each phylogenetic group resistant to the antibiotic that it connects to. Black ribbons correspond to percentages of isolates ≥90%. water but with a lower contribution than from the rivers (Pérez-López et al., 2011). The concentration of metals in the estuarine water is much lower than what is usually found in the Tinto River, which is probably the result of the dilution effect from seawater entry. Also, the higher pH of the water could trigger the precipitation of some metals (Rodgers et al., 2019). The higher concentration of phosphor in H1S versus H2S may be due to the input of phosphogypsum-rich wastewater from a nearby fertilizer plant (Sainz et al., 2004), which recently was shown to contaminate the estuary with heavy metals through a weathering process (Papaslioti et al., 2018). The presence of other contaminants in the estuary was not the aim of this study; however, it is advisable, based on the discharge of effluents from the Huelva WWTP into the estuary, which can provide contamination by antibiotics and other pharmaceuticals as has been detected in the effluents from urban WWTPs around the world. Considering the differences observed in metals composition and the physicochemical characteristics of both sampling zones, a differential impact of these factors on the biodiversity could not be ruled out.

Abundance of Culturable Bacteria
As customarily expected by the marine origin of most waters arriving into the estuary, the abundance of culturable bacteria was higher on the marine medium than on the nutritive medium. Also, in both phases, the abundance of bacteria was higher in H1 than in H2. Our results are also in line with previous reports that have shown a higher bacterial richness in sediments than in the water column of several environments (Ye et al., 2009;Dai et al., 2016;Narciso-da-Rocha and Manaia, 2016). The lower abundances found in the H2 samples on both media might be due to the possible presence of toxic pollutants emitted by nearby industries. Additionally, H2 is closer to the WWTP, which could discharge antibiotics into the estuary that may not have been completely removed during the sewage treatment as reported in other cases (Zhou et al., 2013;Tran et al., 2016). The difference in abundance between the phases indicates the accumulation of bacteria in the sediments with higher stability in the estuarine system. The sediments offer protection against biotic and abiotic stresses like tidal water movements. Bacteria counts of sediments may even be underestimated as some studies claim (Hassard et al., 2017).

Antibiotic Resistance Levels and Diversity in the Estuary
Eleven antibiotics representatives of different mechanisms of action were chosen. The antibiotic concentrations to use for susceptibility testing of environmental bacteria have not been standardized; and the reported studies used different concentrations (D'Costa et al., 2006;Bhullar et al., 2012;Voolaid et al., 2012;Narciso-da-Rocha and Manaia, 2016;Jardine et al., 2017). In this study, we wanted to focus on the ARs that might be relevant for the ecosystem and for clinical microbiology. So, concentrations related to the reported resistance/susceptibility breakpoints for various clinical bacteria according to EUCAST (European Committee on Antimicrobial Susceptibility Testing [EUCAST], 2013) were used. However, the used values were slightly higher than the breakpoints of EUCAST to avoid overestimations of the resistance level.
Resistants to Km, Nx, Sm, Tc, or Tm were more abundant among bacteria growing on the marine medium than on the nutritive medium for all samples. Furthermore, Cc R , Em R , or Rp R were rarely found in bacteria growing on marine media. These results indicated that the population of bacteria growing on one or the other media are not the same and carry determinants for different ARs. The highest difference between media was for Tc R , common in isolates from the marine medium. Several studies have shown a great abundance of Tc R genes in different environments including marine sediments , river estuary (Lin et al., 2015), or wastewater and activated sludge (Auerbach et al., 2007;Zhang and Zhang, 2011;Yang Y. et al., 2013). Moreover, Tc has been reported to be non-biodegradable in activated sludge or different conditions (Prado et al., 2009;Cetecioglu et al., 2014), so it needs special treatment for degradation in water (Wang et al., 2011) and sediments (Chang and Ren, 2015); these factors could help with the selection of Tc R bacteria.

Prevalence of Resistance to Individual Antibiotics Among Ap R Bacterial Isolates
Resistance to Sm was the most prevalent considering all isolates (H), those from liquid (HL) or sediment (HS) phases or H1 and H2 separately. Tc R was prevalent in most of the groups of isolates obtained on the marine medium except isolates from H2S-M and HS-M, which had a slightly higher prevalence of Nx R . Likewise, Tm R was prevalent on those on the nutritive medium except for the isolates from H1S-N and HS-N for which Cz R was more prevalent. These data show that, in general, Sm R , Tc R , and Tm R were the prevalent resistances among the isolates mostly in agreement with data from the direct plating of original samples. Thus, the isolates obtained are a good representation of the estuarine culturable ARB. The quantitative differential presence of each bacterium in the original samples was partially deconvoluted when colonies were selected in the isolation step and this might account for the differences observed.
Some studies have shown that streptomycin can persist in aquatic environments for some time and its degradation products still maintain certain antibacterial activity (Shen et al., 2017). On the other hand, selection of Sm R bacteria in in vitro evolution experiments with low concentrations of antibiotic has been observed even after only 96 h (Spagnolo et al., 2015). Moreover, some studies have shown a relationship between heavy metal and Sm resistances (Matyar et al., 2014) or aminoglycoside resistance genes (Zhao et al., 2019) and plasmids have been found to contain determinants for co-resistance to Sm and Cu (Sundin and Bender, 1996). So, one possible reason for the high prevalence of Sm R in the estuary's isolates could be the co-selection due to the presence of Cu or other heavy metals. Also, correlations of Tc R genes prevalences with the presence of heavy metals, also present in the Huelva estuary, such as Cu, Mn, Ni, and Zn have been reported (Knapp et al., 2017). Recently, in a paleontological study of sediments contaminated with Zn, bacteria were found to share Tm R with Zn resistance, which was related to Zn concentration (Dickinson et al., 2019). The co-selection for Zn R likely might be one of the causes for the prevalence of Tm R in bacteria from the Huelva estuary. In addition to possible effects by other contaminants, different studies have shown the coexistence, and sometimes a correlation, of AR and heavy metals resistances in waters and sediments [reviewed by Nguyen et al. (2019)].
Further molecular studies on the genetics determinants for Sm R , Tc R , and Tm R in the isolates will give insights on the causes of such high prevalence and mechanisms of co-selection.

MDR Levels and Profiles Diversity Among Ap R Bacterial Isolates
Multi-drug resistance levels were calculated for different groups of isolates using the MAR index median as the comparative parameter. The diversity of multi-resistance profiles were assessed as the ratio number-of-profiles/number-of-isolates. These parameters could give insights on the availability of resistance's determinants, some of which might be transmissible, as well as on the diversity of putative groups of them that might be transmitted together.
Based on the diversity of colony morphologies and sample origin, 579 Ap R -bacteria were isolated to determine the level and profiles of MDR. This is one of the highest numbers of isolates analyzed in similar studies using culturable bacteria. The test of resistance to the other 10 antibiotics carried out on this collection showed diverse resistance profiles. 92.7% of the whole set of isolates were resistant to 4 or more classes of antibiotics and the median of MAR indexes corresponded to resistance to 7/11 of the tested antibiotics. These results represent outstanding levels of multi-resistance and as far as we know, one of the highest measured up to now in a non-clinical setting.
Higher MDR levels were found in isolates from sediments than from those in waters when considering all isolates from the estuary (H), or each individual zone (H1 or H2). This difference between phases could be due to the lower mobility of sediments and the possible accumulation of antibiotics, metals, or other pollutants (including ARGs) that may facilitate the selection of resistant bacterial populations, as has been proposed for other systems (Rodgers et al., 2019). Furthermore, the accumulation of resistant bacteria in sediments could facilitate bacterial cell contact; and, therefore, to increase the efficiency of HGT (Zhang et al., 2018).
The diversity of MDR profiles is higher in sediments than in waters when considering all isolates (H) or those from H1 zone alone. When isolates were segregated by the culture medium used for isolation, only those on the marine medium showed higher profile diversity in sediments than in waters.
When comparing whole sets of isolates from the two sampling zones (H1 and H2), no difference was found in the MDR level. However, when isolates were segregated by the culture medium a significant difference appeared between the isolates on the nutritive medium, with the highest value for those of H2.
Thus, in general, the isolates on the nutritive medium from sediments were resistant to more antibiotics and showed less diverse profiles than those from water samples. However, those obtained on the marine medium showed higher levels of multiresistance and also higher profiles diversity when considering the isolates obtained from sediments versus those isolated from water samples.
To try to explain the observed differences between the levels and diversity of MDR for the isolates from sediments of H1 and H2 zones on the nutritive medium, our primary hypothesis was the more direct influence of the effluents of the WWTP on H2. This might be responsible for the higher MDR in the isolates on the nutritive medium since the effects of WWTP effluents on bacterial resistance have been shown in several studies (Guo et al., 2017;Manaia et al., 2018). Nevertheless, we noticed that a number of isolates from the sample H2S-N had the same MDR profile and had similar colony morphology. They were initially selected because they appeared in high numbers in the original selection step. When only one of these isolates was considered, the MAR median for the H2S-N sample decreased to 0.636, with no difference with H1S-N. The ratio profiles/isolates became in this way 64.7 and 46.2% for H2S-N and H2S, respectively; and were more similar to the behavior of H1 isolates. Hence, the overrepresentation of a particular bacterium with a high MAR value was biasing the MAR value of the corresponding sample. This bias stresses the need for deconvolution of the isolates to detect the real levels of resistance and diversity of MDR. We thus concluded that not much difference existed between the MDR levels of isolates on the nutritive medium from the sediments of both sampling zones.

Biodiversity of Culturable MRB From the Estuary
16S rDNA sequences from 345 isolates were firstly compared with those available in GenBank database and assigned to eight phylogenetic classes found in both sampling areas in similar proportions. Alphaproteobacteria constituted around 50% of the isolates in each zone. This percentage suggests that the tidal movement of water in the estuary could play an essential role in the homogenization of the bacterial composition at both sampling sites. However, analysis at the genus level showed higher diversity differences and stressed the need to phylogenetically ascribe down to the lowest possible phylogenetic level to obtain a more accurate picture of diversity distribution differences. A higher difference in the classes' distribution was observed between liquid and sediment phases: Actinobacteria and Betaproteobacteria were more abundant in the sediments; and Gammaproteobacteria and Flavobacteriia in waters.
A more detailed phylogenetic analysis using trees built with sequences at RDP identified 48 genera. Although most of the genera were common to both sampling zones, some were found only in one zone, which perhaps reflects physicochemical differences between them.
The phylogenetic composition of the culturable Ap R isolates found in this study is similar to the diversity of bacterioplankton detected in other estuaries or marine environments when compared at the phylum and class levels (Selje et al., 2005;Jamieson et al., 2012;Du et al., 2013;Sjöstedt et al., 2014;Vaz-Moreira et al., 2014). However, a deeper level of phylogenetic assignation is needed in most studies to be able to make further comparisons.
The majority of the identified species have not been described as pathogens. However, some of them may be considered as opportunistic pathogens or have been found in some sporadic cases of infections (Supplementary Table S4). So, some health concerns should be considered for the bacteria in the estuary due to their high level of multi-resistance. Since the phylogenetic diversity of MDR bacteria found in the estuary is high, the possible dissemination of AR or MDR to diverse other bacteria is also likely.

Intrinsic and Acquired Resistance
Intrinsic resistance is the most common form of resistance in some bacteria (Cox and Wright, 2013). In this study, the putative InR of four phyla, seven classes, and 11 genera represented by eight or more isolates were analyzed and putative InR was found in Bacteroidetes and Proteobacteria. Different results have been found in other studies (Walsh and Duffy, 2013). Hence, a selection of bacteria that depends on the presence of different substances and physicochemical conditions should have occurred in the Huelva estuary.
At the class level, Gammaproteobacteria did not show any InR, while Alphaproteobacteria, Betaproteobacteria, and Flavobacteriia showed InR to Sm. These last two also shared InR to Nx. On the other hand, Flavobacteriia and Sphingobacteriia shared InR to Km. Moreover InR to Em was only found in Betaproteobacteria and Cz in Sphingobacteriia.
Among the genera that did not show any InR, the results for Pseudomonas contrast with the case of clinical P. aeruginosa, which is intrinsically resistant to several antibiotics (Alvarez-Ortega et al., 2011;Breidenstein et al., 2011;Vaz-Moreira et al., 2012). Prevalence of antibiotic resistance in environmental Pseudomonas is related to the phylogenetic species level rather than to the genus level (Vaz-Moreira et al., 2012), and for environmental P. aeruginosa is lower than for clinical isolates (Liew et al., 2019). In a recent study (Luczkiewicz et al., 2015), Pseudomonas spp. were isolated from influents, effluents, and receiving marine-water from a WWTP, recovering fewer isolates from marine water than from the water from the WWTP. Furthermore, P. putida was the predominant species in all samples whereas P. aeruginosa was the second major species in the WWTP effluents (enriched from the influents water), but the latter was found in a much lower proportion in receiving water. This effect could explain the lack of this species among our isolates. Additionally, in the abovereferred report, P. putida and P. aeruginosa isolates were tested for AR to 14 antibiotics, but no InR would have been found under the criterion used by us. Our isolates of the remaining analyzed genera showed putative InR to one; as Bacillus, Mesonia, and Ponticaulis; or several antibiotics. The genus with the highest numbers was Ochrobactrum that had five InRs.
None of the phylogenetic groups showed InR to Cc or Rp, and the number of isolates resistant to these antibiotics was generally low. The percentage of Cc R -isolates reached up to 50% only for Pseudomonas. These low levels of resistants may be a consequence of the limited use of Cc because of its toxicity (Ingebrigtsen et al., 2017); and, therefore, the existence of low selective pressure for this antibiotic in the environment. It has been shown that Rp R in pathogenic bacteria has a high metabolic cost, which significantly affects the fitness of the bacteria making their ARGs uncommon (Hall et al., 2011). The fitness cost might have conditioned most isolates, but more than 80% of the Ochrobactrum isolates showed Rp R and seemed to be able to work around these fitness costs.
Non-intrinsic resistances should be considered as acquired independently of the mechanism by which they were obtained (mutation or HGT).

Correlation of the Presence of Pairs of Antibiotic Acquired Resistances in Phylogenetic Groups of Isolates
In not considering the putative InRs discussed above, the correlations of resistances to pairs of antibiotics were low when all available isolates were considered. These data indicated a low probability of cross or co-resistance to any of the tested pairs of antibiotics in the whole bacterial population. However, some moderate to almost perfect correlations were found at the class level, more frequently in Beta or Gammaproteobacteria, and strong ones at the genus level.
Strong positive correlations were observed at the genus level, particularly in Paenibacillus. In the network built based on κ >0.4 values, Paenibacillus appears central sharing many correlated pairs of resistances with other genera suggesting possible sharing of MGEs. Some of these elements have been described in this genus (Pednekar et al., 2010;Soundarapandian et al., 2013). The presence of these elements indicates that this genus might act as an ARG exchanger among different genera, and perhaps, plays a crucial role in the resistance dissemination in the Huelva estuary. Apart from this, an accumulation of ARGs in this genus has been described but classified as intrinsic (Pawlowski et al., 2016). Our hypothesis about Paenibacillus' role in the ecosystem deserves to be studied in the future because specific published data about this genus suggest a great versatility among its members. First, this genus contains a high number of species and new ones are described frequently (Sáez-Niero et al., 2017), so this genus displays a high intragenus genetic variability. In the second place, Paenibacillus species have been described in a variety of environments, some of them interacting with other organisms (Finkelshtein et al., 2015;Grady et al., 2016;Xie et al., 2016), which suggests their potential for affecting ecology. Also, strains of the same species showed numerous variations in genome sequences with a high rate of recombination, which makes them prone to exchange DNA fragments by horizontal genetic exchange Xu et al., 2017). Likewise, several species of the genus produce diverse antibiotic substances Pasari et al., 2019), which could make them players in the shaping of their ecosystems. Furthermore, Paenibacillus genomes have a plethora of genes for adaptation, some of which considered to be acquired by HGT Xu et al., 2017). Moreover, some studies have considered that the species of the genus evolve not as individual species, but inter-related in a pan-genomic evolution model (Xu et al., 2017). Finally, the resistance elements found in genomes of some Paenibacillus species could be the origin of resistance elements found in other genera (Guardabassi et al., 2005;Liu et al., 2016); and, thus, they show mobilization potential. All this suggests that this genus may play a crucial role in any biological system because of its capabilities to adapt and interact with other organisms. Considering that some Paenibacillus species can infect humans, animals, and plants and are quite ubiquitous, they could be an essential and currently unconsidered threat when present in an ecosystem, either directly or by its putative ARGs exchanger role, which could affect the ecology of the systems where they are present. A recent study by Wright's group (Pawlowski et al., 2018) also suggested a role for Paenibacillaceae bacteria in environmental mobilization and dissemination of resistance genes.

CONCLUSION
In conclusion, the culturable bacterial population from the Huelva estuary contains highly-multi-resistant bacteria some of which have been isolated and characterized for MDR profiles to 11 antibiotics, which show a diversity of combinations of resistances. The analysis of AR of the phylogenetically ascribed isolates has shown the existence of putative InR in most of the genera. However, some other ARs may be acquired and some are perhaps transferrable. Moderated to optimal correlations of pairs of resistances in certain phylogenetic groups suggests that these resistances could be transferred among their members. In particular, Paenibacillus showed a large number of these correlations, many of which are shared by other genera as a network study has shown. From these results, a hypothesis for future studies can be outlined in which Paenibacillus could have a role as a player in the dissemination of resistance in this environment.

DATA AVAILABILITY STATEMENT
The 16S rRNA sequences generated for this study can be found in the GenBank (accession numbers LT601033-LT601378).

AUTHOR CONTRIBUTIONS
JA contributed to the conception and design of the study, sampling, and wrote the manuscript. BE-C, HM-F, and JA performed the experiments and analyzed the data. BE-C and JA performed the statistical analysis. All authors contributed to the manuscript revision, read and approved the submitted version.

ACKNOWLEDGMENTS
We want to thank Dr. J. J de la Cruz Troca from the Department of Preventive Medicine and Public Health and Microbiology of the Universidad Autónoma de Madrid for help in using statistical software. BE-C was a doctorate fellow supported by the Spanish Agency for International Development and Cooperation (AECID) and the work presented here is in part from his Ph.D. thesis (Eduardo-Correia, 2016).