Two different approaches of microbial community structure characterization in riverine epilithic biofilms under multiple stressors conditions: Developing molecular indicators

Microbial communities are major players in the biogeochemical processes and ecosystem functioning of river networks. Despite their importance in the ecosystem, biomonitoring tools relying on prokaryotes are still lacking. Only a few studies have employed both metabarcoding and quantitative techniques such as catalysed reported deposition fluorescence in situ hybridization (CARD‐FISH) to analyse prokaryotic communities of epilithic biofilms in river ecosystems. We intended to investigate the efficacy of both techniques in detecting changes in microbial community structure associated with environmental drivers. We report a significant correlation between the prokaryotic community composition and pH in rivers from two different geographical areas in Norway. Both CARD‐FISH and metabarcoding data were following the pattern of the environmental variables, but the main feature distinguishing the community composition was the regional difference itself. Beta‐dispersion analyses on both CARD‐FISH abundance and metabarcoding data revealed higher accuracy of metabarcoding to differentiate regions and river systems. The CARD‐FISH results showed high variability, even for samples within the same river, probably due to some unmeasured microscale ecological variability which we could not resolve. We also present a statistical method, which uses variation coefficient and overall prevalence of taxonomic groups, to detect possible biological indicators among prokaryotes using metabarcoding data. The development of new prokaryotic bioindicators would benefit from both techniques used in this study, but metabarcoding seems to be faster and more reliable than CARD‐FISH for large scale bio‐assessment.


