Multiple processes acting from local to large geographical scales shape bacterial communities associated with Phormidium (cyanobacteria) biofilms in French and New Zealand rivers

River biofilms dominated by Phormidium (cyanobacteria) are receiving increased attention worldwide because of a recent expansion in their distribution and their ability to produce neurotoxins leading to animal mortalities. Limited data are available on the composition and structure of bacterial communities (BCs) associated with Phormidium biofilms despite the important role they potentially play in biofilm functioning. By using a high-throughput sequencing approach, we compared the BCs associated with Phormidium biofilms in several sampling sites of the Tarn River (France) and in eight New Zealand rivers. The structure of the BCs from both countries displayed spatial and temporal variations but were well conserved at the order level and 28% of the OTUs containing 90% of the reads were shared by these BCs. This suggests that micro-environmental conditions occurring within thick Phormidium biofilms strongly shape the associated BCs. A strong and significant distance-decay relationship (rp = 0.7; P = 0.001) was found in BCs from New Zealand rivers but the Bray-Curtis dissimilarities between French and New Zealand BCs are in the same order of magnitude of those found between New Zealand BCs. All these findings suggest that local environmental conditions seem to have more impact on BCs than dispersal capacities of bacteria.

In order to better understand the causes of cyanobacterial blooms in freshwater ecosystems, numerous studies have investigated the potential relationships between heterotrophic bacteria and cyanobacteria. Based on field studies 1 or cyanobacterial cultures 2,3 it has been shown that bacterial communities (BCs) display dramatic changes in their structure during the development of cyanobacterial blooms. The BCs associated with the biofilm-forming benthic cyanobacteria Phormidium have also been shown to undergo successional changes during their development 4 . These data suggest that the development of cyanobacterial blooms may exerts strong direct and indirect influences on BCs through their modification of chemical parameters of the water (e.g., pH and oxygen concentrations) and on the availability of organic matter that is generated by their proliferation.
There is limited data available on the putative existence of spatial patterns in planktonic or benthic BCs. Studies that have investigated the paradigm "everything is everywhere, but, the environment selects" proposed

