Antibiotic Resistance Genes and Potentially Pathogenic Bacteria in the Central Adriatic Sea: Are They Connected to Urban Wastewater Inputs?

: Despite last decades’ interventions within local and communitarian programs, the Mediterranean Sea still receives poorly treated urban wastewater (sewage). Wastewater treatment plants (WWTPs) performing primary sewage treatments have poor efﬁciency in removing microbial pollutants, including fecal indicator bacteria, pathogens, and mobile genetic elements conferring resistance to antimicrobials. Using a combination of molecular tools, we investigated four urban WWTPs (i.e., two performing only mechanical treatments and two performing a subsequent conventional secondary treatment by activated sludge) as continuous sources of microbial pollution for marine coastal waters. Sewage that underwent only primary treatments was characterized by a higher content of traditional and alternative fecal indicator bacteria, as well as potentially pathogenic bacteria (especially Acinetobacter , Coxiella , Prevotella , Streptococcus , Pseudomonas , Vibrio , Empedobacter , Paracoccus , and Leptotrichia ), than those subjected to secondary treatment. However, seawater samples collected next to the discharging points of all the WWTPs investigated here revealed a marked fecal signature, despite signiﬁcantly lower values in the presence of secondary treatment of the sewage. WWTPs in this study represented continuous sources of antibiotic resistance genes (ARGs) erm B, qnr S, sul 2, tet A, and bla TEM (the latter only for three WWTPs out of four). Still, no clear effects of the two depuration strategies investigated here were detected. Some marine samples were identiﬁed as positive to the colistin-resistance gene mcr -1, an ARG that threatens colistin antibiotics’ clinical utility in treating infections with multidrug-resistant bacteria. This study provides evidence that the use of sole primary treatments in urban wastewater management results in pronounced inputs of microbial pollution into marine coastal waters. At the same time, the use of conventional treatments does not fully eliminate ARGs in treated wastewater. The complementary use of molecular techniques could successfully improve the evaluation of the depuration efﬁciency and help develop novel solutions for the treatment of urban wastewater.


Introduction
Treated and untreated urban wastewaters (i.e., sewage) represent important point sources of microbial pollution for marine coastal areas [1][2][3]. Despite last decades' interventions within regional and communitarian programs, the majority of the objectives addressing the management of urban wastewater have not been achieved or partially achieved, with the Mediterranean Sea still receiving poorly treated sewage. Discharge of untreated or minimally processed wastewater is still a common practice in certain Mediterranean countries, while other regions face chronic malfunctioning of a few wastewater treatment plants (WWTPs) [4,5]. Basic mechanical treatments (i.e., "primary treatments") aim at removing solids, fats, and sand. Still, they do not substantially affect fecal pollution indicators (e.g., fecal coliforms and streptococci, fecal sterols, and potentially pathogenic bacteria-i.e., PPB) nor bacterial community structure [6]. Most WWTPs are designed to remove nutrients and biodegradable organic compounds and reduce the pathogen load through additional aerobic biological processes (i.e., "secondary treatment") and a final disinfection step. However, conventional treatments do not target most emergent contaminants, including antibiotic-resistant genes (ARGs), antibiotic-resistance bacteria, and opportunistic pathogens, with urban WWTPs becoming a major source of microbial pollution [7][8][9][10]. Marine coastal environments, as the WWTP effluent receiving bodies, have the potential of disseminating and promoting the establishment of ARGs and other microbial pollutants [2,11,12], thus facilitating the emerging of new antimicrobial resistance (AMR) pathways [13][14][15] and the exposure of millions of people to disease and opportunistic pathogens [16][17][18]. Recent studies have shown that people's probability of contracting gastrointestinal and non-gastrointestinal illnesses increases in correlation with aquatic recreational activities because of antibiotic-resistant bacteria and other community-acquired microorganisms [19,20]. Furthermore, frequent bathing in coastal waters has been shown to be associated with gut colonization by E. coli harboring ARGs, including, for example, bla CTX-M genes, which provide resistance to clinically relevant antibiotics (cefotaxime) and which are easily mobilized via horizontal gene transfer (HGT), thus increasingly found in natural environments as well as in WWTPs [21][22][23].
The microbiological quality of coastal water is routinely assessed by enumerating Fecal Indicator Bacteria (FIB) (i.e., Escherichia coli and enterococci) through culture-based methods [24,25]. However, FIB association with gastrointestinal and non-gastrointestinal diseases is currently questioned. Several bacteria occurring in natural environments and causing infectious disease are not primarily associated with human or animal fecal contamination [26][27][28]. At the same time, studies have shown that substantial populations of fecal bacteria can persist and regrow in aquatic environments, followed by the discovery of environmentally adapted populations of E. coli and by retrieving fecal bacteria in the absence of obvious fecal sources [29][30][31]. Moreover, most bacteria are or become unculturable because of specific metabolic requirements, physiological conditions, or the disinfection process [7].
It is also worth pointing out that routine monitoring of coastal waters is specifically addressed to the assessment of the quality of bathing waters, whereas the control of the environmental status in areas of treated wastewater discharge is spatially scattered and mainly dependent on local/regional policies. In this respect, information about the diffusion of ARGs and of non-conventional, potentially pathogenic bacteria in the Mediterranean is very scarce. To the best of our knowledge, no dedicated study on the presence of ARGs in Adriatic waters has been performed so far. Driven by this paucity of data, we aimed at depicting the transfer of these types of microbial pollution from WWTPs to seawater, following the hypothesis that the sewage treatment technology affects not only the spreading of classical FIB but also the abundance of other potential pathogens and ARGs. We, therefore, sampled both treated wastewater and seawater in the proximity of discharge points of four WWTPs characterized by different treatment technologies.

