Microbial Biofilms Dynamics and Functionality in an Urban Mycobacterium-Dominated Drinking Water Distribution System

TL;DR

This study explores how biofilms form in urban drinking water systems, focusing on Mycobacterium and their impact on water quality and safety.

Contribution

The study identifies key microbial taxa and novel viral elements involved in biofilm dynamics and antibiotic resistance in a Mycobacterium-dominated system.

Findings

Biofilm communities were more functionally diverse than water communities and correlated with Mycobacterium pathogens and antibiotic resistance genes.

22% of identified microbial taxa were strongly responsible for biofilm formation and matrix remodeling.

New putative Mycobacterium phages were detected that may destabilize biofilms by targeting mycolic acids.

Abstract

Microbial communities in drinking water distribution
systems (DWDS)
develop primarily as biofilms on pipe surfaces. Despite their impact
on water quality, infrastructure maintenance, and biosafety, biofilms
are not routinely controlled. In this study, we investigated the bacterial
community dynamics and functionality in an urban chlorinated DWDS,
dominated by Mycobacterium, through
a multiphasic approach which included 16S rRNA gene metabarcoding,
metagenomics and microscopy. Our results showed that biofilm communities
were more functionally diverse compared to those from water and that
the biofilm maturity was positively correlated with the prevalence
of potential Mycobacterium emerging
pathogens and a broader distribution of antibiotic resistance genes
(ARGs) within the microbial community. The reconstruction of metagenome-assembled
genomes (MAGs) and the corresponding genomospecies allowed the identification
of key microbial taxa involved in the biofilm matrix remodeling, with
22% of them strongly responsible for biofilm formation. A diverse
and novel viral community was detected across the system, including
new putative Mycobacterium phages that
might act against mycolic acids and thus contribute to biofilm destabilization.
Our findings enhance our understanding of DWDS microbial composition
and biofilm formation dynamics, focusing on “who does what”
and then providing a foundation for developing effective biofilm control
strategies in water distribution systems.

Figures

(A) Map of the buildings sampled in the DWDS, for which only one water entry exists. EPDM pipes B and G were located in the same washbasin, with segment G sampled one year after its installation (device G replaced device B). (B) Sampling strategy: after removing the EPDM devices, 50 L of drinking water was sampled and a new EPDM hose was installed. (C–F) Scanning electron microscopy images of biofilm F. Rod-shaped bacteria embedded in the biofilm matrix can be visualized in (C) and (D). Most parts of the EPS matrix observed in (C) and (D) presented a filamentous structure, while the matrix in (E) and (F) was more granular. In each microscopy image, the porosity of the matrix, which allows water to carry nutrients and other compounds to lower the cells through the formed channels, is observed.

(A,B) Alpha diversity metrics (y-axis) derived from the 16S rRNA gene metabarcoding analysis from biofilm and water samples (x-axis). (C) Nonmetric multidimensional scaling (NMDS) analysis based on the relative abundances of the OTUs from biofilm and water samples. Before each sample name, the preffix “BA_” (biofilm amplicon) or “WA_” (water amplicon) was added.

Heatmap of the bacterial phyla (A) and genera (B) that exceeded 1% relative abundance in at least one biofilm or water sample. Those genera present in ≥50% of both biofilms and water samples, with no significant differences in their relative abundances between both types of samples, were considered as “generalist”, while those genera present in biofilms and water samples, with a relative abundance statistically higher in waters were considered as “prevalent in water”. Genera exclusively present in ≥50% water samples were considered as “specific in water”. Genera “prevalent in biofilms” or “specific in biofilms” are not shown as none of them exceeded 1% of relative abundance in any sample. Before each sample name, the prefix “BA_” (biofilm amplicon) or “WA_” (water amplicon) was added.

CARD-FISH images of those phyla with the highest relative abundances in biofilms. Each image is the result of the composite of those cells hybridized with a given probe (in blue, white arrows) and the rest of cells stained with SYBRTM Gold (in gray).

(A) Cumulative relative abundance (TAD80/GE; y-axis) of bacterial phyla in the DWDS metagenomes (x-axis) based on the genomospecies composition. (B) Phylogenomic tree of the obtained 104 genomospecies. In (B), the layers over the genomospecies names are a heatmap of the relative abundance of each genomospecies in the different biofilm stages (incipient, intermediate, and old) and the water samples. For the relative abundances in intermediate and old biofilms, average value were calculated with the data of the 3 metagenomes of each stage. Red dots indicate the presence of, at least, one ARGs in each genomospecies. Branch colors are associated with genomospecies considered as “biofilm producers” (in brown) or the main “biofilm degraders” (in green; those with ≥1% of their genes associated with EPS degradation).

Tables

1. CAZy Families Predicted among the vOTUs

Predicted CAZy families	Predicted CAZy activities	vOTU id	Family	Putative host
CE5	Triacylglycerol lipase (EC 3.1.1.3); Acetylxylan esterase (EC 3.1.1.72); Cutinase (EC 3.1.1.74)	vOTU67*	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU84*	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
GH104	Lysozyme (EC 4.2.2.29)	vOTU41	Caudoviricetes	Genus Bosea (Pseudomonadota)
GH108	Lysozyme (EC 3.2.1.17)	vOTU523	Caudoviricetes	Genus Bosea (Pseudomonadota)
		vOTU112	Caudoviricetes - Crassvirales	Genus ELB16-189 (Bacteroidota)
		vOTU729	Caudoviricetes	Genus Halomonas (Pseudomonadota)
		vOTU503	Caudoviricetes	Genus Hyphomicrobium (Pseudomonadota)
		vOTU391	Caudoviricetes	Genus Sphingobium (Pseudomonadota)
		vOTU529	Caudoviricetes	Genus Sphingobium (Pseudomonadota)
		vOTU172	Caudoviricetes
		vOTU235	Caudoviricetes
		vOTU329	Caudoviricetes
		vOTU390	Caudoviricetes
		vOTU467	Caudoviricetes
GH113	Mannosidase (EC 3.2.1.-; 3.2.1.25; 3.2.1.78)	vOTU30	Caudoviricetes
		vOTU641	Caudoviricetes
GH16	Glucosidase (EC 3.2.1.6,21,39,58,73,151,...); Galactosidase (EC 3.2.1.23,81,83,103,181,...·)	vOTU96	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
GH18	Chitinases (EC 3.2.1.-,14,200,201); Mannosyl-oligosaccharide endo-b-1,4-N-acetylglucosaminidase (EC 3.2.1.96)	vOTU559	Caudoviricetes	-
GH19	Lysozyme (EC 3.2.1.17); Chitinase (EC 3.2.1.14)	vOTU733	Caudoviricetes	Genus 2-12-FULL-53-21 (Desulfobacterota)
		vOTU531	Caudoviricetes	Genus Bosea (Pseudomonadota)
		vOTU417	Caudoviricetes	Genus Bradyrhizobium (Pseudomonadota)
		vOTU646	Caudoviricetes	Genus Reyranella (Pseudomonadota)
		vOTU680	Caudoviricetes	Genus Reyranella (Pseudomonadota)
		vOTU20	Caudoviricetes
		vOTU405	Caudoviricetes
		vOTU477	Caudoviricetes
		vOTU657	Caudoviricetes
		vOTU674	Caudoviricetes