| INTRODUC TI ON
River systems are extremely dynamic ecosystems, providing a variety of services to humans (Arthington et al., 2018). Since ancient times rivers have been exploited in several ways: drinking water sources, water for agriculture, hydropower, cooling systems for industries, recreation, etc. (Poff, 2018;Vörösmarty et al., 2005). These activities have impaired natural flow variability, causing changes to habitats and the biodiversity of adjacent areas (Davies et al., 2014;Poff et al., 1997).
Biodiversity loss is one of the biggest concerns for river ecosystems, as most rivers are being exploited and increasingly losing species (Bunn & Arthington, 2002;Das Gupta, 2008;Merritt et al., 2010). Biodiversity loss and reduced abundance are becoming evident for several biological groups, such as insects, for which recent studies have shown global changes in freshwater (Baranov et al., 2020;Hallmann et al., 2017Hallmann et al., , 2020Sánchez-Bayo & Wyckhuys, 2019).
In contrast, very little is known about how microbial diversity is influenced by anthropogenic stress. Microbial communities play a fundamental part in driving ecosystem processes and play a vital role at the base of riverine food webs (Demars et al., 2020;Robbins et al., 2017). Despite their small physical dimensions, microbes are key drivers of organic matter decomposition in fluvial ecosystems and mineralization of nutrients, making them available for higher levels of the riverine food web (Demars et al., 2020). This functional role is particularly true for the microbial communities living in epilithic biofilms with a high sediment to water phase ratio, producing a large extent of reactive surfaces in streams (Battin et al., 2008).
The role of microbial communities is today often only inferred in studies on ecosystem functioning, for example by using fine mesh bags in leaf litter decomposition studies (Gessner & Chauvet, 2002;Tiegs et al., 2019;Woodward et al., 2012). Mapping microbial diversity, and the functional traits related to key ecosystem processes, has huge potential in increasing our understanding of drivers of ecosystem functioning in the years to come (Besemer, 2015). Moreover, due to their high sensitivity to pollution and fast response to environmental changes, bacterial assemblages could complement the information provided by benthic metazoan communities as indicators of human-induced impacts, but this biological component has not yet been well explored in this regard (Caruso et al., 2016;Szabo et al., 2007). Jackson et al. (2016) strongly argued for use of sequencing of microbial communities as part of next-generation biomonitoring tools. In Europe there is currently a significant lack of prokaryotic indicators in the EU Water Framework Directive (2000/60/EC) and national and international legislations (Heiskanen et al., 2016), although all other relevant biological groups, including microalgae, are included in the assessment of freshwater status (Birk et al., 2012). Obvious reasons for the absence of prokaryotic communities in bioassessment have been the cost of taxonomic analyses and a lack of knowledge of their indicator value in terms of natural variability and human impacts. However, the easy access to high-throughput sequencing (HTS) technologies, allowing the quick taxonomic identification of bacterial assemblages, has led to several attempts to find prokaryotic indicators. Some of the most relevant studies focused on sediments of coastal areas and estuaries (Aylagas et al., 2017;Borja, 2018).
In freshwaters, a first attempt to use prokaryotes as bioindicators focused on quantitative techniques such as real-time qPCR selecting a few prokaryotic strains associated with specific chemicals and water quality parameters (Nzewi et al., 2009). Others have tried to address the lack of prokaryotic bioindicators for freshwater ecosystems, mainly using qPCR to quantitatively analyse specific functional genes (Thompson et al., 2016) and 16S amplicon sequencing to point out changes in the community structure of freshwater prokaryotes impacted by human activities (Salis et al., 2017;Simonin et al., 2019). Martínez-Santos et al. (2018) used both qPCR and 16S amplicon sequencing to analyse the effects of wastewater effluents, on structure and function of the prokaryotic communities dwelling in Deba river sediments. In our study, we focused on epilithic biofilm communities which, in addition to being of key importance for ecosystem functioning in rivers, have been shown to be more sensitive to water quality features compared to those dwelling on high organic matter (OM)-loaded substrates (leaves, roots, wood) (Fazi et al., 2005). Hence, they could provide reliable information regarding human and natural pressures on the environment. We explored the potential of two very different techniques as biodiversity indicators for prokaryotes in epilithic biofilms: a quantitative method, catalysed reported deposition fluorescence in situ.
Hybridization (CARD-FISH; Pernthaler et al., 2002); and a more qualitative method, metabarcoding of 16S rRNA amplicons (Johnson et al., 2019). The two techniques were used to characterize the microbial community structure of rivers in Norway, ranging from river systems affected by natural disturbances, to rivers affected by various types of human impacts, such as dams, hydropower developments and wastewater outlets. Two different geographical regions were selected, characterized by distinct geological and chemical features (Steinnes et al., 1993), probably influencing the prokaryotic community structure of epilithic biofilms.
We hypothesized that: (i) CARD-FISH and metabarcoding will provide similar patterns regarding the overall microbial community structure, but they will give different insights at different spatial scales (regional vs. microscale); (ii) that the community structure of epilithic biofilms would be influenced by both human perturbation and natural conditions such as the geological setting.
The advantage of metabarcoding to CARD-FISH is lower operating costs, which might make it better suited for use in modern biomonitoring networks if the methods yield comparable results in terms of describing prokaryotic communities. There is already evidence for the effectiveness of metabarcoding as a biomonitoring tool, for both eukaryotes and prokaryotes, with results comparable to the traditional methods based on species morphology (Cordier et al., 2019). Metabarcoding of prokaryotic communities would be complementary to the traditionally used bioindicators such as benthic macroinvertebrates, diatoms and fish, which today are the most commonly used group in impact assessment (e.g Birk et al., 2012;Friberg, 2014). It would provide valuable insights into the black box of biodiversity in riverine ecosystems, namely the microbial communities which may be pivotal as early warning indicators of human and natural pressures (Besemer, 2015;Widder et al., 2014).

