Diversity and Functions of Epilithic Riverine Biofilms

This article relates epilithic dry- and wet-seasonal bacterial biofilm composition to water quality along Río de la Sabana near Acapulco, Mexico. Samples were taken from various locations including nearly pristine upland locations, adjacent to residential floodplain developments, and immediately upstream from an estuarine lagoon. Bacterial composition was identified through sequential DNA analysis at the phylum, class, order, and family levels, with most of these categorized as heterotrophs, autotrophs, denitrifiers, nitrogen fixers, pathogens, and/or potential bioremediators based on generalized literature-sourced assignments. The results were interpreted in terms of location by extent of effluent pollution, and by dry versus wet seasonal changes pertaining to biofilm composition, related bacterial functions, and the following water quality parameters: temperature, pH, dissolved oxygen, biological and chemical oxygen demand, fecal and total bacteria counts, methylene blue active substances, electrical conductivity, and nitrite, nitrate, ammonium, sulfate, and phosphate concentrations. It was found that epilithic bacterial biofilm diversity was richest during the wet season, was more varied in abundance along the upland locations, and was dominated by Proteobacteria and Bacteroidetes with bioremediation and pathogen functions along effluent-receiving river locations. Low-abundance families associated with anaerobic and denitrifying functions were more prevalent during the wet season, while low-abundance families associated with aerobic, N2-fixing and pH-elevating functions were more prevalent during the dry season.


Introduction
This article reports on epilithic bacterial biofilm composition, functions, and changes along Río de la Sabana as it flows north to south east of Acapulco, Mexico. This river and its floodplain provide high biodiversity habitats, but these are now affected by rapidly evolving urban encroachments (Arriaga-Cabrera et al. 2008). As previously reported, these encroachments have not only accelerated environmental degradation and water pollution along the river and its coastal lagoon (CONAGUA 2012;Pineda-Mora et al. 2018), but also increased floodinduced health vulnerabilities in low-lying settlements across and adjacent to the Río de la Sabana Valley (CONAGUA 2012;Rodríguez-Herrera et al. 2013).
In general, biofilms are densely populated and self-stabilizing microbiomes (Besemer 2015;Peipoch et al. 2015;Findlay and Battin 2016;Nicholls and Crompton 2017) involving algae, diatoms, fungi, bacteria, and protozoa, all enveloped in an extracellular polymeric substance matrix (EPS) on biotic and/or abiotic surfaces (Nadell et al. 2009;Villeneuve et al. 2013;Nadell et al. 2016). As such, biofilms undoubtedly affect a variety of processes and functions by way of organic matter production and decomposition, nutrient uptake and cycling, and contaminant absorption and processing (Findlay and Battin 2016;Shahsavari et al. 2017;Guasch et al. 2017;Banerjee et al. 2018). In addition, several studies have stressed that microbial biofilm compositions are subject to physical, chemical, and biological changes in water flow and quality (Pesce et al. 2010;Montuelle et al. 2010;Tlili et al. 2011;Baek et al. 2013;Langenheder et al. 2016;Bouchez et al. 2016;Hu et al. 2017). These changes can now be revealed in terms of operational taxonomic units (OTUs) using next generation sequencing (NGS) DNA technologies (Ghiglione et al. 2014;Tan et al. 2015;Pawlowski et al. 2018). With this background, this article reports on changes in bacterial biofilm composition, diversity, and functions as affected by the extent of residential effluent pollution and related seasonal influences along Río de la Sabana. The emphasis is on epilithic biofilms because their bacterial composition would directly be affected by locational changes in water quality, varying flow rates, and desiccating exposure incidences.