GH23	Lysozyme (EC 3.2.1.17, 4.2.2.29); Chitinase (EC 3.2.1.14)	vOTU153	Viruses Uncl.	Genus Mycobacterium (Actinomycetota)
		vOTU307	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU338	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU72	Viruses Uncl.	Genus Mycobacterium (Actinomycetota)
		vOTU81	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU33	Caudoviricetes	Genus Sphingopyxis (Pseudomonadota)
		vOTU68	Caudoviricetes
		vOTU701	Caudoviricetes
		vOTU99	Caudoviricetes
GH24	Lysozyme (EC 3.2.1.17)	vOTU670	Caudoviricetes	Genera Altererythrobacter and Erythrobacter (Pseudomonadota)
		vOTU310	Caudoviricetes	Genus Cereibacter (Pseudomonadota)
		vOTU466	Caudoviricetes	Genus Cereibacter (Pseudomonadota)
		vOTU42	Caudoviricetes	Genus Erythrobacter (Pseudomonadota)
		vOTU737	Caudoviricetes	Genus Erythrobacter (Pseudomonadota)
		vOTU716	Corticoviridae (Tectiliviricetes)	Genus Ga0077527 (Pseudomonadota)
		vOTU734	Caudoviricetes	Genus Reyranella (Pseudomonadota)
		vOTU683	Caudoviricetes	Genus Roseococcus (Pseudomonadota)
		vOTU32	Caudoviricetes	Genus Sphingopyxis (Pseudomonadota)
		vOTU414	Caudoviricetes	Genus Tabrizicola (Pseudomonadota)
		vOTU78	Caudoviricetes	Genus Tropicimonas (Pseudomonadota)
		vOTU28	Caudoviricetes
		vOTU29	Caudoviricetes
		vOTU341	Caudoviricetes
		vOTU394	Caudoviricetes
		vOTU437	Caudoviricetes
		vOTU446	Caudoviricetes
		vOTU491	Caudoviricetes
		vOTU512	Caudoviricetes
		vOTU544	Caudoviricetes
		vOTU630	Caudoviricetes
		vOTU649	Caudoviricetes
		vOTU725	Caudoviricetes
		vOTU361	Caudoviricetes
		vOTU379	Caudoviricetes
GH39	Cellulose glucosidase (EC 3.2.1.74); Arabinofuranosidase (EC 3.2.1.215); Rhamnosidase (EC 3.2.1.-); Xylosidase (EC 3.2.1.37)	vOTU115	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU303	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
		vOTU358	Caudoviricetes	Genus Mycobacterium (Actinomycetota)
GH51	Cellulose glucosidase (EC 3.2.1.4,91); Arabinofuranosidase (EC 3.2.1.55,215); Xylosidase (EC 3.2.1.8,37)	vOTU381	Caudoviricetes
GH73	Lysozyme (EC 3.2.1.-); Mannosyl-oligosaccharide endo-b-1,4-N-acetylglucosaminidase (EC 3.2.1.96)	vOTU661	Caudoviricetes
GT101	Hexosyltransferases (EC 2.4.1.-)	vOTU352	Caudoviricetes
GT102	N-acetylglucosaminyl-PP-decaprenol α-1,3-L-rhamnosyltransferase (EC 2.4.1.289)	vOTU554	Viruses Uncl.
GT2	Glycosyltransferase (2.4.1.-,12,16,117,157,212,···)	vOTU114	Caudoviricetes	Genus Bradyrhizobium (Pseudomonadota)
		vOTU546	Viruses Uncl.
		vOTU619	Caudoviricetes
GT25	Galactosyltransferase (EC 2.4.1.-); Xylosyltransferase (EC 2.4.2.-); Acetylgalactosaminyltransferase (EC 2.4.1.-)	vOTU40	Caudoviricetes
GT32	α-glucosyltransferase (EC 2.4.1.-)	vOTU1	Caudoviricetes
		vOTU550	Caudoviricetes
		vOTU557	Caudoviricetes
		vOTU658	Caudoviricetes
		vOTU675	Caudoviricetes
		vOTU703	Caudoviricetes
		vOTU735	Caudoviricetes
GT51	Peptidoglycan glycosyltransferase (EC 2.4.1.129)	vOTU731	Caudoviricetes	Genus SYSU-D60009 (Pseudomonadota)
PL1	Pectate lyase (EC 4.2.2.2,9); Pectin lyase (4.2.2.10)	vOTU361	Caudoviricetes
		vOTU379	Caudoviricetes
		vOTU409	Caudoviricetes

Full text

Introduction

Microbes are present in drinking water distribution systems (DWDS) as part of the autochthonous water microbiota or because they infiltrate as contaminants due to infrastructural deficiencies. While most microbes in DWDS are harmless, those considered pathogenic might pose significant health risks if they are not adequately controlled. Therefore, maintenance of the water quality, infrastructure integrity, and biosafety are priorities in water control programs. Advances in microbiology research have facilitated the monitoring of microbial communities in DWDS, providing a deeper understanding of how the composition of the drinking water microbiome is influenced by the source of the water, the disinfection treatment employed and the characteristics of the pipe material, which seems to be the principal factor determining the community composition. ?−? ?

Although drinking water typically contains planktonic bacterial numbers between 10^3^ and 10^6^ cells/mL,? up to 90% of the total bacterial biomass in DWDS is found in the biofilms that develop in the interfaces between the water and pipe walls, even in the presence of disinfectants or reduced amounts of organic carbon.? The presence of biofilms in DWDS is a well-known primary cause of water quality problems and the biofilm detachment of cells into the bulk water is a major biosafety concern if released bacteria are pathogenic and/or antibiotic resistant. ?,? Biofilms also participate in pipe corrosion, weakening metal structures, producing undesirable tastes and smells to drinking water, and facilitating the entry of external contaminants into the system. Finally, biofilm formation and growth can lead to pipe obstruction and reduced efficiency of the water treatment processes. This biofouling hinders water flow and requires frequent cleaning or even replacement of system components, resulting in increased costs for water utilities.?