Study Area and Sampling
Four urban WWTPs with outfalls in the central Adriatic Sea were here investigated ( Figure 1) in the framework of the Interreg Italy-Croatia project AdSWiM [32]. The Katalinića Brig and Stobreč WWTPs are part of the Split-Solin public sewage system (Croatia). Their submarine outfalls discharge in the Split Channel, at 1.30 km (43 m depth) and 2.75 km (37 m depth) from the shore, respectively, and are about 5 km away from each other [33]. Katalinića Brig WWTP receives both municipal and rainfall-runoff wastewater (122,000 PE capacity, about 35,000 m 3 /day average flow rate), and it is designed to remove solids by coarse and fine grids (i.e., primary treatment). Stobreč WWTP receives mostly Water 2021, 13, 3335 3 of 20 municipal wastewater (138,000 PE capacity, about 33,000 m 3 /day average flow rate). Like the Katalinića plant, it is designed for the mechanical treatment of sewage, but oil is also removed. Zadar WWTP receives mainly municipal wastewater with minor contributions of other types of wastewaters (i.e., rainfall-runoff and about 5% from septic tank wastewater). It is designed with 100,000 PE capacity, with an average flow rate of 22,600 m 3 /day. The outfall pipe discharges at 2 km from the shoreline (34 m depth). In the Zadar plant, after removing coarse and fine solids, sand, and grease by mechanical treatment (i.e., by grids and gravitational deposition), wastewater is treated in activated sludge tanks for the aerobic biodegradation of organic pollutants and the removal of nutrients (i.e., secondary treatment). Before the unloading, the treated sewage is clarified by gravity deposition and transition in secondary precipitators. The Francavilla a Mare WWTP serves the municipality of Pescara (Italy) for the treatment of municipal wastewater (38,000 PE). The outfall pipe discharges at 2.5 km from the shoreline (15 m depth). As with the Zadar WWTP, Francavilla WWTP includes a biological treatment by activated sludge. However, the sewage inflow is divided into two streams, of which only one is treated in the activated sludge tanks due to plant under-sizing. The other stream undergoes a series of mechanical treatments, including surface aeration in a trickling filter, that decrease the biological oxygen demand (BOD 5 ). The two streams are then reunited for final disinfection by peracetic acid. km (37 m depth) from the shore, respectively, and are about 5 km away from each other [33]. Katalinića Brig WWTP receives both municipal and rainfall-runoff wastewater (122,000 PE capacity, about 35,000 m 3 /day average flow rate), and it is designed to remove solids by coarse and fine grids (i.e., primary treatment). Stobreč WWTP receives mostly municipal wastewater (138,000 PE capacity, about 33,000 m 3 /day average flow rate). Like the Katalinića plant, it is designed for the mechanical treatment of sewage, but oil is also removed. Zadar WWTP receives mainly municipal wastewater with minor contributions of other types of wastewaters (i.e., rainfall-runoff and about 5% from septic tank wastewater). It is designed with 100,000 PE capacity, with an average flow rate of 22,600 m 3 /day. The outfall pipe discharges at 2 km from the shoreline (34 m depth). In the Zadar plant, after removing coarse and fine solids, sand, and grease by mechanical treatment (i.e., by grids and gravitational deposition), wastewater is treated in activated sludge tanks for the aerobic biodegradation of organic pollutants and the removal of nutrients (i.e., secondary treatment). Before the unloading, the treated sewage is clarified by gravity deposition and transition in secondary precipitators. The Francavilla a Mare WWTP serves the municipality of Pescara (Italy) for the treatment of municipal wastewater (38,000 PE). The outfall pipe discharges at 2.5 km from the shoreline (15 m depth). As with the Zadar WWTP, Francavilla WWTP includes a biological treatment by activated sludge. However, the sewage inflow is divided into two streams, of which only one is treated in the activated sludge tanks due to plant under-sizing. The other stream undergoes a series of mechanical treatments, including surface aeration in a trickling filter, that decrease the biological oxygen demand (BOD5). The two streams are then reunited for final disinfection by peracetic acid. The sampling activities were carried out (almost) monthly in spring and summer 2020 (from April to September). We collected extra samples in the Croatian sites in July 2019 and autumn 2019/winter 2020 (Table S1). Samples of treated sewage were taken before unloading, just upstream the injection into the discharging pipeline. Seawater samples were collected in the proximity of each WWTP outfall 1 m above the main diffusion point of the pipelines by means of 5 L Niskin bottles. While samples for FIB enumeration were taken with the same frequency for seawater and treated sewage, seawater samples The sampling activities were carried out (almost) monthly in spring and summer 2020 (from April to September). We collected extra samples in the Croatian sites in July 2019 and autumn 2019/winter 2020 (Table S1). Samples of treated sewage were taken before unloading, just upstream the injection into the discharging pipeline. Seawater samples were collected in the proximity of each WWTP outfall 1 m above the main diffusion point of the pipelines by means of 5 L Niskin bottles. While samples for FIB enumeration were taken with the same frequency for seawater and treated sewage, seawater samples dedicated to DNA extraction for the subsequent analyses of ARGs and microbial community structure were taken at selected time points.

Microbiological Analyses
Sample aliquots of 50 mL, 5 mL, 0.5 mL, and 0.05 mL were filtered, in three replicates, through 0.45 µm pore size sterile membrane filters (47 mm diameter made in hydrophilic mixed cellulose esters, PALL Corporation) and incubated onto selective media for the enumeration of fecal indicator bacteria (FIB, i.e., E. coli and enterococci), in the following conditions: at 44 • C for 24 h onto TBX (Tryptone Bile X-glucuronide) agar plates, and at 37 • C for 48 h onto on Slanetz-Bartley agar plates, for E. coli and enterococci, respectively. Colonies that displayed the characteristic positive colors (i.e., blue-green for TBX-agar medium and reddish-brown for Slanetz-Bartley agar medium) were counted. The result was expressed as the number of colonies forming units (CFU) in 100 mL of sample. In this study, Pseudomonas aeruginosa was selected as a model of emerging pathogens of non-fecal origins and environmental relevance [34][35][36]. Its enumeration was performed by incubating membranes onto Cetrimide-agar plates at 44 • C for 48 h, following the procedure used for fecal indicator bacteria. Analyses were performed in two replicates.