Results
Environmental and biological parameters. Flow velocities at which biofilms were collected in the sampling sites in the Tarn River (see Fig. 1A) ranged between 0.2 ± 0.2 to 0.6 ± 0.2 m s −1 (lowest and highest averages by site ± SD; Supplementary Table S1). New Zealand rivers (See Fig. 1B) were within the same range (Supplementary Table S2 In the Tarn River, the proportion of cyanobacteria varied mostly according to the sampling date ( Fig. 2A,B; Supplementary Table S3). The lowest cyanobacterial proportions were observed in June 2013 and July 2014 ( Fig. 2A,B) shortly after two large flow events ( Supplementary Fig. S1), and the highest were observed in August and September.
In New Zealand rivers, significant differences were observed in cyanobacterial proportions between site WMA and the rest of the sampling sites (ANOVA test, all P < 0.004; Fig. 2C and Supplementary Table S3).
Spatio-temporal variations in bacterial community structure at high taxonomic ranks. Analysis of the 16S rRNA gene sequences showed that biofilms were dominated by cyanobacteria (mean 84.3% ± 3.5 of the total reads per sampling campaign), mainly belonging to Oscillatoriales order ( Supplementary Fig. S2) and to Phormidium (mean 98.5% ± 0.9 of the cyanobacteria reads per sampling campaign). We found the same 16S rRNA dominant Phormidium OTU (with a 97% cut-off) in all rivers (France and New Zealand), which displayed a 100% sequence similarity (Basic Local Alignment Search Tool, BLAST analysis on GenBank TM , KF770970) with Phormidium uncinatum. Excluding cyanobacteria and no-hit (CNH), the most abundant bacterial phyla were Proteobacteria (mean 69.4% ± 0.9) and Bacteroidetes (mean 23.6% ± 1.2) ( Supplementary Fig. S2). Within Proteobacteria, Alphaproteobacteria (mean 26.4% ± 4.4) and Betaproteobacteria (mean 21.3% ± 1.6) and to a lesser extent Gammaproteobacteria (mean 9.0% ± 0.8), were most prominent. Unclassified Proteobacteria also accounted for a high proportion (mean 12.0% ± 4.0) of the reads in this phylum. Within these three classes, most of the sequences were classified in the orders of Rhodobacterales/Sphingomonadales, Burkholderiales and Xanthomonadales, respectively. Bacteroidetes were mainly represented by Flavobacteriales (Tarn River) and  Table S4). At these high taxonomic ranks, the bacterial structure in biofilms displayed a strong similarity during the two sampling years in Tarn River (Supplementary Fig. S2; Bray-Curtis dissimilarity = 0.18 in Supplementary Table S5) and when comparing Tarn River and New Zealand rivers BC (Bray-Curtis dissimilarity = 0.20 in Supplementary Table S5).
When comparing the occurrences and relative abundances of dominant bacterial orders (≥1%; excluding CNH) in Tarn and New Zealand rivers, it appeared that the global structure of all the BC was maintained among sampling sites (see Fig. 3). However, when performing Principal Component Analyses (PCA) on three different data sets (Tarn 2013 & 2014 Site T1, Tarn 2014 and New Zealand rivers), spatio-temporal variations were evident in the structure of the BCs. The PCA performed on the biofilms collected in T1 sampling station from the Tarn River (Fig. 4A), highlighted on the first axis (47.2% of the total inertia) of the projection, a highly significant distinction (PERMANOVA test, P = 0.001; Table 1) between the samples collected in June and those collected from July to September. These seasonal differences in the structure of the BCs were mainly due to the relative importance of Sphingomonadales, Rhizobiales and Rhodobacteriales orders (ANOVA test, all P < 0.001; Supplementary Table S6). On the second axis of this PCA (15.2% of the total inertia; Fig. 4A), a significant difference (PERMANOVA test, P = 0.002; Table 1) was found between samples collected in 2013 and 2014. In concordance with sampling month differentiation observed in site T1, when considering all sampling sites from Tarn 2014 (Fig. 4B), the same temporal pattern was obtained (first axis explaining 36.3% of the total inertia; PERMANOVA test, P = 0.001; Table 1). The PCA performed on the New Zealand rivers highlighted significant differences (Nested PERMANOVA test, P = 0.02; Table 1) between BCs from the North and the South Islands of New Zealand. These differences were mainly due to the relative abundances of bacteria belonging to Burkholderiales and Sphingomonadales orders (Fig. 4C). Finally, significant differences (Nested PERMANOVA test, P = 0.003; Table 1) were found also at the order level, between BCs from New Zealand rivers and the Tarn River (see PCA Supplementary Fig. S3A).    Supplementary Fig. S5B). Shared OTUs (=core OTUs) contained a very high number of reads compared to those only present in biofilms from Tarn or New Zealand rivers ( Supplementary Fig. S5) and there was a clear relationship between the mean number of reads per OTU and the number of samples in which each OTU was found (see Supplementary Fig. S6). At the genus level, the OTUs belonging to the core species (≥1% of total reads) Missing data points in (A) sampling months: June, July in T2 and T4 are due to (i) high flow velocities caused by high rain events, which did not allow us to make the sampling without risk for the operators or/and (ii) the lack of Phormidium biofilm at the site during the sampling campaign. In (C) sites W and WMP are not presented due to missing data (refer to Fig. 1 for site names and locations).  Supplementary Table S7).