Study Area and Field Sampling
Río de la Sabana begins in the north-east portion of the Acapulco municipality, with its headwaters reaching 2000 m above sea level. Its main flow channel is 30.7 km long and widens to about 60 m at its confluence with the coastal Tres Palos lagoon (CONAGUA 2012). Its 466 km 2 basin area supports permanent flow channels with a cumulative length of 727 km, and a channel density of 1.55 km −1 (Villegas-Romero et al. 2009). The steep and mostly pristine terrain of the upper basin feeds the main flow channels with short-duration runoff. In contrast, the lower part of the terrain widens into the broad and densely populated La Sabana Valley floodplain (Pineda-Mora et al. 2018).
This study was carried out in 2017 during the dry (January) and the rainy (July) season, with samples taken from the basin headwaters to the floodplain estuary ( Fig. 1, Table 1) as per the following sampling zonation (Dorigo et al. 2008): 1. Reference Zone (Zone 1; Locations S1, S2): located on the upper part of the sub-basin, meeting national and international standards for the protection of aquatic life, representing oligotrophic water quality conditions 2. Transition Zone (Zone 2; Locations S3, S4): a periurban zone located between the upper and mid-low parts of the basin. 3. Pollution Zone (Zone 3; Locations S5, S6): representing a densely populated area and eutrophic water quality conditions. 4. Recovery Zone (Zone 4; Location S7): a confluence area that combines water from densely populated area west of the river, with water also received from the low population area in the east.
In total, 14 composite biofilm samples were gathered from submerged rock surfaces (Fig. 2): 7 per dry season, and 7 per rainy season using the point contact method (Lyautey et al. 2010). At each location, three rocks of similar type and surface area (about 15 cm in diameter) were selected from a transverse profile approximately 0.3, 0.5, and 0.7 m apart. At least 30 g were collected per sample through scraping the rocks with sterile disposable brushes and using sterile Falcon tubes for collection (50 mL). These samples were stored on ice until ready for centrifugation at 6,000 rpm for 10 min to eliminate excess water. Thereafter, the samples were stored at − 20°C until processed for DNA extraction.

DNA Extraction, Sequencing, and Phylogenetic Analysis
Deoxyribonucleic acid (DNA) was extracted from the biofilms according to the methodology proposed by Ceja-Navarro et al. (2010) and Navarro-Noya et al. ( 2013).

Water Quality Sampling and Analysis
As already reported by Pineda-Mora et al. (2018), water quality was sampled at the same locations during the dry season (January) and rainy season (July) in 2017 (Table 1). Fourteen physicochemical and bacteriological water quality parameters (Table 2, Supplementary Data IV) were measured following established protocols in accordance with Mexican references (NMX) and standardized alternative methods. Temperature, pH, and electrical conductivity were measured in situ. The other parameters were analyzed in a laboratory authorized through Mexican Organization of Accreditation (AC). Temperature, pH, electrical conductivity, ammonium nitrogen, nitrates, nitrites, methylene blue active substances (anionic surfactants), sulfates, and phosphates parameters were measured in triplicate. Dissolved oxygen, biochemical oxygen demand , chemical oxygen demand, total coliforms, and fecal coliforms were determined by way of single measurements per sample.

Statistical Analyses
The bacterial relative abundance data, identified by order, class, family, and unclassified, were compiled into spreadsheets. The BiodiversityR package (Kindt and Kindt 2019) in R studio was used to determine microbial diversity indices by location and season: diversity by type (Richness, Kindt and Kindt 2019), by type and abundance (Simpson 1949), and by diversity of dominants (Berger and Parker 1970).

Location 2 (34 km uplands)
Location 5 (Pte. Alborada) Fig. 2 Examples of biofilm-covered rocks lifted from Location 2 representing oligotrophic upland water conditions (top), and from Location 5 representing eutrophic water conditions near the residential floodplain developments (bottom). Biofilms were found to be thicker, less translucent, and greener along Locations S4 to S7 The sums of the relative taxonomic abundance frequencies contributing to bacterial composition, functions, biodiversity indices, and water quality parameters were all compiled by location and season. Locations S1, S2. S3, S4, S5, S6, and S7 where coded 1, 2, 3, 4, 7, 6, and 5, respectively, to reflect the overall extent of residential effluent pollution. Season was coded 0 when dry and 1 when wet. The resulting variable-to-variable association patterns of these sums were (i) factor analyzed (Principal components, 2 factors, Varimax transformation) and (ii) were subjected to multivariate regression analyses to quantify the extent to which bacterial relative OTU abundances, functions, and water quality parameters relate to each other by location and season. The resulting two-factor analyses plots were sectorized by extent of pollution and season so that Sectors I and II refer to least-polluted (S1, S2, S3) and Sectors III and IV refer to most polluted (S4, S5, S6, S7) locations during the wet and dry seasons, respectively.