However, in spite of the biofilm impact, regulations of water quality and safety of urban DWDS do not consider the control of biofilms and focus solely on determinations of fecal microbial markers, pathogenic bacteria such as Legionella or Clostridium, and total heterotrophic bacteria plate counts as the regular parameters.? Additionally, scientific studies analyzing biofilms from real urban networks are scarce? and, to our knowledge, none of them have employed state-of-the-art metagenomic analyses with the aim of identifying which members of the DWDS microbial communities can be associated to the formation, evolution and persistence of the biofilm. Such information is essential to develop safe and effective treatments to prevent or disrupt biofilms. Indeed, European Regulation EU 2020/2184 requires the implementation of measures that directly impact the formation and control of biofilms. This issue should be addressed within the Water Safety Plans of drinking water distribution systems.

In this work, an in-depth multiphasic approach, including 16S rRNA gene amplicon sequencing (metabarcoding), metagenomics, CARD-fluorescence in situ hybridization and analytical chemistry, has been applied to study the dynamics and functional aspects of the microbial communities in an urban Mycobacterium-dominated DWDS, including virus-host interactions. Our results show how unconnected areas of the same DWDS harbor similar microbial communities and to what extent the degree of biofilm maturity is intimately tied to an increasing abundance of potential pathogens and ARGs, which might constitute targets for the development of molecular diagnostic tools for rapid decision making on water sustainability plans. Statistical analyses allowed distinguishing which bacteria were “specific” or “prevalent” in biofilms or in the planktonic fraction. In addition, we provide insights into relevant functions linked to different biofilm stages based on the reconstruction of metagenome-assembled genomes (MAGs). Thus, this study advances our understanding of the microbial dynamics and interactions in DWDS biofilms, which has implications on water quality and network management by establishing modern water safety plans.

Materials and Methods

Sampling and Physicochemical Determinations of Water Samples

Twelve biofilms and water samples from six buildings, A–F (FigureA), included in an urban DWDS in the Southeast of Spain and supplied with the same chlorinated mixed water (freshwater from aquifers mixed with desalinated seawater) were obtained in May 2021. Biofilm biomass, developed on ethylene propylene diene monomer (EPDM) hoses that connected the DWDS to individual washbasins (FigureB), was collected with a sterile spatula and kept at −80 °C (Table S1). While the time in place for devices A-F was unknown, an additional hose G was sampled one year after its installation. In addition, 50 L of water from each sampling point were filtered through 0.22-μm pore-size polyvinylidene fluoride membranes (Millipore) and filters kept at −80 °C. A physicochemical analysis was performed with the water samples by the company LABAQUA, S.A.U. accredited under ISO 17025. Data were obtained for 20 variables: temperature (°C), conductivity, pH, oxidability, total organic carbon, chloride, sulfur, sulfate, nitrogen, nitrate, ammonia, aluminum, barium, cobalt, copper, manganese, nickel, lead, selenium, and zinc.

(A) Map of the buildings sampled in the DWDS, for which only one water entry exists. EPDM pipes B and G were located in the same washbasin, with segment G sampled one year after its installation (device G replaced device B). (B) Sampling strategy: after removing the EPDM devices, 50 L of drinking water was sampled and a new EPDM hose was installed. (C–F) Scanning electron microscopy images of biofilm F. Rod-shaped bacteria embedded in the biofilm matrix can be visualized in (C) and (D). Most parts of the EPS matrix observed in (C) and (D) presented a filamentous structure, while the matrix in (E) and (F) was more granular. In each microscopy image, the porosity of the matrix, which allows water to carry nutrients and other compounds to lower the cells through the formed channels, is observed.

Cell and Virus-like Particle (VLPs) Counts in Water Samples

Ten mL water samples B, D and G were fixed with 4% formaldehyde (final concentration) and the fixation stopped by adding sterile phosphate-buffered saline (PBS) 1X. Four mL fixed samples were filtered through 0.22 μm pore-size isopore GTTP filters (Millipore), for cell counts, and 0.02 μm pore-size Anodisc (Whatman) filters, for VLP determinations. Subsequently, isopore GTTP filters were stained with 4′,6-diamidino-2-phenylindole (DAPI) and Anodisc filters were stained with SYBR^TM^ Gold (25X). Cells and VLP were counted on a Leica DM4000B epifluorescence microscope. In parallel, cell and VLP concentrations were also determined by flow cytometry. We followed the protocol described by Brussaard.? Two mL of each selected water sample was fixed with glutaraldehyde (0.5% final concentration), frozen in liquid nitrogen, diluted in Tris-EDTA buffer (pH 8), stained with SYBR^TM^ Gold (0.5X final concentration), incubated at −80 °C, and maintained for 5 min at room temperature prior to analysis. The cytometer settings were as follows: the threshold was set in blue fluorescence (300 units), FITC voltage = 500, SSC voltage = 300, and forward scatter voltage = 500, and the flow rate was established as low. Background noise was checked on blanks, composed by TE buffer stained with SYBR^TM^ Gold. Samples were recorded at an event rate of 100–1000 events per second. Cells and VLPs counts were obtained by correcting the noise measured in blanks and expressed as cells/ml or VLPs/ml.

Catalyzed Reported Deposition-Fluorescence in Situ Hybridization

(CARD-FISH) and Scanning Electron Microscopy (SEM)

We followed the protocol described in Pernthaler et al.? and adapted by Escudero et al.? EPDM devices were cut in small pieces, without damaging the structure of the biofilm, and the samples were washed in ethanol 90% and air-dried. Cell permeabilization was determined with lysozyme and different hybridization and washing buffers were used depending on the probe (Table S2). Biofilms were counterstained with SYBR^TM^ Gold (25x), and preparations were mounted onto glass slides and covered with Citifluor (Citifluor, United Kingdom). The hybridized biofilms were observed with a ZEISS confocal laser scanning microscope model LSM 800, coupled with an Axio Imager Z.2 straight microscope and equipped with three laser diode excitation sources: 405 nm (UV), 488 nm (blue), 561 nm (green) and 640 nm (red). Images were collected with 63X oil immersion lens and the software ZEISS ZEN 3.9 was used to project the stacks to 2D images. As for SEM microscopy, first, the sample F was dehydrated with ethanol at different concentrations: 25% (5 min), 50% (5 min), 75% (5 min), 95% (5 min) and 100% (10 min). Then, a thin layer of platinum was sputtered onto the samples to enhance conductivity and prevent charging under the electron beam. The equipment used was a FESEM (Field Emission Scanning Electron Microscope) ZEISS model Merlin VP Compact, and the resolution was achieved at 5 kV.

Characterization of the Biofilms’ Extracellular Polymeric

Substances (EPS)

Biofilms with the highest biomass weight were used (D and F), and two strategies were followed for the detection and identification of sugars and other compounds: (i) total acid hydrolysis with trifluoroacetic acid (TFA) followed by High-Performance Anion-Exchange Chromatography with Pulsed Amperometric Detection (HPAEC-PAD); (ii) total acid hydrolysis with TFA followed by Gas Chromatography Mass Spectrometry (GC-MS).