Spatio-temporal variations in bacterial community
In agreement with the analysis on BCs at the order level, PCA performed at the OTU level reveal the same spatio-temporal variations in BCs structure ( Fig. 4D-F). Significant temporal differences (PERMANOVA test, all P = 0.001; Table 1) in BCs were found when comparing samples collected in the Tarn River at site T1 between June and the other months ( Fig. 4D; first axis explaining 21.1% of the total inertia) and between the two sampling years (second axis explaining 14.1% of the total inertia). The PCA performed on the BC sampled in 2014 in the Tarn River ( Fig. 4E) confirmed the existence of temporal variations with a distinction among samples from June, July-August and September and showed also spatial variations among the different sites (PERMANOVA test, all P = 0.001; Table 1). The structures of BCs associated with Phormidium biofilms from New Zealand rivers located in the North Island display significant differences with those of rivers located in the South Island ( Fig. 4F; Nested PERMANOVA test, P = 0.04; Table 1). In addition, BCs from South Island biofilms were tightly grouped ( Fig. 4F; Bray-Curtis dissimilarity = 0.45 in Supplementary Table S5), whereas those from the North Island were much more dispersed (Bray-Curtis dissimilarity = 0.61 in Supplementary Table S5). Moreover, biofilms from the North Island followed a gradient of ordination from north to south on the second axis of the analysis (Fig. 4F), except for site K. Finally, significant differences (Nested PERMANOVA test, P = 0.001; Table 1) were found at the OTU level, between BC from New Zealand rivers and those from the Tarn River (see PCA Supplementary Fig. S3B).
There was a weak significant relationship (Mantel test, r p = 0.32; P = 0.04; Fig. 5) between Bray-Curtis dissimilarities and geographical distances among the BCs collected in five sampling sites of the Tarn River in August 2013 (Fig. 5). Conversely, this relationship was highly significant among BCs sampled in New Zealand rivers (Mantel test, r p = 0.73; P = 0.001; Fig. 5). Finally, while the geographical distance between France and New Zealand is much greater that the geographical distances between New Zealand rivers, the Bray-Curtis dissimilarity values estimated between BCs from Tarn and New Zealand rivers are in the same order of magnitude than those found between the most distant rivers in New Zealand ( Fig. 5; Average Bray-Curtis dissimilarity values presented in Supplementary Table S5). A similar curve was found after transforming abundance data to presence/absence data using Jaccard index of dissimilarity ( Supplementary Fig. S7).