Bacterial Biofilm Abundances, Phylum Level
The > 1% relative abundances of the DNA-recognized bacterial phylum (p), class (c), order (o), or family (f) in the Río de la Sabana biofilms are presented in Supplementary Data I. Altogether, 11 entries were identified at the phylum level, 27 at the class level, and 68 entries at the family/order level. Of the ≥ 1% entries, 27 were present during the dry and wet seasons, 20 occurred during the dry season only, and 24 occurred during the wet season only.
The bar-charted phylum abundance display in Fig. 3 (i) across all locations and seasons combined, (ii) across all locations by season, and (iii) across the least to most polluted locations by season revealed the following differences: increased pollution promoted the relative abundance on Proteobacteria, Bacteroidetes, Chloroflexi, at the S4-S7 locations but decreased the relative abundances of Actinobacteria, Firmicutes, Cyanobacteria, and Thermi. During the wet season, Proteobacteria and Bacteroidetes further increased their relative dominance at the S4-S7 locations at the expense of Firmicutes and Chloroflexi.
These compositional phylum changes are further summarized in Fig. 4 by the 2-factor analysis plot that differentiates the relative abundance sums of each phylum as they fall into Sector I (Locations S1-S3; least polluted, wet season), II (Locations S1-S3; least polluted, dry season), III (Locations S4-S7; wet, most polluted), and IV (Locations S4-S7; dry, most polluted). Proteobacteria, Chloroflexi, and Synergistetes clustered along the polluted Sector III to IV transitions. Acidobacteria, Actinobacteria, and Planctomycetes dry & wet season S1-S3 S4-S7 S1-S3 S4-S7 dry season wet season Relative phylum abundances % S1-S7 dry season S1-S7 wet season Fig. 3 Relative abundance of biofilm bacteria grouped (i) across locations and seasons, (ii) by season across locations, (iii) from least to most polluted locations, by season, i.e., S1 to S3 versus S4 to S7    and Rhodobacteraceae were present at ≥ 1% at all seven locations and both seasons (Table 3). DNA sequences associated with Moraxellaceae ≥ 1% were absent at S2 during the wet season, and absent at S1, S2, S3 during the dry season.
Most of the other families occurred at < 10% abundances, but their relative abundances were strongly affected by pollution and season, and somewhat related to their sector I to IV phylum assignments as apparent from Fig.4 and Tables 3 and 4. For example, classes and families associated with Acidobacteria, Cyanobacteria, and Planctomycetes occurred for the most part in Sector I and II (least polluted during the wet or dry seasons, respectively). In detail, the phototrophic Acaryochloridaceae, Chamasiphomnaceae, Chroococcales, Phormidiaceae, Pseudanabaenales, Stramenopiles, and Xenococcaceae members of the Cyanobacteria phylum occurred with relative abundances ≥ 1% only at the upland S1 to S4 locations (Sectors I and II), and more so during the dry than the wet season (Table 3, Fig. 5).
Factor analyzing the Tables 3 and 4 entries by extent of pollution and season revealed three clusters in the resulting 2-factor, 4-Sector plot in Fig. 6: one cluster in Sector I (least polluted, wet season), one in Sector II (least polluted, dry season), and one straddling Sectors III and IV (most polluted, with entries somewhat influenced by dry and wet season). In this, residential pollution or lack thereof also impacted the relative abundances of the bacterial OTUs. In detail, Alpha-and Betaproteobacteria occurred in all four sectors; Gammaproteobacteria were exclusively found in Sectors III and IV, while Deltaproteobacteria were considerably less abundant ( Fig. 6; Supplementary Data I). While Proteobacteria dominated in Sectors III and to IV, some were also be associated with Sectors I and II, especially soil-associated Rhizobiales, Hyphomicrobiaceae, P h y l l o b a c t e r i a c e a e , B r a d y r h i z o b i a c e a e , Oxalobacteraceae, Acetobacteraceae, Phyllobacteriaceae, and Syntrophobacteraceae.