Biofilms with the highest biomass were used (D and F). First, 300 mg of biofilm biomass was mixed with 1 mL of ultrapure water and sonicated for 15 min. After centrifuging at 10,000g for 10 min, the supernatant was freeze-dried. Part of the dry biomass (100 mg) was used for Nuclear Magnetic Resonance (NMR) analysis, and the rest of the sample for total acid hydrolysis with trifluoroacetic acid (TFA), following the indications described in Pettolino et al.? For TFA hydrolysis, 100 μL of 2.5 M TFA was added to the dry samples and the mixture was incubated at 100 °C in a thermomixer, at 300 rpm for 4 h. Then, the preparations were cooled and placed at 30 °C, and the TFA was evaporated with a stream of nitrogen. Half of the hydrolyzed samples was used for High-Performance Anion-Exchange Chromatography with Pulsed Amperometric Detection (HPAEC-PAD) analysis, using different standards at 0.05 g/L (such as rhamnose, glucuronic acid, fucose, glucose, galactose, mannose, glycerol, maltose, glucosamine, fructose, sucrose, and N-acetyl-d-glucosamine). The compound separation was achieved at 30 °C for 25 min on a Dionex CarboPac PA-1 column (4 × 200 mm), equipped with its corresponding guard column, equilibrated with 4.5 mM NaOH, and running at a flow rate of 1 mL/min. The other half of the hydrolyzed samples was reduced and acetylated as explained in Pettolino et al.? The final products were dissolved in 500 μL of dichloromethane (DCM) and transferred to a glass vial for further analysis by gas chromatography–mass spectrometry (GC-MS) in a GC-MS-Orbitrap Q Exactive system (Thermo Scientific).

DNA Extractions, 16S rRNA Gene Amplification and Sequencing

DNA was extracted from biofilms and the 0.22-μm pore-size membranes using the FastDNA Spin Kit for soil (MP Biomedicals, Santa Ana, CA, USA), following the manufacturer’s recommendations, and fluorometrically quantified using Qubit 2.0 (Thermo Fisher Scientific) and High Sensitivity Assay (HS) kits. For 16S rRNA gene metabarcoding, the V3-V4 region of the 16S rRNA genes was PCR amplified with primers Pro341F/Pro805R.? The PCR program employed was as follows: denaturation at 94 °C for 3 min, 30 cycles of denaturation (94 °C, 45 s), hybridization (51 °C, 1 min) and elongation (72 °C, 2 min), and a final extension at 72 °C for 10 min. PCR products were electrophoresed in a 1% agarose gel in TAE 1X buffer (40 mM Tris at pH 8.0, 20 mM acetic acid, 1 mM EDTA). PCR products were sequenced at FISABIO (Valencia, Spain) using Illumina Mi-Seq (2 × 300 bp). For metagenomes, extracted DNAs were sequenced at Macrogen (Seoul, South Korea) using a NextSeq500 (2 × 150 bp, 4 Gb/sample).

Processing of 16S rRNA Gene Sequences

The 16S rRNA gene sequences from biofilm and water samples were analyzed with QIIME 2 2021.4.? First, raw reads were demultiplexed and filtered by quality using the q2-cutadapt? and q2-quality-filter? plugins. Then, q2-vsearch,? with the de novo clustering option, was used to generate Operational Taxonomic Units (OTUs) at 99% identity. Chimeras were removed with the option uchime-denovo, and final relative abundances were obtained with the plugin q2-feature-table; the option rarefy in this plugin was used for normalization according to the number of reads in the smallest library (sample WA_C, with 38,650 reads). Rarefaction curves were calculated and plotted with the Vegan Package? in R v4.0.3 (www.r-project.org). Then, alpha diversity parameters (richness and Shannon’s indexes) were calculated using q2-diversity. Finally, the OTUs taxonomic affiliation was obtained by comparing the reference sequence of each OTU against the SILVA v138 SSU 16S rRNA database with the q2-feature-classifier? and classify-sklearn? with a confidence threshold of 70%.

Analysis of Metagenomes

Metagenomic raw reads from water and biofilm samples were first processed with Trimmomatic v0.36 to remove the Illumina adapters and identify pair-end reads (LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36).? Then, the FastQC program was used to analyze the quality of the trimmed reads, the amount of Ns, the frequency of duplicated sequences and the presence of adapters.? Second, the sequencing effort and the diversity of the metagenomes were determined with the Enveomics Nonpareil program.? In parallel, reads belonging to 16S rRNA genes were extracted using SortMeRNA v4.3.6? and taxonomically sorted by BLASTn against the SILVA v138.1 database,? selecting only the best hit with the BlastTab.best_hit_sorted.pl option from Enveomics.? Only those hits with coverages and identities higher than 70% and 83%, respectively, were considered.?

To assemble the metagenomic trimmed reads, several softwares were tested: IDBA-UD, metaSPAdes and MEGAHIT. ?−? ? Contigs obtained with metaSPAdes v3.13.0 presented the best quality; therefore, they were used for subsequent analyses. The prediction of Open Reading Frames (ORFs) and the extraction of the amino acid sequences were performed with Prodigal and the -p meta option.?

For the analysis of metabolic pathways, contigs were annotated with Anvi’o? and the KEGG database. ?,? The completeness of each KEGG module and the number of KO-annotated proteins was obtained with “anvi-estimate-metabolism” and the options “--kegg-output-modes modules” and “--kegg-output-modes kofam_hits”. Only those KEGG modules with a completeness of ≥75% were considered. To verify that the observed differences between biofilms and water metagenomes were not attributed to the higher sequencing depth of the biofilm samples, the same analyses were performed after normalizing the metagenomes to the same number of sequences. In parallel, the standalone tool run_dbcan of dbCAN3? was used to predict Carbohydrate-Active enZymes encoding genes (CAZymes).? Run_dbcan is based on HMMER (e-value < 1e-15, coverage > 0.35) and DIAMOND BLASTp (e-value < 1e-102) analyses. Only those proteins correctly identified by both strategies were considered as CAZymes, as previously recommended.?

In addition, we specifically looked for those ORFs directly related to EPS synthesis and degradation. On the one hand, the detection of metabolic pathways and KOs involved in the production of exopolysaccharides (KEGG pathway map00543) and mycolic acids? was performed. On the other hand, specific enzymes related to the synthesis or degradation of the most common exopolysaccharides detected in biofilms (such as alginate, Pel, Psl, cellulose, Poly-N-acetyl-glucosamine and xanthan, among others?), were searched by BLASTp (BlastTab.best_hit_sorted.pl, coverage > 50%, identity > 40%) against a homemade database containing 3,659 entries of glycosidases (EC 3.2.1.-: glycosidases, i.e. enzymes hydrolyzing O- and S-glycosyl compounds) and 535 entries of enzymes acting on polysaccharides (EC 4.2.2.-: lyases acting on polysaccharides), retrieved from the ENZYME repository? and UniProtKB? (Text S7).