Discussion
Numerous studies in the twenty past years have highlighted the major role microbial biofilm communities play in stream ecosystem functioning (see the review of Battin et al.) 17 . In particular, stream biofilms contain an assemblage of phototrophic and heterotrophic microbial communities, which contribute to primary production and numerous biogeochemical cycles. The relationships between phototrophic and heterotrophic taxa can be positive and negative, and these interactions may also be important in shaping biofilm communities 17,18 . In addition, it is well known that environmental factors can have direct and indirect influences on the microbial communities of stream biofilms [19][20][21] .
In this study, we investigated stream biofilms whose phototrophic component was dominated by cyanobacteria belonging to the Phormidium genus for three reasons. These cyanobacteria are potentially toxic and their proliferations seems to be increasing worldwide 22,23 . Secondly, biofilms dominated by Phormidium are found in similar environmental conditions in streams (i.e., depth, flow velocity, substrate) [24][25][26] . This reduces the number of potential environmental factors that might impact BCs associated with Phormidium and enabled a comparative approach to be performed on these communities between France and New Zealand rivers. Finally, biofilms dominated by Phormidium are characterized by a large thickness (up to several millimetres), which results in specific micro-environmental conditions and cyanobacteria probably constitutes the main source of organic compounds for the BC. Consequently, we anticipated that: (i) BCs associated with Phormidium would display a conserved structure at various spatiotemporal scales, and (ii) that the BCs differ from those living in biofilms dominated by other phototrophic microorganisms, for example diatoms, in the same kind of environments.
In agreement with our hypothesis that the specific micro-environmental conditions provided by thick Phormidium biofilms shape the BCs structure, we showed that all the BCs associated with Phormidium biofilms in the Tarn River and New Zealand rivers display similar patterns in their structure. At higher taxonomic ranks (phylum, class or order), these BCs were dominated by species belonging to Proteobacteria (Alpha-and Betaproteobacteria) and Bacteroidetes phyla, the three most abundant genera being Rhodobacter (Alphaproteobacteria), Flavobacterium (Flavobacteria) and Pedobacter (Sphingobacteria). The dominance of Proteobacteria and Bacteroidetes is a common feature of BCs in biofilms of streams 12,21,27 and lakes 28 . However, when analysing the few studies available using NGS approaches on BCs from river biofilms, it appears that these BCs are dominated by various genera, including Acidovorax (Betaproteobacteria), Flavobacterium and Polynucleobacter (Betaproteobacteria) in the study of Besemer et al. 12 ; Arenimonas (Gammaproteobacteria), Rhodobacter and Betaproteobacteria-undefined genus in the study of Bricheux et al. 29 and Rhodobacter, Exiguobacterium (Firmicutes) and Zymomonas (Alpha) in the study of Peipoch et al. 16 .
Furthermore, our data indicates that habitat type (benthic versus pelagic life) is a key factor acting on BCs associated with cyanobacteria. The importance of habitat on BCs associated with Phormidium biofilms in rivers is further supported by the fact that these communities display marked differences with those associated with planktonic bloom-forming cyanobacteria in lakes, which are commonly dominated by Betaproteobacteria 1,30 (i.e., Comamonadaceae members) 30 or Bacteroidetes 31,32 (i.e., Flavobacterium, Sediminibacterium 31 ), and in a lesser extent by Alpha-(i.e., Pelagibacter) 31 and Gammaproteobacteria 1,31 (i.e., Pseudomonas) 31 .
When comparing BC in biofilms dominated by diatoms and by cyanobacteria in Tarn River, marked differences were found between them at all taxonomic levels. The BCs associated with diatoms in biofilms were largely dominated by Alphaproteobacteria as previously observed in marine environments 33 , in epilithic lake biofilms 34 , and in river biofilms 18 . These findings suggest that BC are at least to some extent dependent on the key phototroph microorganisms in the biofilms, and the resulting micro-environmental conditions. Among the processes that could be involved in the selection of bacterial species inside the biofilms, Becker et al. 35 showed that the dissolved organic matter (DOM) produced by diatoms and cyanobacteria is very different, which potentially has a substantial influence on heterotrophic communities utilising it for their growth. Wagner et al. 36 also found that the quality of the DOM impacted the functional and structural diversity of BCs in hyporheic biofilms. Phormidium biofilms can be several millimetres thick and form very cohesive biofilms 23,37 . A recent study has shown that this creates conditions within the biofilm (e.g., pH, dissolved oxygen, nutrient and metal concentrations) that are very different from those of the overlying water 37 . It is likely that these micro-environmental conditions have a substantial impact on the BC structure, which could partly explain the dominance of particular OTUs into the BCs, regardless of their geographical origin (France or New Zealand). As Phormidium biofilms grow, these within-biofilm conditions probably become more intense and this may explain why marked shifts in BCs were observed during different successional phases 4 .
We identified a positive relationship between the abundance of reads per OTU and the number of samples in which they were retrieved. This indicates that the most abundant OTUs are widely distributed and constitute the core species of the BC associated with Phormidium biofilms in rivers. A similar positive relationship has been previously described by Humbert et al. 38 for pelagic BCs from lakes and by Dougal et al. 39 for BCs in the large intestine of horses. As members of the core species, it is probable that these OTUs play a major role in the basic functioning of biofilm BCs while satellite species distributed in a restricted number of samples, are probably involved in the adaptation to local environmental conditions. Many of the dominant OTUs in the biofilms examined in this study are well known for their ability to degrade complex compounds such as recalcitrant humic substances or organic contaminants (e.g., Sphingomonas 40 and Hydrogenophaga 41 ) or large organic polymeric proteins and polysaccharides (e.g., Flavobacterium 42 and Runella 43 ). Others have been shown to be involved in the nitrogen cycle, including nitrogen fixation and denitrification (e.g., OTUs from Rhodobacter 44 , Sphingomonas 40 , Azonexus 45 , and Hydrogenophaga 41 ). During the early stages of development of the Phormidium biofilms (June and July in Tarn River), numerous reads belonging to Rhodobacter genus were identified. This genus is metabolically diverse, with members that perform anaerobic phototrophy and photoautotrophy to aerobic chemoheterotrophy 44  Rhodocyclales, Rhodospirillales and Rhodobacterales raises questions regarding their potential contribution to primary production in river biofilms 18 .
Distance-decay similarities allow variations in beta-diversity across different spatial scales (e.g. Soininen et al. 46 ; Tan et al. 47 ) to be explored. Beta-diversity depends on the sampling scale and should increase with increasing spatial extent of the study area 47 . In our study, three different scales were taken into account, from local (intra-river), to regional (inter-river) and global (inter-country) levels. Different processes drive beta-diversity at these different spatial scales. For example, topology, altitude, discontinuous habitat, latitudinal gradients in productivity, and climate act as regional scales, while habitat structure and composition act as local filters 48 .
As expected, the distance-decay relationship estimated between BCs associated with Phormidium biofilms collected in the five sampling stations of the Tarn River (intra-river comparison) was much lower than that found among BCs collected in various rivers in New Zealand (inter-river comparison). However, it was unexpected to find that the genetic dissimilarities between BCs from France and New Zealand were of the same order of magnitude as those found between the most distant rivers in New Zealand. We had anticipated that the larger distance between France and New Zealand would limit the dispersion of bacteria, and that other processes such as climate factors and continental isolation which act as large-scale filters 48 , would result in much greater differences among these communities. This result is largely due to the presence of several core OTUs, which are highly abundant in Phormidium biofilms from both France and New Zealand.
Collectively these results indicate that in addition to micro-environmental conditions that shape the BCs associated with Phormidium biofilms, the main processes acting on these BCs seems to be mainly related to local environmental conditions rather than to geographic distance per se. This is supported by the study of Lear et al. 13 performed across New Zealand where they found that stream biofilm BCs were more influenced by catchment land use attributes than by dispersal limitation and also by Peipoch et et al. 16 in USA.
Our data are also interesting to consider in the context of the meta-analysis on stream microbial communities provided by Zeglin 20 and on the meta-community concept, which can be defined as a set of local communities linked by dispersal of multiple interacting species 49 . According to Zeglin 20 , BCs associated with Phormidium biofilms seem to be largely structured by local heterogeneity and by temporal variations, while longitudinal (from upstream to downstream) differences were less significant. The meta-community paradigms as described by Leibold et al. 49 (i.e., patch-dynamic, species-sorting, mass-effect and neutral perspectives) allow us to consider our results in a theoretical framework. The development of thick Phormidium biofilms creates well-defined local environmental conditions (e.g., habitat, chemical gradients, organic matter availability) that may delimit ecological niches leading, by a species-sorting process, to the selection of BCs displaying a similar structure. Inside these biofilms, the local dynamics of species due to stochastic and/or deterministic processes (including seasonal succession) seems to have a marked effect on the BCs associated with Phormidium biofilms as emphasized by the importance of intra-site variations.
In conclusion, our findings obtained on the BCs associated with Phormidium biofilms collected in France and New Zealand have shown that various processes functioning at different temporal and spatial scales, shape these BCs. Among these processes, we found a distance-decay relationship between BCs collected in rivers from NewZealand. Future research should focus on a large scale sampling including other countries. Moreover, this sampling should also take into account the influence of temporal variations in the structure of BCs associated with Phormidium biofilms by targeting them at different phases of their development. Finally, a significant challenge for the future will be also to establish to what extent the recent increase in benthic Phormidium proliferations in river worldwide is due to perturbations of heterotrophic BCs in biofilms, or whether it is a result of changes in environmental factors and processes promoting the development of Phormidium, with subsequent consequences on the BCs.