R h i z o b i a l e s ( N 2 -f i x i n g m i c r o s y m b i o t i c
Alphaproteobacteria) were also present during the dry and wet season, but with a strong preference for the less polluted upper Río de la Sabana locations during the dry season.

Bacterial Biofilm Diversity
The bacterial diversity assessment results based on the Table 3 entries are listed in Table 5, by location and Most frequent bacterial order and family occurrences within Río de la Sabana biofilms: relative abundances averaged across S1 to S7 during the dry (left) and wet (right) season season as coded. Factor analyzing these results identified pollution and season as predictive factors for bacterial richness by type (Richness index, higher during the wet than the dry season), diversity abundance (Simpson index, increasing with decreasing pollution level), and dominance (Berger index, increasing with increasing pollution level). Expressed in regression terms, the resulting Factor 1 scores were related to pollution extent, and the Factor 2 scores related to season, as follows: Factor 1 scores related to pollution  6 Factor-analyzed pattern for correlational bacterial family and order abundance pattern, colored by phylum association, with additional reference to α, β, γ, and δ Proteobacteria orders

Bacterial Biofilm Functions
The results of summing the relative abundances of bacterial families, classes, orders, or phylum levels according to their presumed bioremediation, pathogen, aerobic, anaerobic, N 2 -fixing, and denitrifying functions (see Supplementary Data III) are presented Fig. 7, by season and location. These plots reveal higher relative abundances of bioremediators, pathogens, anaerobes, and denitrifiers at Locations S4 to S6. In contrast, aerobes and N 2 -fixers were nearly evenly distributed across the S1 to S7 locations, and registered low relative abundances during the dry and wet seasons.