| Study design and samplings
Two different geographical regions in Norway were selected for the present study, characterized by distinct geological and chemical features, in particular in terms of acid neutralizing capacity (Steinnes et al., 1993) (Figure 1). A total of 16 sites was sampled, embracing an array of environmental gradients (acidity of water, geology, human impacts). In one area (in the Oslo region in southeast Norway), bedrock of the rivers is dominated by lower Palaeozoic sedimentary rocks (limestone and shales) (Calner, 2013). In the other area (Vest and Aust-Agder, southwest Norway) the bedrock mainly comprises magmatic and metamorphic rocks (Slagstad et al., 2018). The 16 sampling sites were situated in four catchments, two in each region: Lysakerelva and Sandvikselva in the Oslo region; and Arendalvassdraget and Mandalselva in the Agder region. Within these four catchments we selected four rivers impacted by hydropower and dams (Lysaker, Lomma, Nidelva and Mandal) and four free-flowing rivers (Iselva, Finnsåna, Haugedøl and Stigvasselva) which were our control sites (Table S1). For each control site, we selected only one sampling site, while for the impacted rivers we selected three sampling sites, one upstream from the dam or hydropower plant, so that this site could be comparable to our control sites as being virtually unimpacted. The second site was always set immediately below the dam, where we expected impact to be highest.
The third site was further downstream, where the water was mixed, and the effect of the dam was not evident.
When establishing a sampling site, it was georeferenced using a Water samples were collected at each site approximately at 10 cm below the surface. A 500-ml water sample was taken at each site for chemical analysis by using polyethylene bottles. Contextually another 60-ml water sample was taken for the analysis of metals by using polyethylene bottles pretreated with a 1% HNO 3 solution.
The samples were immediately placed at 4°C in the field and brought back to the laboratory for analysis.
Biofilm samples were collected at each site after taking the water sample. Five rocks (average individual surface 118 ± 22.4 cm 2 ) were randomly taken within a 50-m 2 area and brushed in the field with a sterile toothbrush to collect the epilithic biofilm. The biofilm brushed from the five rocks was pooled together to give a single sample for each site (16 biofilm samples in total). The pooled biofilm samples were suspended in 65 ml of ultrapure MilliQ water. An aliquot of each sample (50 ml), to be used for hybridization in situ (CARD-FISH), was added to 50 ml of pure ethanol to prevent ice formation and consequent cell lysis; the remaining 15 ml, to be used for DNA extraction, was placed in a 15-ml Falcon tube. We had 16 subsamples for CARD-FISH in total and 16 subsamples for sequencing in total.
Both biofilm samples were kept cool at 4°C until arrival ae the laboratory where they were stored at −20°C until further processing.

| Physicochemical water parameters
Water temperature, pH and electrical conductivity (EC) were measured in situ with a multiparameter portable meter (WTW ProfiLine Multi 3320).

| Biofilm biomass quantification
For biomass quantification, we used the ash free dry mass (AFDW) content of the biofilm samples. Two replicates subaliquots (~1 g wet weight) were taken from each sample and preserved for hybridization in situ. The subaliquots were dried at 60°C in a thermostatic oven for 72 hr to obtain the dry weights. Subsequently the dried aliquots were pooled together and burned in a muffle oven at 550°C for 3 hr to obtain the ash weights. Subtracting the ash weights from the pooled dry weights we were able to measure the biomass content of our samples.