Methods
Sites description. Samples were collected from nine wadeable rivers in France and New Zealand. Previous studies had identified the presence of potentially toxic benthic cyanobacteria in these rivers [50][51][52] . In France, samples were collected at five sites (T1 to T5; see Fig. 1A) in the Tarn River between June and September 2013 and 2014, which is the most favourable period for the growth of Phormidium biofilms. Sites T3 and T5 were only sampled in August 2013 as part of an expanded monitoring campaign. In New Zealand, ten sites located in eight rivers were sampled in February 2013 (see Fig. 1B). The bed substrates at all sites in France and New Zealand were dominated by cobble and boulder (6 to 26 cm length), except at the Kuratau River in New Zealand (site K; Fig. 1B), where sand was the most abundant substrate. Nitrate (NO 3 -N) and total phosphorus (TP) concentrations were higher at the Tarn River (mean 1.4 ± 0.6 mg L −1 and 0.02 ± 0.01 mg L −1 respectively; data extracted from the French database SIE Adour-Garonne 53 ); compared to values from New Zealand rivers (mean 0.4 ± 0.4 mg L −1 and 0.01 ± 0.01 mg L −1 respectively; mean values collected at each sampling site according to previously described 4 ). At each sampling point, water temperature, pH, flow velocity and depth were measured. In The Tarn River, variations in flow rates during the sampling period at two survey stations (Bedoues located 28 km upstream, and Mostuéjouls located 38 km downstream from Sainte-Enimie) were obtained from hydro France database 54 . Sites TP and TW were both located in the Tukituki River (North Island) and sites WN and WMP in the Wakapuaka River (South Island; Fig. 1B). Rivers are very dynamic systems and thus biofilms developing in these ecosystems are subjected to changes occurring in water flow/velocities. Missing data points in Fig. 1 are a consequence of (i) high flow velocities caused by high rain events, which did not allow us to make the sampling without risk for the operators or/and (ii) the lack of Phormidium biofilm in the site during the sampling campaign. single field of view of an underwater viewer (707 cm 2 ), the percentage of coverage of the river substrate colonized by Phormidium biofilms was visually estimated (independently by two operators). Phormidium cover (%) for a single site was consecutively calculated by averaging the estimated values from the 10 sampling points. A single cobble with a Phormidium biofilm was sampled. Biofilms were removed using sterile tweezers. Subsamples were taken for DNA and chlorophyll-a (Chl-a) extraction, and a known area (1 cm 2 ) for later microscopic identification and enumeration (preserved immediately with Lugol's iodine solution). Samples were stored chilled in the dark, and subsequently frozen (−20 °C), except for the samples for microscope analysis which were stored at 4 °C.
For site T1 in June 2013, where Phormidium development was low (percentage coverage < 5% of the river substrate) and biofilms were thinner, sampling was undertaken at three points along a transect parallel to the water's edge and positioned at two meters from the shoreline. At each point, all cobbles (5 to 20 cm length) visible in a single field of view of the underwater viewer were collected. The cobbles were scrubbed and the biomass collected in 150 mL of river water. Aliquots (5 mL) were filtered for Chl-a (GF/C Whatman) and DNA extraction (Polycarbonate 0.2 μm GTTP Millipore). Samples were stored chilled in the dark, and subsequently frozen (−20 °C) for later analysis in the laboratory. Subsamples (1 mL) were fixed with Lugol's iodine solution and stored at 4 °C for later microscopic identification and enumeration.
Enumeration and biomass estimation of the photosynthetic community. Lugol's iodine preserved samples were homogenized briefly (Ultra-Turrax T25 IKA, Germany, 3 × 2 s, 9.5 min −1 ) to break up filaments, but avoid cell damage. The samples were diluted in Milli-Q water and identification and enumeration performed using a light microscope (200× magnification, Nikon Optiphot-2, Japan), and a Malassez chamber (Marienfeld, Germany). All cells contained in 25 squares of the Malassez chamber were counted and analysis of triplicate aliquots undertaken. Cell identification and biovolumes were calculated as outlined previously 52 .
Chl-a concentrations were used as a proxy of biofilm biomass (µg of Chl-a per cm 2 ). Biofilm samples and glass fiber filters with concentrated biomass (for site T1 in June 2013) were placed in Falcon tubes (15 mL) covered with aluminium foil, frozen until analysis, and Chl-a was extracted with methanol. Chl-a concentrations were estimated according to the protocol described previously by Echenique-Subiabre et al. 52 .

Molecular analysis of biofilm communities.
Triplicates subsamples from each site were lyophilized (FreeZone2.5, Labconco, USA) and DNA extracted from ca. 50 mg using a Power Biofilm ® DNA Isolation Kit (MOBIO, USA) following the manufacturer's instructions. Nucleic acid concentrations were determined by spectrophotometry (Nanodrop 1000, Thermo Fischer Inc, USA) and samples stored at −20 °C. A region of the 16S rRNA gene (~394 base pair (bp)) including the variable region V4-V5 was amplified and sequenced using the 515 F and 909 R universal primers 55 56 . Briefly, after denoising (i.e., removing of short sequences and singletons) and chimera checking, the remaining sequences were clustered using UPARSE pipeline 57 at 97% similarity threshold. The centroid sequences of each cluster were run against a database of high quality sequences derived from the NCBI database using the USEARCH global alignment algorithm 58 . The number of read per sample ranged from 604 to 57,391. Of the 95 biofilm samples sequenced, three were removed from our analyses because an insufficient number of reads was obtained. Rarefying (randomly down-sampling using software R version 3.1.0 59 ) resulted in a total of 820,456 16S rRNA sequences (8,918 sequences per sample), which clustered in 2,198 operational taxonomic units (OTUs). In order to study the BCs associated to Phormidium biofilms, a second analysis was performed after removing cyanobacteria and no-hit (CNH) sequences. The number of reads ranged between 32 to 30,560 per sample and a total of 402 sequences per sample were selected for further analysis. Thus, only 82 samples were retained (at least two replicates per sampling site). Rarefying resulted in 32,964 16S rRNA sequences, which clustered in 1,004 OTUs.
OTUs were classified as abundant when they comprised ≥1% of the total sequences, intermediate when between <1% and >0.01%, and rare when ≤0.01% as previously defined by Mangot et al. 60 .
Statistical analysis. The R software (version 3.1.0 59 ) was used for statistical analyses of the data. Differences were considered significant when P < 0.05. The effect of sampling site and date on cyanobacteria biovolume proportions and abundance of microbial taxa were analysed using two-way ANOVA, followed by multiple comparisons using Tukey's honestly significant difference (HSD) post-hoc test.
Rarefaction curves were calculated using package 'vegan' 61 . Hellinger transformations were undertaken on Illumina sequencing data. Beta-diversity was calculated to measure the variation in BC among sites and dates using Bray-Curtis dissimilarities performed at order and OTU levels (package 'vegan'). Principal component analyses (PCA) were performed on rarefied and transformed data (Hellinger transformation) using the ' Ade4TkGUI' package 62 in order to visualise samples distribution based on their order and OTU structures. Permutational Multivariate ANOVA (PERMANOVA) 63,64 was performed on Bray-Curtis dissimilarities to test for differences in BC among sites and among sampling times (countries, sites and dates) including factorial and nested designs using package 'vegan' and 'BiodiversityR' 65 respectively. Finally, the relationship between geographical distance and Bray-Curtis dissimilarity was performed using the package 'fossil' 66 , longitude and latitude coordinates were transformed to kilometres as measured in a straight line. In order to remove temporal variability (month and year), to analyse exclusively Phormidium biofilm samples, and to capture the totality of the sampling sites, only samples from August 2013 were considered for the Tarn River. Therefore, this analysis was performed on biofilms collected at the five sampling sites for the Tarn River, and nine sampling sites from New Zealand rivers (samples Scientific REpoRtS | (2018) 8:14416 | DOI:10.1038/s41598-018-32772-w from site WMA were excluded because biofilms were dominated by diatoms). Mantel statistic based on Pearson's product-moment correlation (package 'vegan') was performed between dissimilarity matrix and permutation (999) used to evaluate statistical significance.