Water Quality Related to Microbial Biofilm Diversity and Functions
Factor analyzing the water quality measures as listed in Supplementary Data IV in combination with the bacterial diversity results in Table 5 and the relative abundances by bacterial functions as plotted in Fig. 7 led to the following observations by way of Fig. 8: 1. High fecal coliform counts, high electrical conductivities, and elevated temperatures were-for the most part-associated with the polluted locations. 2. These locations were also associated with high pathogen and bioremediator occurrences, which-in turn-registered elevated biological and chemical oxygen demands (   Factor analysis plot combining the general association patterns between (i) relative bacterial biofilm taxa abundances pertaining to bioremediators, pathogens, aerobes, anaerobes, denitrifiers and N 2 -fixers (white white squares), (ii) water-quality parameters ( Table 2; colored dots), and (iii) the bacterial diversity indices (black circles), by Sectors I to IV, with each signifying one of the four dry and wet season combinations pertaining to low and high pollution extent, as in Fig. 4 the freshwater locations during the dry season. In contrast, the diversity index for dominance (Berger) fell into Sector III in association with increased effluent pollution and increased abundances of biofilm-resident denitrifiers, anaerobes, aerobes, bioremediators, and pathogens. These locationand season-specific diversity changes can also be discerned through inspecting the Supplementary Data I OTU bar charts, and most clearly so by focusing on the OTU class compilation. 5. Aerobes prevailed more consistently during the dry than the wet season (Sector IV), anaerobes and denitrifiers dominated during the wet season (Sector III), while N 2 -fixers prevailed in Sector II.
By way of regression analysis, BOD and COD varied by location and season as follows: where locations and seasons were coded again as per Table 5. Correspondingly, dissolved oxygen levels decreased with increased biological and chemical oxygen demand, and this decrease was more pronounced during the wet season, i.e., Note that the relative abundance values for the aerobic and anaerobic biofilm families did not enter this formulation in a significant way. This likely related to the general expectation that aerobes and anaerobes function best when and where DO remained steady at high or low levels, respectively.
In reference to seasonal water quality changes, anaerobes including denitrifiers were prevalent under neutral pH condition. In contrast, N 2 -fixers and organicmatter consuming aerobes were associated with increasing pH levels likely due to (i) converting N 2 to NH 3 and (ii) heterotrophic release of CO 2 from organic acid groups. Relating the reported pH values to season, location, and bacterial biofilm functions generated the where N 2 -fixers and bioremediators were entered by their sum of relative abundance per location and season.
In terms of electrical conductivity, which indexes overall effluent concentrations, the regression analysis produced the following result: Hence, EC levels were not only lower during the wet season (likely due to higher flow rate and subsequent dilution) but were also lower where DO levels remained high, and where pollution loads were low to absent, i.e., at Locations S1 to S3.
The best-fitted scatterplots for Eqs. 3, 5, 6, and 7 are presented in Fig. 9 to further illustrate how water quality was affected by location and season. In particular, the changes on BOD (and therefore COD as well) were categorically relatable to season, and linearly relatable to location from least to most polluted. In contrast, the changes in DO were both categorically related to location and season, while the changes in pH and EC were-for the most part-categorical by season, and clustered by location.
Remarkably, increasing pollution levels rendered Proteobacteria, Bacterioides, Firmicutes, and Chloroflexi to dominate all the other biofilm bacteria at > 98% relative abundance during the dry and the wet season. In addition, residential effluents rendered Proteobacteria and Bacterioides to become even more dominant during the wet season through diminishing the relative abundances of Firmicutes and Chloroflexi. Among the Proteobacteria, there was a further pollution-induced dry-to-wet seasonal shift towards greater Betaproteobacteria than Gammaproteobacteria dominance. This observation parallels the relative abundance of Proteobacteria in active sludge as reported by Jiang et al. (2008): Betaproteobacteria 39.8% > Gammaproteobacteria 2 0 . 1 5 % > A l p h a p r o t e o b a c t e r i a 6 . 7 9 % > Deltaproteobacteria 4.85%.
Hence, Alphaproteobacteria would dominate biofilms in low nutrient waters such as those found at the S1, S2, and S3 locations, and heterotrophic Betaproteobacteria and Gammaproteobacteria would dominate and co-dominate polluted waters such as those found at the S4, S5, S6, and S7 locations.
At the less polluted S1 to S3 locations, Actinobacteria and Planctomycetes became more prevalent from the dry to the wet season, while the opposite occurred with Fermicutes and Thermi. These changes were likely due to (i) higher detrital availabilities that would increase Actinobacteria and Planctomycetes activities during the wet season while (ii) Fermicutes and heat-resistant Thermi would be more adapted to persist on rocky surfaces during the dry season when exposed to periods of intermittent desiccation and heat stress (Marxsen et al. 2010).
The low presence of Cyanobacteria at Locations S4 to S7 was likely caused by pollution-enhanced and light-blocking green algae growth, and this would especially be so under continuous flow conditions, as reported by Sabater et al. (2016). In addition, lower Cyanobacteria abundance at the polluted locations may also be related to high molar N:P ratios, as reported by Peipoch et al. (2019). Determining whether this also Best-fitted DO  Fig. 9 Actual versus best-fitted values for log 10 BOD (Eq. 3), DO (Eq. 5), pH (Eq. 6), and EC (Eq. 7) along the seven Río de la Sabana biofilm sampling locations, by season. Dry-season sampling results are placed on gray-tone background uptake of dissolved organic matter from an agricultural stream, versus 4% from a forest stream (Graeber et al. 2019). While there are significant correlations between bacterial biofilm diversities, functions and relative phylum, class, order, family abundances by location and/or season (Figs. 4,6,8,9,10), it is difficult to quantitatively relate specific water quality measures to bacterial biofilm functions. In part, this refers to the fact that biofilms function synergistically as self-stabilizing microbiomes with their own overall input, output, and life-supporting requirements (Peipoch et al. 2015;Toyofuku et al. 2016;Zeglin 2015). To illustrate, one would have expected that the measured changes in water pH and DO would have influenced bacterial abundances at the individual phylum, order, and family levels, but only Bacteriodetes were so correlated. In part, Bacteriodetes may be more sensitive to increasing pH and DO levels (Eq. 9) because of their anaerobic requirements and inability to produce spores. In general, however, biofilm-internal pH and DO values likely did not correlate with biofilmexternal pH and DO measurements since the former would be the result of balanced microbial biofilm coexistences. For example, biofilm internal DO levels would result from combined O 2 producing and reducing activities, and biofilm internal pH levels would depend on (i) the organic acidity of biofilm structures and cellular exudates, and (ii) the mix of cellular nutrient cation and anion uptake, respectively.

Actual EC
To discern which effects on water quality were brought about by specific family functions by location and season would be difficult due to wide phenotypic variabilities per bacterial family. In this regard, Comamonadaceae represent 29 genera within an overall 16S rRNA gene sequence similarity of 93-97%. Yet, these genera differ by function as aerobic organotrophs, anaerobic denitrifiers, Fe 3+ -reducers, hydrogen oxidizers, photoautotrophs, photoheterotrophs, and fermenters (Willems 2014a). Rhodobacteraceae are even more diverse, containing 100 genera involving aerobic photoand chemo-heterotrophs, with some performing photosynthesis under anaerobic conditions. Characteristically, members of this family facilitate sulfur and carbon biogeochemical cycling, often in symbiosis with aquatic micro-and macro-organisms (Pujalte et al. 2014). Moraxellaceae include 5 genera with non-fermentative saprophytic functions, with some species recognized as p a t h o g e n s ( T e i x e i r a a n d M e r q u i o r 2 0 1 4 ) . Sphingomonadaceae, known for their important roles as bioremediators, involve 13 genera with facultative to anaerobic fermentation functions in soils, freshwater systems, sludge, rhizospheres, and phyllospheres (Glaeser and Kämpfer 2014). Rhizobiales are primarily involved with N 2 -fixation, with some forming root nodules while others dominate in tap-water biofilms, with some being pathogenic (Potter 2013;Wang et al. 2018), and some being pathogen controlling. Pathogen-controlling species are associated with Actinobacteria (Wink et al. 2017), with Acidimicrobiales, Actinomycetales, and Nocardioidaceae occurring with ≥ 1% abundances at the S1, S2, and S3 locations of this study.
To reveal further details about positive and negative bacterial interactions in addition to those surmised above along Río de la Sabana locations, it would be necessary to conduct multiple location-and seasonspecific DNA sampling studies. In addition, such

Best-fitted relative abundance sum
Relative abundance sum Proteobacteria Fig. 10 Actual versus best-fitted relative abundance sum plots for Bacteroidetes (Eq. 9) and Proteobacteria (Eq.10) by Locations 1 to 7 and dry (white square) versus wet (filled black circle) season an effort should be expanded to include other biofilm constituents, e.g., green algae and diatoms (Kouzuma and Watanabe 2015;Ramanan et al. 2016). The focus on epilithic biofilms could also be expanded to biofilms residing on other mineral surfaces (e.g., sediments) and submerged macrophytes. Also, valuable would be to analyze how dissolved and total elemental concentrations in biofilms relate to corresponding changes in water quality.

Conclusions
The above results show how bacterial biofilm composition, functions, and water quality measures varied along Río de la Sabana from the least to most polluted locations and by season. The main effect of residential effluent discharge expressed itself through biofilm thickening caused by green algae growth and a bacterial dominance shift towards heterotrophic Proteobacteria (Comamonadaceae, Moraxellaceae, Rhodobacteraceae) and Bacteroides (Weeksellaceae). This dominance shift towards common wastewater and sludge OTUs was most pronounced during the wet season and occurred on the expense of the ≥1% relative abundance of most other families. By bacterial biofilm function and water quality associations, Locations S5 and S6 registered the highest levels of bioremediators and pathogens. Water-sampled fecal coliform counts, biological and chemical oxygen demands, electrical conductivity, and NH 3 , NO 3 − , NO 2 − , PO 4 3− and SO 4 2− , concentrations were also generally highest at these locations but more so during the dry than the wet season. Dissolved oxygen levels registered highest at the least pollution locations (S1 and S2) during the wet and dry season. In contrast, pH remained near neutral during the wet season but increased to alkaline levels during the dry season. Apart from general factoranalyzed location and seasonal association patterns between bacterial composition, diversity, functions and water quality measures, correlational changes in relative abundance for any particular OTU could only be discerned for Bacteriodetes with changing pH and DO levels in sampled stream and river water.