DNA Extraction
Seawater and sewage samples were filtered onto 0.22 µm hydrophilic PES membrane filters (PALL Laboratory) and stored at −80 • C until processing. Filtration was endured until clogging of membrane pores. Filtered volumes were in the range of 1-2 L for seawater samples and 40-500 mL for sewage samples, respectively.
Genomic DNA was isolated from filters using the DNeasy PowerWater Kit (Qiagen) following the manufacturer's protocol, with few modifications aimed at increasing the extraction yield [37]. Filters were cut in half and were extracted separately and pooled together after completion of the extraction procedure. After adding the lysing buffer provided by the manufacturer, samples underwent three cycles of incubation at 70 • C for 5 min followed by 2 min vortexing at the maximum speed. We then proceeded as described in the kit handbook. The DNA concentration was measured using a Qubit Fluorometer (ThermoFisher). DNA extracts were dispensed into aliquots and stored at −80 • C until further processing.

Microbial Community Characterization by High-Throughput Amplicon Sequencing
The V4-V5 hypervariable region of the 16S rRNA gene was amplified using universal primer pair 515F-Y/926R, which allows differentiating taxa unresolvable with the 515F/806R primer pair for the V4 region (Parada et al., 2016). Libraries were prepared following the Illumina protocol [38] and then pooled in equal molar ratios for pair-end sequencing (250 bp) on an Illumina MiSeq System (Illumina, San Diego, CA, USA) at the ARGO Open Lab Platform (Area Science Park, Trieste, Italy). Samples were sequenced in two runs. Since the library size of the second run was up to 10× higher than the first one, we subsampled the second run reads prior to any bioinformatics analysis (i.e., we filtered reads with the sum of extended errors ≤ 1 out and then evened the library sizes to 250,000 using "ShortRead" package for random sampling [39]). Sequences were deposited in the NCBI Sequence Read Archive (SRA) with bioproject numbers PRJNA771212 and PRJNA771227. Reads processing was performed using the "DADA2" package (v 1.20.0) in R 4.1.0 [40,41], following the package instructions. After primer removal and visual inspection of the read quality profiles, forward and reverse reads were truncated at positions 220 and 190, respectively. Before running the dada algorithm (pseudo-pooling method), we enforced monocity of the estimated error rates (to better cope with binned scored quality and, thus, preventing false positives) by manual correction. Chimeric sequences were identified and removed in consensus mode (chimeras were calculated to make up about 6.3% of the merged amplicon sequence variants-i.e., ASVs). Non-target-length sequences and any singletons that arose after reads merging and chimera removal were filtered out. Taxonomy was assigned using the native implementation of the naive Bayesian classifier method [42] provided by the DADA2 package and using SILVA SSU reference database version 138 [43,44]. The subsequent downstream analyses (i.e., inspection, visualization, and pre-processing) were performed using "phyloseq", "ggplot2", and "ranacapa" packages [45][46][47]. The whole prokaryotic diversity was captured by the sequencing effort, as shown by the rarefaction curves in Supplementary Materials ( Figure S1). ASVs with unassigned domains and those taxonomically assigned as chloroplast, mitochondria, and Eukaryota were removed. Phyla with both extremely low prevalence and read count were filtered out in preliminary inspection phases. Subsequently, processed sequences were imported into QIIME2 [48], aligned using MAFFT [49], and then used to infer an unrooted phylogenetic tree based on the FastTree [50].