Finally, proteins related to antibiotic resistance genes (ARGs) were annotated with the Resistance Gene Identifier (RGI) tool, based on the reference data from the Comprehensive Antibiotic Resistance Database (CARD), ?,? while Queuosine (Q) precursors genes were evaluated following the indications of Díaz-Rullo and González-Pastor,? using the COG database, where these genes are referred as queD (COG0720), queE (COG0602), queC (COG0603), queF (COG0780 and COG2904), tgt (COG0343), queA (COG0809) and queG/queH (COG1600 and COG1636).

Reconstruction and Analysis of Metagenome-Assembled Genomes

(MAGs)

To generate metagenome-assembled genomes (MAGs), optimization of the binning process was first carried out. For the optimization, two different programs (MaxBin2 v2.2.7,? and MetaBAT2 v2.15?) and several minimum contig sizes (1 kb, for MetaBAT2 only; 1.5 kb, for MaxBin2 only; 2 kb, 2.5 kb, 3 kb, 5 kb and 10 kb) were used. Completeness and contamination were determined for each MAG by CheckM2? and for the subsequent analyses, the optimal contig size for each program and metagenome was selected based on the highest number of MAGs with completeness above 80% and contamination below 10%, considering that if the same number of MAGs with two different contig sizes were obtained, the largest contigs were selected. The selected MAGs were processed with the DAS Tool,? that refines and selects nonredundant MAGs of good quality from a set of MAGs previously reconstructed with other binning algorithms. Then, in order to eliminate low quality MAGs, only those whose Q50 value, defined by [%completeness] – (5 × [%contamination]), was >50 were selected for further analysis.? In addition, a manual refinement of the best quality MAGs (completeness

80%) based on protein taxonomy of each contig (BLASTp vs NR database?) and their sequencing depth (coverage > 70%, identity > 95%) was performed, as previously described.? Finally, MAGs were clustered in genomospecies at 95% ANI (Average Nucleotide Identity)? with minimum aligned fraction of 30% (default value) using the dRep program,? and functional annotations were performed as described above for metagenomic contigs, with the representative MAG.

The relative abundance of each genomospecies was estimated using the highest quality MAG within each genomospecies as the representative. A MAG was considered as detected in a given sample based on the value of breadth coverage calculated using samtools,? assuming a positive detection if the value was greater than 10%. ?,? Then, for each metagenomic data set, the sequencing depth was estimated per position (Bowtie; BEDTools) ?,? and truncated to the central 80% values (i.e., removing the top and bottom 10% positions by depth; BedGraph.tad.rb),? a metric hereafter termed TAD80 (truncated average sequencing depth). Relative abundances were then calculated as TAD80 normalized by the genome equivalents of the metagenomic data set (MicrobeCensus),? resulting in units of community fraction. For the taxonomic classification of the genomospecies, the GTDB-tk v2.4.0,? which is based on the alignments of 120 single-copy tagged genes, was used with the representative MAGs. Phylogenomic trees were constructed using PhyloPhlAn 3.0? (function supertree_aa.cfg), and the corresponding amino acid identity (AAI) matrixes obtained (function supermatrix_aa.cfg). The final trees were visualized and drawn with iTOL.?

Finally, CRISPRCasFinder was implemented to detect CRISPR-Cas systems in the genomospecies, keeping only hits with levels 3 and 4 of confidence.?

Detection and Analyses of Viral Contigs

Putative viral contigs ≥3 kb were extracted from the metagenomes with Vibrant? and VirSorter2.? Then, putative viral contigs were analyzed using the CheckV tool,? which give a quality classification in addition to an estimation of the number of genes associated to phages or cells. Based on the obtained results, only contigs with more viral genes than host genes, and within the quality classification of “complete”, “high-quality”, “medium-quality” and “low-quality” were considered as “putative viral” contigs. Also, the viral genes of each putative viral contig had to represent at least 8% of their total predicted genes.?

Putative viral contigs which fulfilled these criteria were considered as viral contigs bona fide and were clustered into viral OTUs (vOTUs) based on their average nucleotide identity (ANI, ≥95%) and the alignment fraction (AF, ≥85%).? The vOTUs were functionally annotated as described for contigs or MAGs, and taxonomically classified using the geNomad tool.? The presence and relative abundances of the vOTUs in the different DWDS samples were assessed through read mapping of the trimmed reads from each metagenome by BLASTn, with a cutoff of 70% query coverage, e-value < 10^–1^ and the “best hit” option. Only those vOTUs mapped with identities ≥95% along ≥75% of the vOTUs reference sequence were considered.? The in silico putative host assignation was conducted by the iPhoP tool,? as well as with the protospacers of those CRISPR-Cas systems detected in the MAGs of this study against the vOTUs by BLASTn.?

Viral OTUs assigned to Mycobacterium were compared by BLASTn against all of the Mycobacterium isolated phages present in the NCBI Viral Genomes Resource. In addition, vConTACT2? was implemented to assess the relationships between the Mycobacterium-associated vOTUs and all the Mycobacterium isolated phages. The tool utilizes a gene content-based network approach to cluster viral genomes into related groups based on shared protein families. The resulting viral genome similarity network was visualized using Cytoscape (version 3.9.1).?

Statistical Analyses

Statistical differences in physicochemical characteristics between the different water samples were calculated with PERMANOVA? (permutations: 4999, method: “bray”) after checking the homogeneity of variance with “betadisper” from the Vegan Package? in R v.4.0.3, while statistical differences in Alpha diversity parameters (Richness and Shannon’s index) in the whole data set were calculated using the Stats Package in R and plotted using ggplot2.? Normality was examined using the Shapiro-Wilk test, and nonparametric tests (Mann–Whitney U Test and Kruskal–Wallis Test) were applied to the Alpha diversity values of the different sample types (biofilm or water).

Statistically significant differences of the OTUs distribution in waters and biofilms were also evaluated using PERMANOVA? (permutations: 4999, method: “bray”) and represented using a using nonmetric multidimensional scaling (NMDS) based on a Bray–Curtis distance calculation with the Vegan Package? in R. SIMPER analysis, included in the Vegan Package,? was used to determine “prevalent” bacteria within the data sets, those with a relative abundance statistically higher (p-value < 0.05) under a certain condition. As for metabolic traits differences, the SIMPER method was also employed to identify any of the detected functions from the various databases used for the metagenomes and MAGs annotation that showed significant differences among the conditions (sample type, biofilm stage).

Results and Discussion

Sampling Characteristics and Physicochemical Properties