| Total prokaryotic abundance and single cell hybridization (CARD-FISH)
The total prokaryotic cell abundance was assessed by 4′-6-diamidino-2-phenylindole (DAPI, Vector Laboratories) staining, following extraction and detection procedures described in Amalfitano and Fazi (2008). Microson XL2000 ultrasonic liquid processor with 1.6-mm-diameter microtip probe, Misonix) to detach the cells from the organic matter. The resulting slurry was left overnight at 4°C, which allowed coarse particles to settle. Thereafter, 1 ml of the resulting slurry was transferred to a 2-ml Eppendorf tube, and 1 ml of the density gradient medium Nycodenz (Nycomed) was placed underneath using a syringe needle. High-speed centrifugation was performed in a swing-out rotor for 90 min at 4°C. Nycodenz-purified subsamples (375 µl) were filtered on 0.2-µm polycarbonate membranes (47 mm diameter, Nuclepore) by gentle vacuum (< 0.2 bar) and washed with 10-20 ml of sterile ultrapure water. One section of each filter was stained for 10 min with DAPI (1 µg ml -1 final concentration) and then fixed to a glass slide to be analysed by epifluorescence microscopy. The remaining filter was stored at −20°C for further CARD-FISH analysis.
To quantify the community composition, CARD-FISH was used.
The relative abundances for the domains of Bacteria and Archaea, four subphyla of the Proteobacteria (Alpha-, Beta-, Gamma-and Delta-Proteobacteria) and the phylum Firmicutes were obtained.
In situ hybridization was carried out following the protocol of BET42a and GAM42a served as competitors for each other; for further details on probes see probeBase (Greuter et al., 2016). In addition, the abundance of photosynthetic picoplankton cells (Cyanobacteria) was estimated by their autofluorescence signal as described in Tassi et al. (2018).
The stained filter sections were observed on a Leica DM LB30 epifluorescence microscope (Leica DM LB 30, at 1,000× magnification). At least 300 cells were counted in 10 microscopic fields randomly selected across the filter sections. The relative abundance of hybridized cells was estimated as the ratio of hybridized cells to total DAPI-stained cells. (1.6 < A 260 = 280 < 1.8 and A 260 = 230 > 2) was performed by using a Nanodrop 3300 (Thermo Scientific were checked for quality by using fastqc software (version 0.11.4; Andrews, 2010) to inspect the overall quality of the sequences and look for primers, adapters and Ns content.

| DNA extraction, library preparation and sequencing
The inspection revealed no presence of Ns in the sequences, no presence of adapters and good overall quality of the sequences. By looking at the overrepresented sequences, we found out that the FWD primer was in the R2 files, while the REV primer in the R1.
The demultiplexed sequences were processed in r 3.5.1 (R Development Core Team, 2008) by first using cutadapt 1.14 (Martin, 2011) to trim the primers from the sequences.
The primers were identified by creating two objects, one for the FWD and one for the REV primer. Subsequently a function was created to detect all possible orientations for the primers.
Next, a function was applied to check the number of times the primers appeared in the forward and reverse read, while considering all possible primer orientations. Finally, the FWD, REV and their complements were trimmed off the sequences by using cutadapt. To ensure a good outcome of the trimming step, the primer count was run again on the sequences processed with cutadapt and no primers were found in all possible sequence orientations. Once the primers were trimmed, we used dada2 (version 1.10.1) (Callahan et al., 2016) to construct amplicon sequence variants (ASVs).
Taxonomic assignment to the ASVs was made by using the "as-signTaxonomy" function, which is based upon the naive Bayesian classifier method (Wang et al., 2007). The input for this command is the set of ASVs to be classified and a training set of reference sequences with known taxonomy; we used the "silva_nr_v132_train_ set.fa" (Callahan, 2018).
After taxonomic assignment, we ran the "assignSpecies" command to assign species-level taxonomy with more accuracy by using the "silva_species_assignment_v132.fa" database (Callahan, 2018).
As stated in Edgar (2017)

| Statistical analysis
Statistical analyses were performed in r, version 3.5.1 (R Development Core Team, 2008).
The environmental parameters were tested for normality using Shapiro-Wilks; only NH 4 was log-transformed (logNH 4 ) to meet normality. We tested for multicollinearity using the correlation matrix and computed the variance inflation factor (VIF) and the tolerance statistic. The analysis led us to select a few noncollinear parameters: To analyse the prokaryotic community structure three tables were created, one with ASV abundances and, from the CARD-FISH results, one with absolute abundances and one with relative abundances. To achieve equal sampling depth, we rarefied (randomly subsampling) the ASVs to the same library size number (n = 10,130, minimum number of total sequences found).
From the rarefied and standardized (by using the "decostand" function and Hellinger method) ASV abundances a Bray-Curtis dissimilarity matrix was created using the "vegdist" function. For visualization of the prokaryotic community distribution, nonmetric multidimensional scaling ordination (NMDS) was performed using the Bray-Curtis dissimilarity matrix. Starting from an initial configuration we produced 100 configurations, using the "global" model (Liu et al., 2008), and 200 as the maximum number of iterations. Unreliable distances (B-C > 0.9) were replaced by geodesic distances using a step-across method to calculate the shortest distance on any kind of "underlying nonlinear structure" (Williamson, 1978).
We extracted the two best solutions, those with the lowest stress value, and then scaled the axes of both the solutions to half change units and varimax rotation by using the "postmds" function.
To assess the fit between the two best NMDS found, we used the Procrustes comparison analysis and the "protest" function. The protest statistics (Sums of Square Difference [SSD] =1.144e-11; r = 1; permutation test [999] =0.001) confirmed the fit between the two best NMDS found. We then used the "envfit" function to fit the environmental parameters, used to produce the PCA, to see which variable was driving the community composition of microbes most.
The ordination diagram was then built with the best solution overall, with the fitted values for the water physicochemical parameters. To test for differences in the microbial community structure, we performed an analysis of similarities (ANOSIM) using the catchments as the factorial variable. This type of analysis provides statistical information on the difference between microbial communities according to the grouping variable.
We also performed a Beta dispersion analysis to test if the prokaryotic abundances from CARD-FISH and 16S sequencing were homogeneously dispersed among groups of two different factorial variables: rivers and regions. By using the "adonis" function from the vegan package (Oksanen et al., 2013), we also tested the species compositional difference between the factorial variables.
Mantel tests were run between Bray-Curtis similarity matrices for ASV abundances, CARD-FISH relative and absolute abundances, and the Euclidean distance matrix for the standardized environmental variables to detect similar patterns and thus the driving variables for the bacterial community composition.
From the ASV abundance tables at class and genus levels, we calculated the prevalence (occurrence for sampling site) and coefficient of variation (standard deviation of taxon abundance divided by the mean) to detect taxa which might be suitable as biological indicators.
These two parameters were plotted against each other in a scatter plot to find taxa with the highest prevalence and highest variation.
Taxa with high prevalence and variance could be used as indicators of environmental gradients. We performed a redundancy analysis (RDA) with the ASV abundance matrix at class and genus levels and the environmental variables used for the PCA to detect any relationship among the taxa with highest prevalence and variance and the environmental variables. The class-level ASV abundance matrix was transformed using the Hellinger method to reduce the effect of large abundances.

| Water physicochemical characteristics
The water chemistry parameters for the two regions are shown in Table S3. The two regions showed different patterns (Figure 2), but the main difference was pH, which in the Oslo region was on average 7.53 ± 0.42, while the mean value for Agder was 5.70 ± 0.40.
Mean conductivity, measured in the rivers from the Oslo region, was 86.86 ± 37.77 µS cm -1 , considerably higher than the mean values recorded in Agder (19.20 ± 9.57 µS cm -1 ). Another marked difference between the two regions analysed was the mean value for sulphate, showing higher concentration in the rivers from Oslo Measurements for TN and TP showed a similar pattern, with higher concentrations in the catchments from the Oslo region (on average 574.3 ± 20.7 µg L -1 TN, 9 ± 2.7 µg L -1 TP). The average concentration for TN for Agder was 345.6 ± 181.6 µg L -1 . The average TP concentration was 5.4 ± 3.1 µg l -1 for Oslo and 3.6 µg L -1 for Agder. As confirmed by the t test result (p < .001), pH was the variable demarcating the two regions. Results for t tests were significant also for TN (p < .05), TP (p < .05) and TOC (p < .05) but not for logNH 4 (p = .3).
Further detail are given in Table S4.
Detailed information on the abundances of the specific bacterial groups analysed is presented in Table S5. By using the data from only those rivers impacted by hydropower and dams, we Bacteroidetes (5.6%).
The NMDS plot shows a clear clustering of the sampling sites according to geographical distribution, dividing the samples from the region around Oslo from those from Agder ( Figure 5). This is consistent with the cluster analysis and the PCA conducted on the environmental variables. The Mantel test (Table 1)  To visualize the underlying trends in our ASV abundances, we performed an "envfit" analysis to fit the environmental features, highlighted in the PCA for water chemistry, to the NMDS created from the sequencing data. The correlation results confirmed that community composition and pH were closely associated (r =.9, p <.001) (Table S6).

| Comparison among methods and bacterial indicators
The results of the ANOSIM performed on the Bray-Curts matrix obtained from ASVs showed a significant association with specific catchment (r =.8, p <.001). This was confirmed by Beta Dispersion Analysis where we detected variation in the species composition at both, regional (PERMANOVA, r 2 =.274, p <.001) and catchment scale (PERMANOVA, r 2 =.461; p <.001). In comparison, the CARD-FISH results showed very high variability, and thus no significant association with specific regions or catchments could be observed.
Analysis of prevalence and the coefficient of variation at the class level showed that some taxa (at class and genus level) were distributed across all sampling sites (Figure 6a,b; Figure S1a,b), showing large variability in their abundances. The highly variable taxa with wide prevalence across the catchments and two geographical regions showed distinct patterns as revealed by RDA (Figure 6c). For example, Bacilli was positively associated with TP and Bacteroidia with ammonium, as well as Alphaproteobacteria with low pH and Gammaproteobacteria with high values for TN and TOC. Similarly, using abundance data at the the genus level revealed associations of Janthiniobacterium with TP and Sphingomonas with low values of TN and TP ( Figure S1c). As such these taxa may provide biological indicators for the status of the river system with regard to nutrients and acidification.

| DISCUSS ION
This study casts new light into the prokaryotic community structure of epilithic biofilms dwelling in rivers affected by natural and anthropogenic impacts. The combination of two techniques, Quantitative results such as those obtained by CARD-FISH allow us to detect variations in actual cells numbers and activity of specific taxa, which would be otherwise lost by analysing only the number of gene copies provided by sequencing (Fazi et al., 2020).
The CARD-FISH results obtained in our study revealed a great variation at small spatial scales, such as in biofilms belonging to the same river, and thus seems to detect in-system variability in the microbial community composition to a greater extent than metabarcoding. However, this extreme microscale variability might mask the overall effects of the main drivers for the whole microbial community composition. With regard to analysing communities occurring across large spatial scales, the sequencing methods used in our study have proven their validity. Metabarcoding provides a huge amount of data with high taxonomic resolution, which can be related to the physicochemical parameters of the environment (Ligi et al., 2014). It enables the exploration of large-scale patterns in relation to environmental conditions and a finer taxonomic resolution than hybridization methods (Bouvier & del Giorgio, 2003;Corte et al., 2013).
Both techniques revealed a dominance of Proteobacteria across all the samples corroborating most previous studies on freshwater epilithic bacterial communities (Battin et al., 2016;Besemer et al., 2012;Wilhelm et al., 2013). Overall, the total prokaryotic cell abundances obtained by DAPI staining were an order of magnitude lower than those found by Fazi et al. (2005), but comparable to those found by Zoppini et al. (2010) in similar freshwater systems. Beta-and Alphaproteobacteria were the most abundant classes according to the CARD-FISH results, in line with the results of studies on microbial communities in urban streams (Araya et al., 2003) and freshwater mesocosms (Lupini et al., 2011).
Gamma-and Deltaproteobacteria were less abundant, similar to the findings of Webster et al. (2004), where biofilms at different stages of development were analysed by FISH (fluorescence in situ hybridization). From the absolute abundances of bacterial groups obtained by CARD-FISH, we were also able to detect peaks of bacterial cell numbers, which may be related to pollution sources and impacts that would not have been identified by using sequencing alone (Freixa et al., 2016;Bakenhus et al., 2019). Nevertheless, the variability between sampling sites, within the same river, was too large to identify any associated variables providing potential explanations. We speculate that cell numbers might be affected by microscale ecological features that we did not measure, such as interactions with other biological or physicochemical components varying at the microscale. Variability might also be caused by random events such technical or sampling biases.

F I G U R E 5
To analyse the distribution of the biofilm community structure and its relationship with the environmental parameters we plotted the envfit analysis produced by using the best GNMDS out of 100 iterations, performed on the Bray-Curtis matrix of ASV abundances and the data frame for the variables used in the PCA.  (Lee at al., 2008), and the latter is also common in soils and freshwaters (Williams et al., 1996). In terms of read numbers, the second most important bacterial class was the Alphaproteobacteria with Sphingomonas being the most common genus. Sphingomonas has been previously found to be an important player in biofilm structural composition because of the high production of expolysaccharides, a major constituent of microbial biofilms (Johnsen et al., 2000).
Betaproteobacteria were the third proteobacterial group to be highly represented in the sequencing data, with Massilia and Janthiniobacterium being the most common genera for this class.
Betaproteobacteria are the group most associated with freshwater ecosystems, including important functional groups of bacteria such as ammonia oxidizers, which are vital in the global nitrogen cycle (Barberán & Casamayor, 2010;Sekar et al., 2004;Zhang et al., 2012). Massilia and Janthiniobacterium, both belonging to the order Burkholderiales, are typical of freshwater environments (Gołębiewski et al., 2017). Massilia is a ubiquitous genus, often present in soils and in biofilms and exhibiting unique properties including expolysaccharide production, incredible adhesive force and hydrophobicity, making biofilms more resistant .
According to our findings, the microbial community structure is profoundly dependent on the physicochemical features of the region, confirming previous results on microbial communities from sediments in coastal areas, estuaries and rivers (Freixa et al., 2016;Aylagas et al., 2017;Borja, 2018;Fazi et al., 2020). While several environmental characteristics were associated with the epilithic community dynamics, the driving environmental parameter appears to be the acidity of water (as confirmed by the pH results of the envfit analysis, r 2 = 0.9, p <.001), which is considerably lower in the southwestern region of Norway. It is well known that pH can influence microbial communities favouring certain strains such as members of the Alphaproteobacteria (Bragina et al., 2012;Dedysh, 2009;Goffredi et al., 2011), which were dominant in the region of Agder, where rivers had on average lower pH. In addition to the more acidic environment, the nutrient load in the rivers from Agder was generally much lower compared to the rivers flowing through the Oslo area. This is due to different anthropogenic pressures in the two regions (Peder Flaten, 1991;Nordeidet et al., 2004;Reimann et al., 2009;Johannessen et al., 2015). This characteristic might also have affected the ratio between Alpha-and Gammaproteobacteria, which in the Oslo region displayed similar cell numbers, whereas in Agder the ratio was consistently different. Overall, study regions confounded the relationships between microbial community structure and environmental variables because of their distinct differences in water chemistry. So, while our study was able to show a strong response to a number of environmental variables, we are not able to disentangle this from regional effects which would need inclusion of more regions and more sampling sites.
Our study had some limitations, and true replicates are needed to get indicator taxa. However, the comparison between the two  techniques, besides showing corresponding patterns, provide different insights into the complexity of the prokaryotic community structure of riverine epilithic biofilms. Sequencing allowed us to detect the deep diversity among the microbial taxa dwelling in different river systems, with higher taxonomic resolution than with CARD-FISH. CARD-FISH provided absolute cell numbers for specific prokaryotic groups, which is the only quantitative way, based on absolute cell numbers, to assess the composition of microbial communities (Bakenhus et al., 2019;Corte et al., 2013). CARD-FISH showed high variability at the microscale, highlighting patterns between the bacterial groups analysed that were not evident from the metabarcoding results.
Overall, our results suggest that sequencing is better suited than CARD-FISH to assess overall community dynamics. On the other hand, hybridization in situ is extremely valuable in later stage studies, aiming to analyse target taxa (i.e., indicators for pollution, diseases, eutrophication, etc.). Consequently, the use of a specific technique parallels the experiences gained from other biological groups, such as macroinvertebrates, where methods, including taxonomic resolution and enumeration, differ depending on the type of bio-assessment or scientific aims (Friberg, 2014).

| Future perspectives
Here, we show how new microbial indicators can be provided by looking at the ratios between coefficients of variation and prevalence of prokaryotic taxa detected by 16S rRNA sequencing and by using absolute abundances from CARD-FISH as a conversion factor to correct for the relative read abundances ( Figure S2). By associating ASVs at specific taxonomic levels with environmental properties we might also be able to detect prokaryotic biological indicators to be used in setting environmental quality thresholds in aquatic and terrestrial environments. Our results also indicated that communities could differ substantially between geologically distinct regions, emphasizing the need to use a reference conditions approach (sensu Water Framework Directive [WFD]) in future biomonitoring with microbial indicators. While the scope of our study was too limited to establish generalized relationships between environmental variables and microbial indicators, it strongly implied that such relationships indeed exist and could be the backbone of powerful bioindicator tools for the future, filling in the black box that currently exists with regard to large parts of the microbial communities in rivers.

CO N FLI C T O F I NTE R E S T S
The authors declare that they have no known competing financial interests or personal relationships that influenced the work reported in this paper. Exception: Alexander Eiler is the owner and founder of eDNA solutions AB specialized on bioinformatics and eDNA research and innovation. He was also a part-time employee of eDNA solutions AB parallel with this research project conducted as part of his appointment as a Professor at the University of Oslo.

AUTH O R CO NTR I B UTI O N S
LP and NF conceived the study and the experimental design. LP worked on the acquisition of data. LP analysed the data. LP, AE, SF and NF discussed the results. LP, AE, SF and NF wrote the manuscript.