Quantification of Antibiotic Resistance Genes (ARGs)
DNA extracts were diluted in nuclease-free water (i.e., 10× in the case of seawater samples, 10-100× for treated sewage samples, based on total DNA quantification) and analyzed for the quantification of 16S rRNA gene (as housekeeping gene) and 8 ARGs genes conferring resistance to commonly used antibiotic families, such as tetracyclines (tetA), sulfonamides (sul2), macrolides (ermB), fluoroquinolones (qnrS), beta-lactams (bla TEM , bla CTX-M , and bla OXA ), and colistin resistance gene (mcr-1). qPCR assays were performed in SYBR green chemistry (SsoAdvanced universal SYBR Green supermix, Bio-Rad), using the CFX96 Touch Real-Time PCR Detection System (Bio-Rad) in the analytic conditions previously described [51]. Before quantification, a potential PCR inhibition effect was tested by dilution method (i.e., 10×, 100×, 1000×), as described by Di Cesare et al. [52]. Synthetic consensus sequences from the Comprehensive Antibiotic Resistance Database (CARD; [53]) cloned into the pEX-A128 plasmid vector (Eurofins) were used as positive controls. Standard calibrations curves were built using the purified, quantified, and diluted amplicons of each positive control, thus with linearized PCR products to avoid issues due to supercoiled circular configuration [54]. Details about the primer pairs used in this study are shown in the Supplementary Materials (Table S2) with the respective annealing T [51,52,[55][56][57][58][59][60][61]. Quantification was performed in two technical replicates. The specificity of the qPCR assay was investigated by melting curve test. The limit of quantification (LOQ) of the standard curves, as defined by [62], and the other parameters for gene quantification are shown in Supplementary Materials (Table S3). For the samples in which the threshold cycle was below the LOQ but above the limit of the qPCR (theoretically three copies per PCR; Bustin et al., 2009), the analyzed gene was considered present but not quantifiable [63]. ARG abundances were expressed as ARG copies per copy of 16S rRNA gene to prevent differences among sample cell abundances from becoming a confounding factor.

Statistical Analyses
Statistical analyses were carried out in R 4.1.0 [41]. The correlation between the two FIB (i.e., E. coli and enterococci) and between FIB and P. aeruginosa has been tested in both types of samples by Spearman rank correlation analysis due to the unbalanced design, heteroscedasticity of data, and presence of potential outliers, as assessed by data exploration [64].
As regards the 16S sequencing data, we applied a preliminary filtration to remove rare ASVs (total sum > 1 × 10 −4 %) prior to performing any statistical analysis or ordination. The resulting dataset comprised 8,693,932 reads from the 45 analyzed samples, for a total of 14,454 prokaryotic ASVs. Data were not rarefied as the number of reads in samples did not vary by much, while rarefaction could result in information loss and a high rate of false positives [65]. A natural log(1 + x) transformation was applied prior to data ordination of either overall microbiomes or data subsets, using Bray-Curtis distance to calculate the dissimilarity matrices. Taxa that most contributed to beta-diversity were identified by pairwise comparisons (Bray-Curtis distances) of categorical variables of interest by SIMPER analysis, followed by Kruskal-Wallis one-way ANOVA on ranks to test for significance with p-values correction by false discovery rate [66].
As the microbial signature of fecal pollution, a data subset of fecal indicator taxa was obtained extracting ASVs belonging to the traditional fecal indicator taxa (i.e., the family Enterobacteriaceae, which includes the genus Escherichia, and the family Enterobacteriaceae, that includes the genus Enterococcus), and to alternative fecal indicator taxa, according to the approach proposed by Newton et al. [67]. In particular, five feces-associated bac-terial families (i.e., Bacteroidaceae, Porphyromonadaceae, Clostridiaceae, Lachnospiraceae, and Ruminococcaceae) were used as signatures of fecal (both human and non-human) contamination. Three sewer infrastructure-associated bacterial genera (i.e., Acinetobacter, Arcobacter, and Trichococcus) were used as signatures of sewage contamination [31]. The diversity and abundance of potentially pathogenic bacteria were investigated by creating a data subset of 60 genera, as selected from a combination of public databases ( [68]; https: //www.bode-science-center.com/center/relevant-pathogens-from-a-z.html (accessed on 15 September 2021)). For this study, we included Shewanella, as strains affiliated with this genus are considered emerging marine pathogens [28,69].
Differences between sampling areas and between types of treatment (i.e., with or without secondary treatment) for the composition of fecal bacterial communities or in the distribution of PPB were assessed by ANOSIM (analysis of similarity) based on a Bray-Curtis dissimilarity matrix. Due to the high sparsity of PPB in the seawater sample dataset, we performed a differential abundance analysis of those genera using a generalized linear model based on the negative binomial distribution, with p-values corrected by Benjamini-Hochberg corrected p-values using the R package DESeq2 [70].
Statistical differences among groups were examined by multiple comparisons with the Kruskal-Wallis test, followed by a post hoc test using the criterium Fisher's least significant difference, as in the R packages agricolae v1.3.5.
P. aeruginosa counts in the two Split WWTPs showed a peak in April 2020 (1.1 × 10 8 and 0.9 × 10 8 CFU/100 mL), with average values of 2.2 and 2.5 × 10 7 CFU/100 mL and min values of 1.8 × 10 5 and 9.0 × 10 4 CFU/100 mL, for Katalinića and Stobreč plants, respectively. Like FIB, P. aeruginosa in the treated sewage from Zadar WWTP was 10 2 -fold lower than the two Split WWTPs. At the same time, it was often undetected in the treated sewage from Francavilla WWTP (max value: 1.1 × 10 1 CFU/100 mL).
A positive correlation was found between the abundance of E. coli and enterococci both in treated sewage samples (S = 580, ρ = 0.857, p < 0.001) and in seawater (S = 258, ρ = 0.936, p < 0.001). In contrast, there was no significant correlation between the abundance of FIB and Pseudomonas in both seawater and treated sewage. P. aeruginosa counts in the two Split WWTPs showed a peak in April 2020 (1.1 × 10 8 and 0.9 × 10 8 CFU/100 mL), with average values of 2.2 and 2.5 × 10 7 CFU/100 mL and min values of 1.8 × 10 5 and 9.0 × 10 4 CFU/100 mL, for Katalinića and Stobreč plants, respectively. Like FIB, P. aeruginosa in the treated sewage from Zadar WWTP was 10 2 -fold lower than the two Split WWTPs. At the same time, it was often undetected in the treated sewage from Francavilla WWTP (max value: 1.1 × 10 1 CFU/100 mL).
A positive correlation was found between the abundance of E. coli and enterococci both in treated sewage samples (S = 580, ϱ = 0.857, p < 0.001) and in seawater (S = 258, ϱ =

Bacterial Community Structure
Bacterial community structure in treated sewage was mainly dominated by Proteobacteria (29.0-75.9%), Firmicutes (1.4-49.2%), Bacteroidota (2.9-40.4%), and Actinobacteriota (1.2-20.8%), as shown in Figure S2. Ordination based on Bray-Curtis distances (i.e., PCoA, Figure S3) showed that treated sewage samples obtained from Katalinića and Stobreč WWTPs (i.e., primary treatment only) were very similar in bacterial community structure and different from those obtained by Francavilla and Zadar WWTPs, which formed two separated clusters. In the presence of secondary treatment (i.e., Francavilla and Zadar WWTPs) Firmicutes decreased significantly (SIMPER, fdr corrected p-value < 0.001) by 24.7% on average, in favor of Actinobacteriota, Chloroflexi, Patescibacteria and Planctomycetota (SIMPER, fdr corrected p-value < 0.01 for Actinobacteriota, p-value < 0.001 for the others). Francavilla WWTP effluent had an average percent contribution of Actinobacteriota lower than the Zadar WWTP one by 10.1% (SIMPER, fdr corrected p-value < 0.01). Deinococcata, Acidobacteriota, and Myxococcota contribute to minor (<5%), yet significant, differences between the two WWTPs. Abundances of ASVs affiliated with traditional and alternative fecal indicator bacteria in treated sewage are shown in Figure 3. Traditional fecal indicator bacteria in treated sewage from WWTPs with primary treatment only (i.e., Katalinića and Stobreč) accounted for relative abundances of 0.6-3.6%, while in the presence of a secondary treatment (i.e., Francavilla and Zadar), they were statistically lower (ANOSIM, R = 0.5852, p-value < 0.001) and ranged between 0.03-1.6%. In both cases, Enterobacteriaceae exceeded Enterococcaceae by about an order of magnitude. Alternative feces-associated families (i.e., Bacteroidaceae, Porphyromonadaceae, Clostridiaceae, Lachnospiraceae, and Ruminococcaceae) covered a very high proportion of the total diversity (up to 34.5%, in Katalinića, treated sewage), although Clostridiaceae and Porphyromonadaceae were detected in low abundances (<0.5%). Like the traditional fecal indicator families Enterobacteriaceae and Enterococcaceae, the relative abundances of alternative feces-associated families were statistically higher in treated sewage from WWTPs with primary treatment only (ANOSIM, R = 0.4685, p-value < 0.001). Sewage-associated genera (i.e., Acinetobacter, Arcobacter, and Trichococcus) were main components of treated sewage as well (i.e., up to 21.7% in Katalinića treated sewage), with higher relative abundances in the absence of a secondary treatment (ANOSIM, R = 0.4839, p-value < 0.001) but significant differences between Francavilla and Zadar treated sewage (ANOSIM, R = 0.2058, p-value < 0.05).

Effect of WWTPs on Marine Community Structure
Seawater samples collected next to WWTP outfalls showed a mild differentiation based on sampling site or type of treatment performed in the corresponding WWTP (Figure S3). The bacterial community structure was mainly dominated by Proteobacteria (45.1-

Effect of WWTPs on Marine Community Structure
Seawater samples collected next to WWTP outfalls showed a mild differentiation based on sampling site or type of treatment performed in the corresponding WWTP ( Figure S3). The bacterial community structure was mainly dominated by Proteobacteria (45.1-50.2%), Bacteroidota (16.6-21.6%), Cyanobacteria (9.2-22.6%), and in a lower contribution by Actinobacteriota (2.0-4.1%) and Planctomycetota (1.3-4.5%). Firmicutes accounted for relative abundances less than 1%, except for samples collected next to the outfall pipe of Katalinića WWTP (i.e., average relative abundance 4.9%) and for two other seawater samples ( Figure S2). None of the phyla that contributed the most to the within-group similarities (i.e., >70% as assessed by SIMPER) varied in a statistically significant fashion.
Abundances of ASVs affiliated with traditional and alternative fecal indicator bacteria in seawater samples collected next to the discharging points of WWTPs are shown in Figure 3. Traditional fecal indicator bacteria accounted for relative abundances up to 0.1% (i.e., the sum of taxa per sample), with differences between the type of sewage treatment performed by the WWTP discharging in the area and with Katalinića seawater samples with the highest relative abundances over time (ANOSIM, R = 0.278, p-value < 0.05). However, we did not find any read affiliated with traditional fecal indicator bacteria in a few seawater samples. Enterococcaceae were significantly lower than Enterobacteriaceae (as observed for WWTP samples) and absent in the majority of the seawater samples. On the contrary, alternative fecal indicator taxa (both feces-and sewage-associated) were found in all seawater samples collected in this study. As observed for the treated sewage, Bacteroidaceae, Lachnospiraceae, and Ruminococcaceae were the main feces-associated families contributing to fecal pollution ( Figure 3B). The distribution of feces-and sewage-associated taxa was not affected by the type of sewage treatment performed by the WWTP in each area, but it was site-specific (ANOSIM, R = 0.265, p-value < 0.05), with the Katalinića seawater samples with the highest contribution of both feces and sewage-associated families (i.e., average overtime was 4.3% and 2.7%, for feces-and sewage-associated taxa, respectively).
We screened the treated sewage samples for a list of 60 genera of PPB (679 ASVs), and then we looked at their presence and distribution in seawater samples. Out of the 42 PPB genera found in the dataset (Figure 4), Alcaligenes, Propionibacterium, Ralstonia, Leptospira, and Serratia were not found in the marine samples. Acinetobacter, Coxiella, Prevotella, Streptococcus, Pseudomonas, Vibrio, Empedobacter, Paracoccus, and Leptotrichia were among the most abundant. Conversely, FIB Escherichia-Shigella and Enterococcus ranged between 0.002-0.03% and 0-0.003%, respectively. The overall distribution of investigated PPB genera differentiated seawater samples based on the type of treatment performed in WWTPs discharging next to seawater sampling points (ANOSIM 0.2216, p-value < 0.05), with Katalinića seawater samples characterized by the highest relative abundances for several genera of PPB. However, despite differences in PPB distribution between samples, one-way ANOVA (on ranks, Kruskal-Wallis test, fdr correction applied) did not highlight any statistical difference associated with a specific site or type of treatment (performed in the WWTP discharging next to seawater sampling points). On the contrary, a differential abundance analysis based on the negative binomial distribution suggested that Acinetobacter was strongly associated with the seawater affected by the discharge from WWTP performing only primary treatment of sewage (Benjamini-Hochberg corrected p-values < 0.001), while Sphingomonas, Legionella, Vibrio, Brevundimonas, and Pseudomonas were associated with the presence of a secondary treatment (Benjamini-Hochberg corrected p-values < 0.001). The large majority of ASVs belonging to the PPB genera found in the seawater samples were also found in treated sewage (including Pseudomonas, which had only two ASVs found in seawater samples only). However, most ASVs belonging to Coxiella were found only in the seawater samples (96 out of 104). Similarly, several ASVs out of the 21 affiliated with Vibrio were exclusive for marine samples.

ARGs in Treated Sewage and Receiving Marine Water
ermB, qnrS, sul2, and tetA were found in the highest concentrations in all WWTPs effluents here investigated, with tetA lower than the other three genes by one order of magnitude or more ( Figure 5). Katalinića and Stobreč WWTPs were characterized by the presence of those four genes in all the surveyed months, unlike Francavilla and Zadar ones. ermB was found in concentrations statistically higher in Katalinića and Stobreč effluents than Francavilla and Zadar ones (H = 32.50548, p-value < 0.001). qnrS did not vary in a statistically significant manner among all tested treated sewage (H = 5.160097, p-value = 0.16). On the contrary, sul2 was significantly higher in Zadar and Francavilla effluents (H = 7.509246, p-value < 0.05), while there were no differences in the content of tetA among treated sewage from Zadar, Francavilla and Katalinića WWTPs (H = 7.509246, p-value = 0.0573), but tetA was significantly lower in Stobreč effluent (H = 6.706772, p-value < 0.01). The large majority of ASVs belonging to the PPB genera found in the seawater samples were also found in treated sewage (including Pseudomonas, which had only two ASVs found in seawater samples only). However, most ASVs belonging to Coxiella were found only in the seawater samples (96 out of 104). Similarly, several ASVs out of the 21 affiliated with Vibrio were exclusive for marine samples.

ARGs in Treated Sewage and Receiving Marine Water
ermB, qnrS, sul2, and tetA were found in the highest concentrations in all WWTPs effluents here investigated, with tetA lower than the other three genes by one order of magnitude or more ( Figure 5). Katalinića and Stobreč WWTPs were characterized by the presence of those four genes in all the surveyed months, unlike Francavilla and Zadar ones. ermB was found in concentrations statistically higher in Katalinića and Stobreč effluents than Francavilla and Zadar ones (H = 32.50548, p-value < 0.001). qnrS did not vary in a statistically significant manner among all tested treated sewage (H = 5.160097, p-value = 0.16). On the contrary, sul2 was significantly higher in Zadar and Francavilla effluents (H = 7.509246, p-value < 0.05), while there were no differences in the content of tetA among treated sewage from Zadar, Francavilla and Katalinića WWTPs (H = 7.509246, p-value = 0.0573), but tetA was significantly lower in Stobreč effluent (H = 6.706772, p-value < 0.01). Genes conferring resistance to beta-lactams, blaCTX-M, blaTEM, and blaOXA, were found in concentrations from 10 to 100 times lower than the other ARGs. blaCTX-M and blaOXA in treated sewage from Francavilla and Zadar WWTPs were negative or <LOQ, with the sole exception of 1 positive sample (each) in the Zadar WWTP. blaOXA was found in quantifiable concentrations only in the Katalinića treated sewage, while in the Stobreč one, it was always positive but <LOQ. blaTEM also was basically absent in Francavilla treated sewage, while it was detected in all the samples from the other three WWTPs, despite being not quantifiable for Katalinića and Stobreč ones due to primer-specificity issues. The colistinresistance gene, mcr-1, was only found in four samples of treated sewage, all <LOQ, but never found in samples collected at Francavilla WWTP. Figure 6 shows ARG presence/absence distribution (% of samples positive for a given gene) in marine sites investigated here due to the low number of observations as quantifiable abundances (i.e., we found several samples to be <LOQ or positive but not quantifiable due to specificity issues, as shown in Figure S4). As observed for WWTP samples, Genes conferring resistance to beta-lactams, bla CTX-M , bla TEM , and bla OXA , were found in concentrations from 10 to 100 times lower than the other ARGs. bla CTX-M and bla OXA in treated sewage from Francavilla and Zadar WWTPs were negative or <LOQ, with the sole exception of 1 positive sample (each) in the Zadar WWTP. bla OXA was found in quantifiable concentrations only in the Katalinića treated sewage, while in the Stobreč one, it was always positive but <LOQ. bla TEM also was basically absent in Francavilla treated sewage, while it was detected in all the samples from the other three WWTPs, despite being not quantifiable for Katalinića and Stobreč ones due to primer-specificity issues. The colistin-resistance gene, mcr-1, was only found in four samples of treated sewage, all <LOQ, but never found in samples collected at Francavilla WWTP. Figure 6 shows ARG presence/absence distribution (% of samples positive for a given gene) in marine sites investigated here due to the low number of observations as quantifiable abundances (i.e., we found several samples to be <LOQ or positive but not quantifiable due to specificity issues, as shown in Figure S4). As observed for WWTP samples, ermB, qnrS, sul2, and tetA were found in higher concentrations and more constantly than the other ARGs. On the contrary, the most abundant gene in marine samples was sul2, with concentrations up to 1.6 × 10 −3 copies/16S rDNA copy (Stobreč, February 2020). Interestingly, mcr-1 was rarely found and in concentration usually <LOQ, although few seawater samples collected next to Francavilla and Zadar WWTP discharging points were higher than LOQ, with average values of 1.9 × 10 −4 and 1.4 × 10 −4 copies/16S rDNA copy, respectively.

respectively.
Seawater samples collected next to Francavilla WWTP discharging point were often negative or <LOQ for most ARGs, except for qnrS and mcr-1 that had an average content of 1.4 × 10 −4 and 1.9 × 10 −4 copies/16S rDNA copy, respectively. Seawater samples collected next to Katalinića and Stobreč WWTP discharging points were characterized by the highest concentrations of ARGs, with sul2 average content of 2.4 × 10 −4 and 4.5 × 10 −4 copies/16S rDNA copy, respectively, qnrS average content of 1.5 × 10 −4 and 1.2 × 10 −4 copies/16S rDNA copy, respectively, and with ermB positive but not quantifiable in Katalinića seawater samples and having an average content of 6.7 × 10 −4 copies/16S rDNA copy in Stobreč ones. tetA was undetectable or <LOQ in Stobreč seawater samples, while it was found in an average content of 7.8 × 10 −5 copies/16S rDNA copy in Katalinića once. While seawater samples collected next to Zadar discharging point were often positive for the majority of the ARGs (Figure 6), concentrations were frequently <LOQ. In particular, tetA and qnrS were always undetectable or <LOQ, while ermB and sul2 were quite variable, with max values of 3.2 × 10 −4 and 2.0 × 10 −5 copies/16S rDNA copy, respectively. Figure 6. Sample prevalence ratio (% of samples positive for a given gene) of ARGs in seawater samples. As the content of the ARGs investigated in this study was often not quantifiable, observations were converted as presence/absence data and then plotted as prevalence ratio. Seawater samples collected next to Francavilla WWTP discharging point were often negative or <LOQ for most ARGs, except for qnrS and mcr-1 that had an average content of 1.4 × 10 −4 and 1.9 × 10 −4 copies/16S rDNA copy, respectively. Seawater samples collected next to Katalinića and Stobreč WWTP discharging points were characterized by the highest concentrations of ARGs, with sul2 average content of 2.4 × 10 −4 and 4.5 × 10 −4 copies/16S rDNA copy, respectively, qnrS average content of 1.5 × 10 −4 and 1.2 × 10 −4 copies/16S rDNA copy, respectively, and with ermB positive but not quantifiable in Katalinića seawater samples and having an average content of 6.7 × 10 −4 copies/16S rDNA copy in Stobreč ones. tetA was undetectable or <LOQ in Stobreč seawater samples, while it was found in an average content of 7.8 × 10 −5 copies/16S rDNA copy in Katalinića once. While seawater samples collected next to Zadar discharging point were often positive for the majority of the ARGs (Figure 6), concentrations were frequently <LOQ. In particular, tetA and qnrS were always undetectable or <LOQ, while ermB and sul2 were quite variable, with max values of 3.2 × 10 −4 and 2.0 × 10 −5 copies/16S rDNA copy, respectively.
bla TEM was relevant only in samples collected next to Katalinića discharging point, with all samples positive despite not being quantifiable due to specificity issues (tested by melting curve) except for one <LOQ. bla OXA and bla TEM were found only in two out of the five seawater samples collected next to Zadar WWTP and with concentrations <LOQ. bla CTX-M gene was negative in all areas investigated.

Discussion
In this study, we sampled treated sewage from four WWTPs (two performing only basic mechanical treatments and the other two performing bio-oxidation in activated sludge tanks) over at least six months, and we investigated their role as a source of microbial pollution in marine coastal areas. Treated sewage that underwent a primary treatment exclusively was characterized by a higher content of culturable FIB (including the opportunistic pathogen Pseudomonas aeruginosa) than treated sewage that underwent an additional secondary treatment, as well as by microbial communities with higher relative abundances of fecal indicator taxa (both feces-and sewage-associated). As expected, in the presence of a final disinfection step (i.e., with peracetic acid, Francavilla WWTP), the number of culturable E. coli, enterococci, and P. aeruginosa decreased further and significantly, often up to 0 CFU/100 mL values. However, a comparison between plate numbers and relative abundances obtained by NGS showed that several reads affiliated with genera Escherichia-Shigella, Enterococcus, and Pseudomonas could still be found in the treated sewage from Francavilla WWTP and may also account for numbers similar to the ones obtained for the other WWTP with secondary treatment (i.e., Zadar WWTP; Figure 7). blaTEM was relevant only in samples collected next to Katalinića discharging point, with all samples positive despite not being quantifiable due to specificity issues (tested by melting curve) except for one <LOQ. blaOXA and blaTEM were found only in two out of the five seawater samples collected next to Zadar WWTP and with concentrations <LOQ. blaCTX-M gene was negative in all areas investigated.

Discussion
In this study, we sampled treated sewage from four WWTPs (two performing only basic mechanical treatments and the other two performing bio-oxidation in activated sludge tanks) over at least six months, and we investigated their role as a source of microbial pollution in marine coastal areas. Treated sewage that underwent a primary treatment exclusively was characterized by a higher content of culturable FIB (including the opportunistic pathogen Pseudomonas aeruginosa) than treated sewage that underwent an additional secondary treatment, as well as by microbial communities with higher relative abundances of fecal indicator taxa (both feces-and sewage-associated). As expected, in the presence of a final disinfection step (i.e., with peracetic acid, Francavilla WWTP), the number of culturable E. coli, enterococci, and P. aeruginosa decreased further and significantly, often up to 0 CFU/100 mL values. However, a comparison between plate numbers and relative abundances obtained by NGS showed that several reads affiliated with genera Escherichia-Shigella, Enterococcus, and Pseudomonas could still be found in the treated sewage from Francavilla WWTP and may also account for numbers similar to the ones obtained for the other WWTP with secondary treatment (i.e., Zadar WWTP; Figure 7).  On the contrary, the 16S rRNA relative abundances and the plate count data of Katalinića, Stobreč, and Zadar treated sewages followed comparable patterns. Although metagenomics investigations do not provide evidence of activity or for bacteria to be alive, regrowth experiments demonstrated that disinfection processes could inactivate or damage bacteria while not necessarily killing them [63,71,72]. Moreover, cultivation methods can overestimate the efficiency of disinfection treatment as damaged cells can lose culturability for short or long periods [7,73]. The discrepancy between NGS data and plated count in Francavilla treated sewage was particularly evident for the Pseudomonas genus, an emerging pathogen of environmental relevance that can survive and adapt to a diversity of environments [34,74]. Alternative fecal indicator taxa (i.e., five bacterial families constituting up to 85% of total sequences in human and animal feces, and three genera very abundant in urban WWTPs while not prevalent in human feces) covered high proportions of the bacterial diversity in all treated sewages here investigated, despite differences between types of treatments.
Seawater samples collected next to WWTP discharging points revealed a marked fecal signature, despite evident mitigation due to the dilution of the treated sewage plumes. Alternative fecal indicator taxa provided a better understanding of the actual fecal bacterial load in marine samples compared to traditional bacterial indicators that could not even be detected into the microbiome. In addition, NGS showed the presence of potentially pathogenic bacteria, especially Acinetobacter, Coxiella, Prevotella, Streptococcus, Pseudomonas, Vibrio, Empedobacter, Paracoccus, and Leptotrichia. The large majority of PPB phylotypes were shared between treated sewage and seawater samples, except for Vibrio and Coxiella, suggesting that WWTPs may represent the primary source for PPB in marine areas investigated in this study. Indeed, the presence of Vibrio, Pseudomonas, and Coxiella phylotypes in marine waters is not uncommon due to the ubiquitous nature of representatives of these genera [75][76][77]. However, deeper investigations and the use of long-read sequencing approaches may reveal differences that could not be appreciated with the methodologies chosen for the purposes of this study. Seawater samples collected next to Katalinića WWTP discharging points showed the highest content in traditional and alternative fecal indicator bacteria, while statistically significant differences with all the other sites, suggesting that the local hydrographic characteristics could contribute to a scarce dilution of the sewage plume.
WWTPs in this study represented continuous sources of ermB, qnrS, sul2, and tetA genes for coastal marine water. The three Croatian WWTPs (both with and without secondary biological treatment) were also important sources of the bla TEM gene, which codes for one of the most common beta-lactamase enzymes in Gram-negative bacteria, thus conferring resistance to ampicillin [78]. Contrary to what was expected, the bla CTX-M gene was found in low concentrations in the Katalinića WWTP effluent and in <LOQ or null concentrations in the other WWTPs. It has been hypothesized that urban WWTPs are the principal source of bla CTX-M for aquatic environments [23,51] because bla CTX-M harboring strains isolated in rivers and wastewater were found to be related to nosocomial strains of the same geographical area [79,80].
Although the two WWTPs with only primary treatments (i.e., Katalinića and Stobreč) were characterized by a more constant presence of ARGs over time, differences between genes were often not associated with the type of treatment (i.e., primary vs. secondary), except for sul2 (significantly higher in the presence of secondary treatment) and ermB (significantly lower in the presence of secondary treatment). Differences observed in this study were gene-specific and were likely related to the geographical origin of the corresponding wastewater streams. Francavilla WWTP was characterized by the absence or by low concentrations of ermB, bla TEM , bla OXA , and bla CTX-M . Although sewage in Francavilla WWTP undergoes final disinfection with peracetic acid, the disinfection by chemicals (e.g., chlorination or with peracetic acid) is known to be effective in minimizing the release of pathogens, but they would not affect the concentration of ARGs in treated sewage (when normalized to the 16S copy number) [7]. On the contrary, disinfection may favor the co-selection of antibiotic-resistant bacteria, although results of previous investigations were sometimes contrasting [7,63,71,81,82]. In particular, it was observed that chlorination could promote HGT events [83,84], which may favor AMR spreading and emerging of new AMR mechanisms. Due to its small size (i.e., 38,000 PE) Francavilla WWTP receives a lower inflow volume than the other WWTPs investigated in this study, which might relate to the observed lower concentrations for some ARGs.
One of the aims of this study was the assessment of a possible ARG chronic pollution into marine coastal waters affected by WWTP effluent discharge. To our knowledge, this is the first study investigating the presence over time of ARGs in the Adriatic Sea and one of the few ones in the Mediterranean Sea. Di Cesare and co-authors [85] investigated the presence of ermB, qnrS, tetA, and bla TEM in marine sediments collected from the Pula Bay (Croatia, Adriatic Sea), one of the most significantly polluted sites along the Croatian coast, and found that only ermB and bla TEM were detectable in some of the sampling points and that they generally remained below the quantification limits. Šamanić and coauthors [86] cultured chlorine-resistant bacteria from surface water collected from Croatian public beaches and found isolates resistant to clinically significant β-lactamases, such as bla CTX-M-15 , bla CTX-M-3 , bla SHV-12 , and bla TEM .
As observed for the treated sewage, ermB, qnrS, sul2, and tetA were found at sea in higher concentrations and more constantly than the other ARGs, despite the evident mitigation due to dilution. Only Katalinića seawater samples were always positive to ermB, qnrS, sul2, tetA, and bla TEM , probably because of the hydrographic characteristics of the area. We have not found bla CTX-M in seawater samples, likely due to its low levels in the treated sewage. mcr-1 was positive in few samples collected next to the Zadar, Stobreč, and Francavilla WWTP discharging points. Interestingly, the few positive Zadar and Francavilla seawater samples had higher concentrations than their corresponding effluents. That could suggest that other environmental sources are relevant for mcr-1 (and, thus, other ARGs) in marine environments. The current mcr-1 distribution has been achieved through multiple translocations. A likely driver for the global spread is the trade (i.e., livestock animals and meat), but poorly treated sewage and aquaculture (including integrated aquaculture systems) have been suggested as probable responsible for the worldwide dissemination of the mcr gene variants in aquatic environments and as major hotspots for genetic recombination and horizontal gene transfer of mcr and mcr-like genes [9,87,88]. To our knowledge, this is the second finding of the plasmid-mediated colistin resistance gene mcr-1 in the Adriatic Sea and one of the first for the Mediterranean Sea. A very recent study [86] has detected the mcr-1 gene in one bacterial strain isolated from a surface water sample collected in a Croatian public beach.

Conclusions
In this study, we investigated WWTPs over time as continuous point sources of microbial pollution for coastal water of the central Adriatic Sea. In the absence of secondary treatments, marine coastal waters received high inputs of fecal indicator bacteria and potentially pathogenic bacteria (including the emerging opportunistic pathogen Pseudomonas aeruginosa). Moreover, WWTPs under study represented continuous sources of genetic elements conferring resistance to major classes of antibiotics, with potential consequences on the emerging and the establishment of AMR pathways in the environmental microbial communities. NGS-based approaches provided a broader detection of specific taxa than classic microbiological techniques, especially in the presence of a final disinfection step, which could lead to an underestimation of fecal pollution. In particular, the exploration in the microbiomes of alternative fecal indicator taxa and potentially pathogenic bacteria allowed a more thorough understanding of emerging and potential pathogens spread into the environment from either low efficiency or conventional WWTPs. Overall, our results highlight the need for an improvement of wastewater treatment technologies, especially for those plants which perform primary treatment only, and stress the importance of dedicated monitoring programs at the WWTP discharge points at sea.