Seven EPDM rubber hoses (A–G) that connected the main DWDS from six buildings to individual washbasins were sampled (FigureA,B). While the time of use of segments A–F could not be determined (not even by considering certain biofilm characteristics such as the amount of recovered biomass, because of inorganic deposits mixed with microbes), segment G was sampled one year after its installation. The chlorinated water circulating through these devices showed no statistically significant physicochemical differences (Table S3), including the water sample taken one year later. Cell and VLP concentrations were in the range of 10^4^/mL (Table S4), consistent with previous studies.?

Biofilms with the highest microbial biomass (D and F) were visualized under SEM (FigureC–F) and were also used to chemically characterize their EPS matrix. No sugars could be identified by Mass Spectrometry, and although various peaks were observed by HPAEC-PAD, only glycerol could be clearly detected (Figure S1). Albeit glycerol itself is not a backbone component of exopolysaccharides, it can be involved in the biofilm matrix as a metabolic byproduct or as part of glycerol-containing lipids, such as the Mycobacterium mycolic acids.?

Metataxonomic Study of the Selected DWDS: Biofilms and Planktonic

Fractions

Amplicons of 16S rRNA genes were obtained for all water and biofilm samples, yielding 1,384,422 high-quality reads which were clustered into 1,321 OTUs at 99% identity (to refer to the amplicon sequences of each sample, prefixes “BA”, biofilm amplicon, or “WA”, water amplicon, were added before the sample names). Rarefaction curves reached a plateau in all cases, indicating that the sequencing effort captured most of their microbial diversity (Figure S2), and Good’s coverage, which calculates the proportion of nonsingletons,? was above 99.8%.

Alpha diversity was computed using the number of different OTUs (richness) and the Shannon diversity index (FigureA,B), with no statistically significant differences between biofilms and water samples. Shannon’s indexes were also used to determine the degree of biofilm maturity (FigureB), as in previous works, ?,? since no other characteristics of the samples could be used. As expected, the lowest Shannon value corresponded to biofilm G, whose biomass was collected one year after the EPDM hose placement. Consequently, we catalogued biofilm G as “incipient”, biofilms A, C and D as “intermediate” and B, E and F as “old” biofilms. Regarding beta diversity, an NMDS ordination based on OTU relative abundances of the OTU showed that most samples clustered according to their origin (biofilm or water) (FigureC). This separation was further supported by PERMANOVA analysis (p = 0.019). Indeed, only 24.4% of the OTUs were detected in both sample types. At the genus level, a higher proportion of OTUs (43%) belonged to genera shared between water and biofilm samples (Tables S5 and S6). However, these shared genera accounted, on average, for up to 65% of the reads in biofilm samples but only 26% in water samples. Given that no significant differences in the physicochemical parameters of the water were observed among the studied DWDS samples, even one year later, the observed differences in beta diversity between waters and biofilms might be explained by the presence of planktonic bacteria that do not participate in biofilm formation as well as bacteria detached from biofilms developed upstream on different materials.

(A,B) Alpha diversity metrics (y-axis) derived from the 16S rRNA gene metabarcoding analysis from biofilm and water samples (x-axis). (C) Nonmetric multidimensional scaling (NMDS) analysis based on the relative abundances of the OTUs from biofilm and water samples. Before each sample name, the preffix “BA_” (biofilm amplicon) or “WA_” (water amplicon) was added.

Overall, water and biofilm samples were dominated by members of Bacteria (Table S5), with Pseudomonadota, Actinomycetota, Cyanobacteriota, Bacteroidota and Planctomycetota as the phyla with the highest averaged relative abundances (FigureA). Statistical analyses were performed to distinguish between “generalist” bacterial genera in the DWDS (present in ≥ 50% of both biofilm and water samples with no significant differences in their relative abundances between both types of samples), “specific” genera (present in ≥ 50% of biofilm or water samples), or “prevalent” genera (those with statistically higher relative abundances in biofilms or water) (FigureB; Table S6). Only Mycobacterium, Chryseolinea and Amphiplicatus, all of them previously detected in other water facilities, ?,?,? fell into the category of “generalist”. They were also “abundant”, since their relative abundances were ≥1% (in an analogous way to the definition of abundant biosphere by Pedrós-Alió?) in ≥50% of the samples (FigureB). Mycobacterium, able to survive and develop in treated waters and considered as an emerging pathogen of concern for water utilities,? was not detected in the incipient biofilm G, but it was present in all water samples and in the rest of biofilms with high relative abundances. To our knowledge, there are no studies describing Mycobacterium as a surface colonizer; however, mycobacteria are known to promote biofilm formation and to dominate established biofilms. ?,?

Heatmap of the bacterial phyla (A) and genera (B) that exceeded 1% relative abundance in at least one biofilm or water sample. Those genera present in ≥50% of both biofilms and water samples, with no significant differences in their relative abundances between both types of samples, were considered as “generalist”, while those genera present in biofilms and water samples, with a relative abundance statistically higher in waters were considered as “prevalent in water”. Genera exclusively present in ≥50% water samples were considered as “specific in water”. Genera “prevalent in biofilms” or “specific in biofilms” are not shown as none of them exceeded 1% of relative abundance in any sample. Before each sample name, the prefix “BA_” (biofilm amplicon) or “WA_” (water amplicon) was added.

In water samples, 5 and 9 genera were specific or prevalent, respectively, and commonly found in aquatic environments, including DWDS. ?,? Among the specific genera, only the cyanobacterium Candidatus obscuribacter exceeded 1% relative abundance in at least one sample (FigureB). In biofilms, 3 genera were specific and 2 were prevalent; in terms of relative abundances, all of them showed values < 1%. However, in the incipient biofilm G high relative abundances of some generalist genera were found, such as Erythrobacter, Sphingomonas, Aquabacterium and Blastomonas, which include biofilm-forming species (FigureB). ?−? ? ? As the biofilm’s maturity increased, 5 generalist genera (Amphiplicatus, Dongia, Gemmata, Ketobacter and Piscinibacter) increased their relative abundances, exceeding 1% in the oldest biofilms. Remarkably, the genus Pseudomonas, ubiquitous in DWDS, ?,?,? was only detected in two biofilms with <0.1% of relative abundance. Its low relative abundance could be due to the pipe material since Pseudomonas biofilms are preferably developed on ductile iron pipes rather than polyvinyl chloride, polyethylene or steel materials. ?,?

In Situ Distribution of Phyla in DWDS Biofilms

CARD-FISH was used on the biofilms with the highest amount of biomass (B-F) to obtain information about the spatial organization of the main bacterial phyla, in terms of relative abundance, within the biofilm matrix (Figure). Our results showed no specific phyla aggregation or distribution, but a general presence all over the biofilms suface. Detected bacteria were estimated from around 3.6 × 10^3^ cells/mm^2^ (for Actinomycetota and Bacteroidota) to 1.0 × 10^3^ cells/mm^2^ (for Planctomycetota and Cyanobacteria). Members of the classes Alphaproteobacteria and Gammaproteobacteria were observed with a frequency around 2.0 × 10^3^ and 1.2 × 10^3^ cells/mm^2^, respectively.

CARD-FISH images of those phyla with the highest relative abundances in biofilms. Each image is the result of the composite of those cells hybridized with a given probe (in blue, white arrows) and the rest of cells stained with SYBRTM Gold (in gray).

Functionality in the Selected DWDS through Metagenomics

The microbial functional patterns from water and biofilm communities were analyzed in all the biofilms and 3 water samples by metagenomics (DNA amounts extracted from water samples A, D, E and G were not enough for sequencing). All metagenomes were of high quality, based on Nonpareil data (Table S7) and the taxonomic composition of the metagenomic reads, as well as the diversity parameters, were consistent with those values obtained in the metataxonomic analysis (Table S8).

Open reading frames (ORFs) were predicted from contigs ≥1 kb and protein sequences were annotated with the KEGG database? (Table S9). The percentage of genes associated with KEGG Orthologs in each metagenome ranged from 40.16% to 48.69%. No statistically significant differences were observed between biofilm and water samples at the different KEGG levels; however, the number of metabolic pathways reconstructed was higher in biofilms (211) than in water (187). For the three biofilm stages that number was in concordance with their degree of maturity (86 pathways in the incipient biofilm, 194 in intermediate biofilms, and 205 in old biofilms) (Table S9).

To obtain information about biofilm dynamics, we searched for genes related to biofilm evolution and virulence, such as those involved in queuosine synthesis, EPS metabolism, and ARGs (Table S9). Queuosine is a modified nucleoside which is incorporated into tRNAs to enhance the expression of Q-genes, those enriched in NAU codons and that control various processes in bacteria, such as cell adhesion, biofilm formation or virulence.? In the studied DWDS, genes for the synthesis of queuosine were detected at significantly higher abundances in biofilms than in the water samples and also in the intermediate/old biofilms compared to the incipient one (Table S9). With respect to the EPS, the biosynthesis of complex exopolysaccharides involves the action of different glycosyltransferases (GTs).? Genes encoding GTs from the families GT2 and GT4, predicted to act in the synthesis of specific polysaccharides such as those contained in colanic acid, alginate, Psl, Pel, PNAG, xanthan, succinoglycan, and cellulose, were found in the metagenomes of all biofilms, with relative abundances which were proportional to the biofilm maturity (Figure S3; Table S9). However, since these CAZymes families are not specific enough, we searched for the genes involved in the synthesis of the abovementioned exopolysaccharides. In the incipient biofilm, no complete pathways could be reconstructed (Figure S3), consistent with the fact that these functions become more complex as biofilm gains in microbial diversity. ?,? KOs for the synthesis of Mycobacterium mycolic acids, important in the formation of hydrophobic extracellular matrices,? were detected in all intermediate and old biofilms, where Mycobacterium was abundant (Table S9).

Antibiotic resistance genes (ARGs) were also detected in all metagenomes, except in the incipient biofilm, with their highest proportions in old biofilms (Figure S3; Table S9). Our observation was consistent with the literature? and highlights the importance of frequent DWDS treatments, such as Ice Pigging, to mitigate the biofilm development which can converge toward an increase in the proliferation and dispersion of antibiotic resistance mechanisms. Indeed, although the Drinking Water Directive does not set specific parameters for antibiotic resistance bacteria (ARB) or ARGs as contaminants with limit values, it does introduce new requirements that are relevant to the fight against antimicrobial resistance. Specifically, the Directive “adopts a One Health approach where all the emerging contaminants should be analyzed, including ARB or ARGs”.

Among the detected ARGs, the most abundant was the gene adeF, which is a membrane protein of the AdeFGH multidrug efflux complex. This complex is related with resistance against tetracyclines, among other antibiotics.? However, we must be cautious since the complete operon was not detected and, in any case, the role of the AdeFGH complex in antibiotic resistance seems to depend on its level of expression, rather than in its mere presence.? Other relevant detected ARGs were those coding for the OXA-50 protein (a beta-lactamase)? and an erythromycin-resistance rRNA methylase.? In the old biofilms, complete pathways for the synthesis of kanosamine (an antibiotic which can also be part of more complex antibiotics such as kanamycin?) were detected, as well as genes for the synthesis of pyrrolnitrin (a secondary metabolite with antifungal and antibacterial properties?) (Table S9).

Metagenome-Assembled Genomes: “Who” Does “What”

Community metagenomics does not easily assign functions to specific taxa, so, to unveil the role of certain microbes within the DWDS, a total of 185 high- and medium-quality MAGs (≥80% completeness, < 5% contamination) were reconstructed and clustered in 104 genomospecies (ANI ≥ 95%) (Table S10). All genomospecies accounted for 30.9–61% of the reads in each metagenome (FigureA) and represented 22.4% of the OTUs obtained in the metataxonomic study using an identity threshold above 98.6%? (Table S10), so they were considered as good representatives of the total genomic diversity of the samples. For biofilm G and the water sample B, with the lowest diversities, most reads were recruited by only one genomospecies. Indeed, only one MAG could be assembled from incipient biofilm G, whose reads did not map against any other genomospecies.

(A) Cumulative relative abundance (TAD80/GE; y-axis) of bacterial phyla in the DWDS metagenomes (x-axis) based on the genomospecies composition. (B) Phylogenomic tree of the obtained 104 genomospecies. In (B), the layers over the genomospecies names are a heatmap of the relative abundance of each genomospecies in the different biofilm stages (incipient, intermediate, and old) and the water samples. For the relative abundances in intermediate and old biofilms, average value were calculated with the data of the 3 metagenomes of each stage. Red dots indicate the presence of, at least, one ARGs in each genomospecies. Branch colors are associated with genomospecies considered as “biofilm producers” (in brown) or the main “biofilm degraders” (in green; those with ≥1% of their genes associated with EPS degradation).

MAGs were classified in 12 bacterial phyla (FigureB). Remarkably, 95% of the genomospecies belonged to unclassified taxa below the family level, indicating that the DWDS is a source of taxonomic novelty. Six genomospecies, affiliated to Mycobacterium gadium (Gsp09), genus Bradyrhizobium (Gsp71) and families Parvularculaceae (Gsp62), Sphingomonadaceae (Gsp76), Obscuribacteraceae (Gsp25) and Hyphomonadaceae (Gsp56), were ubiquitous and abundant in the system, as also determined in the metataxonomic study (FigureB).

Seventy-eight percent of the genomospecies harbored at least one ARG. The resistances against a broader range of antibiotics were present in members of Cyanobacteriota, Actinomycetota, Acidobacteriota and Gammaproteobacteria, while some taxonomic groups, such as Gemmatimonadota, the UBA11222 order within Alphaproteobacteria or the Sphingorhabdus and Immundisolibacter genera, among others (FigureB; Table S11), did not present any ARG (or no contigs with ARGs were binned for these MAGs). As previously mentioned, adeF was the most abundant ARG, while an RNA-polymerase binding protein (RbpA), which confers resistance to rifampin,? and the Mycobacterium tuberculosis gene murA, that confers intrinsic resistance to fosfomycin,? were found in the genomospecies affiliated to Mycobacterium. The wide distribution of ARGs among our genomospecies is consistent with the correlation between chlorine pressure and antibiotic resistance emergence.?

Genes for queuosine synthesis were detected in 65% of the genomospecies (Q-sources), mainly members from Pseudomonadota, and some genomospecies were associated with the biosynthesis of antibacterial compounds or with the degradation of some potential contaminants. The syntheses of beta-lactams (cephamycin C), fosfomycin, and pyrrolnitrin were the most represented pathways, although incomplete, within members from Pseudomonadota, Acidobacteriota, Actinomycetota and Bacteroidota. Remarkably, 34 genomospecies from 8 phyla, well distributed among the different biofilm stages, presented genes merA and/or merB, which participate in the methylmercury detoxification (Table S11).? Also, complete metabolic pathways for the degradation of other contaminants that can impact the quality of the water were found, such as catechol, phenylacetate, and phthalate, among others. Phthalates are widely recognized as endocrine-disrupting compounds? and phenylacetate is a chlorination byproduct which is associated with undesirable effects.? Catechol is not commonly studied as a direct contaminant in DWDS, but it can contribute to biofilm formation and act as a precursor for various disinfection byproducts.? Such degradation pathways were generally found in genomospecies, which were more abundant in old biofilms and classified within Pseudomonadota, Actinomycetota and Gemmatimonadota.

All genomospecies seemed to be involved in the evolution of the EPS matrix. Those genomospecies with more ORFs related to EPS formation rather than EPS degradation were considered as “biofilm producers” (22% of the genomospecies), while those showing the opposite trend (72%), were considered as “biofilm degraders”. Biofilm producers mainly affiliated with Pseudomonadota and Cyanobacteriota (FigureB; Figure S4) and harbored at least one gene related to the synthesis of alginate, Pel, Psl, succinoglycan or colonic acid and, as expected, KOs for mycolic acid production were detected in the three Mycobacterium genomospecies (Table S11). Classes Bacteroidia and Planctomycetes, well-known polysaccharides degraders, ?,? contained genomospecies with ≥1% of their genes associated with EPS degradation (FigureB), with Gsp20 (order Cytophagales) harboring the highest proportion of genes associated with polysaccharides degradation (2.3%) and half of them directly related to the degradation of alginate, Pel and Psl. In general, most “biofilm degraders” and “biofilm producers” increased their abundances along with biofilm maturity.

In summary, the formation and evolution of biofilms are associated with heterogeneous microbial communities, which pose a challenge to the routine maintenance of DWDS, as they contribute to water quality deterioration, increased disinfectant demand, and health risks. The presence of potential pathogens and ARGs may represent suitable targets for the development of molecular diagnostic tools to support rapid decision-making within Water Safety Plans (WSPs; Environment, Climate Change and Health [ECH], 2022). Additionally, strategies aimed at preventing biofilm formation and thereby reducing the release of microbial cells into the water flow are critical for maintaining water quality.

Insights into the DWDS Viral Population Diversity and Dynamics

The diversity and impact of viruses in DWDS remains underexplored despite their role either killing biofilm microbial producers and potentially pathogenic and antimicrobial resistant bacteria, or disrupting the biofilm matrix to better expose microbes to the effect of viruses or disinfectants.?

In this study, more than 1,000 viral contigs extracted from water and biofilm metagenomes were clustered in 712 vOTUs, with more than 80% of the viral contigs from biofilms and an increase in their recovery as biofilms mature (Table S12). The vOTUs recruited up to 6.11% of biofilms reads and most of them (96%) were assigned to the class Caudoviricetes (Figure S5), which contain tailed viruses. None of them could be classified at the genus/species level, underscoring the taxonomic novelty of the DWDS viral assemblages.? Considering all of the reads that mapped against the vOTUs (the set of detected viral reads within each metagenome), Caudoviricetes recruited 98 and 99.9% of them in biofilms and water metagenomes, respectively.

The IPhoP software assigned 171 vOTUs to 55 bacterial genera detected in the studied DWDS, including Mycobacterium, Bradyrhizobium, Bosea, Sphingopyxis or Erythrobacter, among others (Table S12). In addition, the search of 1,602 CRISPR-Cas systems spacers, extracted from 29 genomospecies, against the vOTUs, allowed for the association of 2 vOTUs to Gsp85 (Alphaproteobacteria - UBA11222) and Gsp89 (Aquabacterium commune) (Figure S6). In general, the relative abundances of those vOTUs and hosts in the abundant biosphere were positively correlated (p-value < 0.05). The correlation was very clear in the case of Mycobacterium and their associated viruses. Mycobacterium dominated the microbial community in intermediate biofilms (21.9% of 16S rRNA reads in average) and was linked to vOTUs that recruited >60% of the viral reads in those samples. Half of the viral reads were recruited by vOTUs with integrases (putative prophages). Remarkably, although Mycobacterium is among the bacterial genera with the best described virosphere,? only 5 vOTUs presented a direct gene-sharing network with described mycobacteriophages, which reinforced the presence of an important viral novelty in this type of environment (Figure S7).

Since phages are also known for their roles as “EPS degraders”,? some specific genes, such as viral depolymerases, were looked for in our set of vOTUs. Putative CAZyme encoding genes were predicted in 11% of the vOTUs. In most of them (57 vOTUs) they were related to peptidoglycan degradation (GH19, GH23, GH24, GH73, GH104, GH108), while 11 vOTUs harbored enzymes that could act against alginate compounds (Table). In addition, lysB genes (Lysin B) were detected in three vOTUs, with two of them (vOTU67 and 84) assigned to Mycobacterium viruses. Lysin B is a mycolylarabinogalactan esterase that presents antimycobacterial activities and is predominantly associated with mycobacteriophages. This enzyme facilitates the host lysis by compromising the integrity of the mycobacterial outer membrane.? Also, Lysin B has shown promising activity against multidrug resistant strains of M. tuberculosis, suggesting a potential role as a therapeutic agent in combination with other treatments.? Although phage therapy remains limited by regulatory barriers, research on phage-host interactions within DWDS biofilms offers new insights into its potential effectiveness for biofilm control. The application of phages or viral enzymes is emerging as a promising approach that could be integrated into the development of modern Water Safety Plans.

1: CAZy Families Predicted among the vOTUs

Supplementary Material