We sampled 218 European freshwater lakes and ponds from sites in Norway, Sweden, Germany, Poland, the Czech Republic, Slovakia, Hungary, Romania, Austria, Italy, France, Spain and Switzerland. Site selection focused on the European orogens, specifically the Alps, the Pyrenees, the Apennine, the High Tatras, the southern Scandinavian mountains and the connecting flatlands (Suppl. materials 1, 2). We selected natural lakes (and reservoirs) which were typical for the respective area. Lakes which were known to be extremely deviate such as e.g. acidic mining lakes were excluded. Plankton samples were taken at a depth between 0.2 and 0.8 m near the shoreline and immediately processed after sampling. All samples were taken in August 2012. Temperature, pH and conductivity were measured directly in the field using a portable instrument. Climatic data including photosynthetically active radiation (PAR) and UV irradiation were partly measured on site and partly taken from databases of the Deutsche Wetterdienst and associated weather services. For molecular analyses, water samples were filtered onto 0.2 µm nucleopore filters (47 mm diameter) until the filters clogged, i.e. several hundred milliliter up to 600 ml water per filter. Subsequently, filters were air dried and then immediately frozen in liquid nitrogen (Cryoshippers). The filters were stored at -80 °C in the laboratory until further processing.

Genomic DNA was extracted using the my-Budget DNA Mini Kit (Bio-Budget Technologies GmbH) following the protocol of the supplier with the following modifications: Filters were homogenized in 800 µl Lysis Buffer TLS within lysing Matrix E tubes (MP Biomedicals) using the FastPrep instrument (MP Biomedicals). Homogenization was run three times for 45 seconds each at a speed setting of 6 m/s and then incubated for 15 min at 55 °C. The next steps followed the standard protocol supplied by Bio-Budget Technologies GmbH. The quality and quantity of the DNA was checked using a Thermo Scientific NanoDrop® ND-2000 UV-Vis spectrophotometer (Thermo Fisher Scientifics). PCR amplifications targeted the SSU V9 region and ITS1. Briefly, in order to cover a broad taxonomic spectrum, two primers with different wobble positions were combined in a ratio of 10% : 90%: 5’-GCTGCGCCCTTCATCGKTG-3’ (ITS2_Dino; 10%) and 5’-GCTGCGTTCTTCATCGWTR-3’ (ITS2_broad; 90%). We used the Amplicon-Duo pipeline (Lange et al. 2015) to separate true biological variation from sequence artefacts. For each sample, two technical replicates of the extracted DNA were independently amplified using primers with different sample identifiers [see Lange et al. 2015 for details]. The concentrations of the PCR reaction were as follows: 1 μl of DNA template (depending on the concentration, different dilutions of 1:1, 1:10, and 1:100 were used) in 25 μl PCR reaction with 0.4 units of Phusion Taq (Biozym), 0.25 μM primers, 0.4 mM dNTPs and 1 x Phusion buffer (Biozym). The PCR-cycling conditions included an initial denaturation step at 98 °C for 3 min followed by 35 cycles each including a denaturation step at 98 °C for 30 s, annealing step at 52 °C for 75 s, and an elongation step at 72 °C for 60 s. At the end the PCR was completed by a final extension step of at 72 °C for 10 min. Equimolar subsamples were pooled and commercially sequenced using paired-end HiSeq 2500 sequencing, applying 2x 300 bp reads using the “rapid run” mode on the Illumina platform of a sequencing provider (Fasteris, Geneva, CH).

Adapter and quality trimming was performed by the sequencing company. Samples were demultiplexed by the sequencing company (Fasteris) using MID sequences. Base quality of the sequence reads was checked using the FastQC software (Andrews 2015). A split-sample filtering protocol for Illumina amplicon sequencing was used as described in Lange et al. (Lange et al. 2015). In brief, the raw sequences were quality filtered to remove reads with an average Phred quality score below 25 using PRINSEQ-lite [v0.20.4, (Schmieder and Edwards 2011)]. Additionally, all reads with at least one base with a Phred quality score below 15 were removed provided that the sequences contained Ns after base calling. The paired-end reads were assembled and quality filtered with the tool PANDASeq [v2.10, (Masella et al. 2012)]. Reads with uncalled bases, an assembly quality score below 0.9, a read overlap below 20, or a base with a recalculated Phread-score below 1 were discarded. All reads were dereplicated, whereupon chimeras were identified using UCHIME [usearch v7.0.1090, (Edgar et al. 2011)] with default parameters. Finally, sequences that were not found in both sample branches using AmpliconDuo (Lange et al. 2015) were discarded.

Reads remaining after the AmpliconDuo filter were clustered via the software SWARM (version 2.1.9; Mahé et al. 2014), then clustered by identical V9 sequences (150 basepairs, ident = 100 %; R-Script "V9_Clust.R" by (Jensen 2017) and subsequently taxonomically assigned by searching the NCBI database using BLASTn. The taxonomic resolution therefore roughly corresponded to a 100% SSU V9 sequence identity. All reads assigned to Metazoa and Viridiplantae as well as higher fungi were excluded from further analyses, as protists are the target organisms of this study. The remaining read counts ranged between 19,267 (Lake Estany de Trebens) and 2,330,071 (Lake Millstätter See) per sample. OTU richness and abundance were computed from not transformed (raw) data; community structure indices were also calculated based on the raw read abundance of OTUs. Rarefaction curves were plotted by the function “rarecurve” of the R-package vegan (version 2.4-1; Oksanen et al. 2017). In order to achieve a balanced dissimilarity matrix of all sites (accounting for differences in sequencing depth between the sites, etc.), we first applied Hellinger normalization. For comparison of community structures of the different sites we subsequently used the percentage difference dissimilarity matrix (Bray). We finally checked that the resulting dissimilarity matrix was euclidic (Dray et al. 2017). The resulting dissimilarity matrix was used for all further analyses. We calculated the local contribution of every site to beta-diversity (LCBD) (Dray and Dufour 2007) using again the percentage difference dissimilarity of the Hellinger-transformed data (see Suppl. material 3 for the Bray distance matrix). For further insights into the processes determining high betadiversity we also analysed the relative contribution of replacement of OTUs between lakes and abundance changes using the function beta.div.comp of the R-package adespatial (Legendre 2014b). By the parameters coef="S", quant=T this analysis was again based on Hellinger transformed (euclidic) Bray-distance data. Geographic groups of community compositions were calculated by geo-constrained hierarchical clustering analysis (Legendre 2014) resulting in 14 geographic clusters. The R-package tripack (Renka et al. 2016) was used for construction of the geographic neighbour matrix based on Delauney triangulation. The effect of geographic distance was investigated separately for lowland and high mountain lakes (see supplementary meterial 4). We calculated the distances [km] of all site pairs as well as the community distances of all site pairs (submatrices of the community dissimilarity matrix) and analyzed the effect of geographic distance by linear regression after conversion of geographical latitude and longitude data into distances [km] (Chambers 2013). The effect of mountain ranges (in particular of the Alps and the Pyrenees) as geographic barriers was also analyzed by linear regression based on splitting the sample set into three groups: an alpine group including the Alps, the Pyrenees and the High Tatras, a northern group (cis) and a southern group (trans). The mountain group was identified by altitudes plus constrained clustering (see above). By exclusion of the alpine group we could compare subsamples, i.e. all subsamples “cis-cis” (northern), all subsamples “trans-trans” (southern) and all cis-trans comparisons. Mean differences of the cis-cis, trans-trans, cis-trans groups were tested by t-tests. All figures were prepared with R and Adobe Illustrator.

Sequencing of the European freshwater samples resulted in 323,430,198 cleaned V9 forward and reverse read pairs. We applied the AmpliconDuo filter (Lange et al. 2015) for separating true biological diversity from PCR and sequencing artefacts. This novel cleaning step based on technical replication effectively discarded sequence artefacts, resulting in a total number of 100,516,809 reads. Operational taxonomic units were constructed with the clustering algorithm SWARM, resulting in 74,713 distinct unique V9 sequence SWARMs (i.e., clusters). 64,969 V9 SWARMs (95,592,035 reads) of these unique sequences were assigned to protists. The remaining 13 % of V9 SWARMs (4.9 % of reads) which were excluded from further analysis were affiliated with Metazoa (1,102 V9 SWARMs; 638,687 reads), Embryophyta (599 V9 SWARMs; 291,307 reads), the fungal lineages Basidiomycota, Ascomycota and Glomeromycota (5,110 V9 SWARMs; 3,856,310 reads) or not affiliated with the domain Eukarya (2,933 V9 SWARMs; 138,470 reads).

The applied sequencing depth was sufficient to approach saturation of eukaryotic taxon richness both at the local and European scale (Fig. 1) and resulted in 2 x 10⁴ to 2 x 10⁶ cleaned reads per sample (average 438,495 reads per site). At the European scale, the rarefaction slope dropped below 10^-4 above 63,693 OTUs corresponding to 67,000,000 reads. Towards the end of the rarefaction curve the slope was 4.44 x 10^-8, correspondingly one new OTU can be expected for 22.5 million new reads. Using the same saturation criterion, i.e. a rarefaction slope of 10^-4, local richness approached saturation between several hundred and 4,000 OTUs, corresponding to approximately 0.5 – 6 % of the European scale richness. At the level of taxonomic groups, ciliates (34.6 %) followed by dinoflagellates (12.7 %), green algae (12.0 %), cryptophytes including katablepharids (8.1 %), chrysophytes (6.6 %) and diatoms (4.2 %) contributed most to richness (Fig. 2). The individual OTUs contributing most strongly to overall eukaryotic reads were affiliated with an uncultured katablepharid, the green alga Desmodesmus insignis, and several ciliates.

Figure 1.

Planktonic protist ribosomal diversity. (A) V9 rDNA OTUs rarefaction curve. (B) Saturation slope versus number of V9 rDNA reads for all of the 218 samples analyzed herein. A slope of 10^-6 indicates that one novel sequence read can be recovered if one million new reads are sequenced. Overall, the saturation values indicate that the extensive sampling effort uncovered the majority of eukaryotic ribosomal diversity both at the European and the local level.

Figure 2.

Phylogenetic affiliation of reads and OTUs to eukaryotic taxonomic groups. The area of the portion of the pie chart represents the number of reads, whereas the angle represents the relative contribution to OTU richness. Alveolates and in particular Ciliophora and Dinoflagellata contributed nearly 50 % to global OTU richness. Cryptophyta, Viridiplantae (without Embryophyta), as well as Chrysophyceae and Diatomea accounted for the majority of the remaining OTUs. Unlabeled pieces within the eukaryotic supergroups comprise all other taxa affiliated with that group.

The OTU pool showed a good fit to the log-normal distribution (Fig. 3a). Variation of V9 OTU richness between distinct lakes was high, with a mean of 438,496 ± 423,896 V9 reads per lake. Community overlap between lakes was generally low, with only a few widespread taxa (Fig. 3b, Suppl. material 1). Taxon inventory differed with elevation, and high mountain lakes had a particularly low taxon overlap with lowland lakes. High mountain lake communities thus distinctively differed from their surroundings. Accordingly, the strongest community dissimilarities between neighbouring lakes occurred along mountain ranges (Fig. 4, Suppl. material 2). Furthermore, different high mountain areas such as the Pyrenees, the Alps, and the High Tatras decisively differed in taxon inventory from each other. High mountain areas are thus biogeographic islands to protist distribution.

Figure 3.

Abundance distribution and ubiquity of OTUs. A) Global OTU abundance distribution and fit to the Preston log-normal model. Quasi-Poisson fit to octaves (red curve) and maximized likelihood to log2 abundances (blue curve). Most OTUs were represented by 3 to 32 reads. Approximations fit to the Preston log-normal model. The maximized likelihood to log2 abundances (blue curve) fit was superior and subsequently used to calculate the Preston veil, which infers the number of OTUs that we missed during our sampling. Preston veil was 8,111, confirming that we captured most of the diversity as also indicated by the rarefaction slopes. Preston vail and rarefaction slopes confirm that holistic and general patterns of eukaryotic freshwater diversity can be extracted from our data. The unit of the y-axis is the number of OTUs B) Distribution of OTUs among European lakes. Few OTUs were ubiquitous (occurring in the majority of lakes) whereas most OTUs were specific to one or only a few lakes. Unit of y-axis as in a).

Figure 4.

Delaunay triangulation plot based on the investigated lakes. The red squares indicate triangles of high dissimilarity (i.e. larger than mean + one standard deviation) between the three corner sampling sites. Community dissimilarities were calculated based on wheighted distances using the distance function x=((1-distance)/maximal triangle leg). The areas of high dissimilarity occur nearly exclusively along mountain ranges, in particular along the Alpes, the Pyrenees, the High Tatras, the Carpathian mountains and the Sierra Nevada.

Beyond a shift in taxon inventory with elevation, richness pronouncedly decrease around 1400m (Suppl. material 4) and V9 OTU richness in high mountain lakes above 1,400m (456 ± 439) was generally lower than in low elevation lakes below 500m (843 ± 653). Most V9 OTUs were specific to a few sites, and thus lakes with a high proportion of OTUs with restricted distribution considerably contributed to beta-diversity. Contribution to beta-diversity was particularly high in high mountain regions (Fig. 5a), indicating a higher degree of habitat isolation due to either ecological factors, i.e. habitat heterogeneity, or geographic factors, i.e. isolation due to geographic barriers. The high beta-diversity data and the high values within the distance matrix indicated that different lakes harbour primarily different OTUs. The analysis of beta-diversity for all lakes indicated that 75.5 % of beta diversity was due to replacement and 24.5 % of beta diversity was due to abundance differences. When the analysis was restricted to lowland lakes of the northern area (cis) the contribution of replacement was 73.3 % and for the lowland lakes of the southern area (trans) to 79.3 %. For lakes above 1400 m beta-diversity was increased and 79.7 % of this beta diversity was due to replacement but only 20.3 % due to abundance differences.

Figure 5.

OTU distribution and protist biogeographic regions in Europe. A) Richness and contribution to beta diversity of individual lakes. Richness is illustrated by grey scale of the symbols: dark shades indicate a high richness whereas light shades indicate a low richness. The size of the symbols indicates the relative contribution to beta diversity, with large symbols indicating a high contribution. Richness of individual lakes varied between 15 and 4092. Lakes in high mountain regions tended to contribute at an above-average level to beta diversity. B) Protist biogeographic regions as calculated by geo-constrained hierarchical clustering. High mountain ranges appear as biogeographic islands; the geographic regions largely comprising high mountain protist communities are depicted in yellow. Furthermore, the lowland areas are biogeographically structured which is mainly due to an east- west-gradient. Even though mountain ranges affect protist community similarity between sites, they do not necessarily separate biogeographic regions.

When restricting the analysis to presence/absence data most of the beta-diversity was due to replacement and only a minor part due to loss of species (nestedness), i.e. 4,0 % for all lakes. However, the nestedness was somewhat higher within the northern low elevation (cis) group, i.e. 5.4%, possibly (even though hypothetical) to incomplete recolonization since the last glaciation events whereas in high elevation lakes which possibly act as refuge the nestedness is very low, i.e. 2.1 %.

On the European scale environmental factors slighty co-varied with distance for lowland lakes but did not systematically change with distance for high mountain lakes (Fig. 6A). In contrast, community similarity generally decreased with increasing distance both for highland and lowland sites (Fig. 6B). A co-inertia analysis using the R-package ADE4 of the normalized community structures (Hellinger-Bray transformed) versus the scaled environmental variables pH, temperature, and log conductivity indicated no strong correspondence between community composition and environmental factors (RV value of 0.22) thus corroborating a contribution of geographic distance and historic factors to patterns community composition.

Figure 6.

Effect of geographic distance on community dissimilarity. A) Habitat dissimilarity for high mountain lakes (above 1,400 m) and for lowland lakes (below 500 m). The latter were either separated (cis-trans pairs) or not separated by mountain ranges (cis-cis pairs, i.e. north of the Alps, Pyrenees and High Tatras and trans-trans pairs, i.e. south of the Alps, Pyrenees and High Tatras). Regression lines are shown for highland (blue), lowland cis-cis and trans-trans pairs (yellow), and lowland pairs separated by mountain ranges (orange). Habitat dissimilarity increased with distance for lowland lakes but did not change for high mountain lakes. B) Bray-Curtis community dissimilarity for the same groups as in A. Community dissimilarity in high mountain regions increased with distance even though habitat dissimilarity did not change with distance. For lowland lakes both, habitat characteristics and community composition, changed with distance. However, the effect of distance on habitat dissimilarity was independent on whether the lakes were separated by mountain ranges or not. In contrast, community dissimilarity was higher at low and intermediate distance for lakes separated by mountain ranges. Accordingly the slope of the regression line for cis-trans pairs was significantly lower as for cis- cis and trans-trans pairs (p = 0.0097). The effect of geographic barriers levels off for longer distances above approximately 1,500 km.

Beyond the general effect of geographic distance, geographic barriers such as mountain ranges contributed significantly to community dissimilarity between lowland sites on a European scale. Mean community similarity was higher for lakes within a given lowland area as compared to community similarities of lakes between areas which were separated by a mountain range (Fig. 6B). For instance, separate analyses of sites north of the Alps and south of the Alps resulted in comparatively high community similarities over a short distance, and a pronounced decrease of community similarity with geographic distance. In contrast, community similarity was generally low for habitat pairs across a mountain range already for short distances and geographic distance had a minor effect on community similarity in the latter case. Thus biogeography, in particular geographic distance and geographic barriers, shapes protist distribution pattern. Interestingly, distant lakes still share a fraction of their taxon inventory, clearly demonstrating that geographic distance and geographic structures such as mountains are not absolute barriers to protist dispersal. The effects of geographic distance and geographic barriers are, however, not strictly additive. In contrast, both factors add to community dissimilarity up to a certain level of background community similarity shaped by environmental variables.

Despite a generally low community similarity between different lakes, groups of lakes sharing a higher level of community similarity can be identified. Lakes sharing a higher degree of community similarity form clusters reflecting biogeographic provinces on a European scale (Fig. 5b). Notably mountain ranges are borders to protist biogeographic provinces as for animals and plants as well. This pattern was consistent in particular for the Alps, the Pyrenees and the High Tatras. The taxon inventory of high mountain regions strongly deviates from that of lowland regions, and lowland areas are separated from each other based on their taxon inventory into distinct provinces. These provinces do not correspond to physicochemical patterns or to the distribution of soil types or geology in Europe, suggesting that factors related to geographic distance and dispersal are key reasons for this separation. It must be noted that towards the edges of the investigation area the delimitation of such biogeographic regions become more uncertain due to the shortage of neighboring sampling sites or sites in between. For instance, whether the lakes indicated by the blue squares in eastern Spain and northern Italy represent a connected biogeographic province or rather two separate areas which share a high similarity by chance remains uncertain.

The decreasing effect of mountain ranges in separating protist communities with distance between sites indicates that geographic barriers act only as relative dispersal barriers for protists; they do not inhibit dispersal, but rather slow it down and thereby largely act on short to intermediate distances. For larger distances, the separating effect of geographic barriers adding to that of geographic distance alone largely vanishes, as demonstrated by cis-cis and trans-trans (i.e. within a distinct lowland area) versus cis-trans (i.e. lowland sites separated by a mountain range) comparisons of habitat dissimilarities (Fig. 6B).

Protist diversity is undoubtedly tremendous (Bass 2004, Countway et al. 2005, Countway et al. 2007, de Vargas et al. 2015, Moreira and López-Garcı́a 2002). Considerable research effort has been focused on marine diversity (Countway et al. 2005, Countway et al. 2007, Campo et al. 2016, Lovejoy et al. 2006, López-García et al. 2001, Massana and Pedrós-Alió 2008, Worden et al. 2006). Despite the global importance of freshwater resources and the heterogeneity of freshwater systems, their biological diversity has been comparatively neglected (Filker et al. 2016, Grossmann et al. 2016a, Kammerlander et al. 2015, Lepère et al. 2007, Simon et al. 2015, Taib et al. 2013, Triadó-Margarit and Casamayor 2012, Slapeta et al. 2005). Our study revealed a high diversity of microbial freshwater eukaryotes and a high level of regional separation. The majority of OTUs was affiliated with alveolates, indicating a high diversity within this group in particular. Within alveolates the Ciliophora showed the highest richness and comprised both, taxa with a restricted distribution and taxa with a broad distribution. In contrast, dinoflagellates and chytridiomycetes comprised a high number of OTUs specific to one or a few lakes indicating a high degree of specialization, whereas in Cryptomonadales and Katablepharidae most OTUs were generalist taxa occurring in many lakes.

Although comparatively smaller, freshwater systems are usually more heterogeneous than marine systems and offer a larger array and diversity of ecological niches (Simon et al. 2015). Freshwater communities are therefore expected to be much more differentiated between sites than marine communities (Barberán et al. 2011, Logares et al. 2009). In fact, we found a low taxon overlap between lakes and a correspondingly high contribution to beta-diversity of individual habitats. In accordance with our findings, a large genetic divergence between distinct lake microbial communities has been reported in recent deep-sequencing campaigns, in particular for the rare biosphere (Filker et al. 2016, Grossmann et al. 2016b, Triadó-Margarit and Casamayor 2012). In particular, distribution patterns of rare taxa seem to vary widely and to deviate from the rather uniform patterns in abundant protists (Grossmann et al. 2016a, Schiaffino et al. 2016). This confirms recent claims that distribution patterns are almost certainly differential, with some taxa being ubiquitous whereas others are rather endemic, as proposed for instance by the moderate endemicity model.

Our survey aimed to cover all eukaryotic lineages, but universal primers may miss some relevant taxonomic groups (Massana et al. 2015). Distribution patterns are further obscured by partly inappropriate species delimitation in protists (Mann and Vanormelingen 2013, Scoble and Cavalier-Smith 2014). Appropriate species delimitation (and thus species distributions) are even more nebulous for the vast number of (so far) uncultivatable protist lineages (Massana et al. 2013). Assertions on the absolute species numbers and species inventory are therefore imprecise and subject to numerous assumptions. It is therefore important to keep in mind that some taxonomic groups are underrepresented and may exhibit patterns deviating from the general trend. Keeping this in mind, the community composition corresponded to those known from other studies (e.g. Lange et al. 2015, Kammerlander et al. 2015, Filker et al. 2016).

Our analysis based on 218 European lakes yielded in 74,713 distinct unique V9 sequence SWARMs. Rarefaction analysis as well as Preston vein indicated that we captured most of the diversity. However, molecular diversity inventories based on high throughput sequencing techniques must be interpreted with caution as high throughput sequencing platforms are error prone and in particular for sequences exclusively present in one or very few samples PCR artefacts and sequencing errors must be taken into account (Degnan and Ochman 2011, Lange et al. 2015). The AmpiconDuo filter largely removed such sequences by splitting the samples followed by two independent PCRs and sequencing steps for each sample. Only errors that occurred independently in both branches of the split sampling protocol may pass the filter. Even though it is highly unlikely that high numbers of artificial sequences pass the filter it is known that some artefacts such as formation of some chimeric sequences may, in part, occur systematically (Fonseca et al. 2012) and thus escape the filter. In fact, particularly some chimeric sequences pass the AmpliconDuo filter and not all of these sequences are recognized by chimera filters such as UCHIME (Lange et al. 2015). Therefore it is likely that some artificial sequences passed the filter steps, even though the share of sequence artefacts is probably lower as compared to past standard filtering algorithms. However, erroneous sequences should mostly be restricted to a single sample and thus not or hardly affect our analyses of community pattern across samples. The high number of OTUs restricted to a single lake may, nevertheless, be overestimated. It is likely, however, that due to the strict filtering steps the majority of these OTUs represent true biological diversity.

For some samples the number of OTUs was considerably lower as expected. This may be due to several factors: First, our filter strategy, including the AmpliconDuo filter and SWARM reduced the number of OTUs considerably. In particular, AmpliconDuo filtered out many presumably artificial sequences most of which had low to medium read numbers. Therefore the total number of OTUs should be expected to be somewhat smaller as compared to studies which do not apply this or a comparable filter (Lange et al. 2015). Second, the comparatively long amplicon (V9 + ITS1) and the long primer constructs also may miss some sequences, as it is generally known for long primer constructs and long amplicons (Huber et al. 2009). Third, some samples were strongly dominated by one or very few taxa due to algal blooms or strong dominance of a distinct bacterivore. For instance, Lake Pfordter See which was strongly dominated by a katablepharid yielded in as few as 15 OTUs. In such communities a lower taxon coverage was to be expected which may even have been self-inforcing as the dominance of one distinct sequence may cause a negative selection of rare sequences during PCR thereby further decreasing detection limits. Even though the total number of OTUs was certainly lower than the actual number of species present, the data structure as well as observations indicate that biological reasons, i.e. community structure, is presumably largely responsible for low OTU numbers in these lakes. In sum, even though in some lakes OTU number was surprisingly low, data structure and observations during sampling indicate that a strong dominance by few OTUs in these lakes is presumably largely real. Additional taxa, even though presumably extremely rare, are certainly to be expected. However, for the comparative analyses the applied Hellinger transformation balanced potentially deviating distributions and the analysis of distribution patterns is much less affected when based on defined units of biological diversity such as the herein-used V9 OTUs.

Ecological factors, in particular abiotic factors such as nutrient concentrations, temperature and light and UV intensity are well known to affect protist community composition (e.g Sommaruga 2001, Sonntag et al. 2010, Tolotti et al. 2003, Tolotti et al. 2006, Triadó-Margarit and Casamayor 2012). However, recent data provide increasing evidence that beyond the ecological factors also geographic distance and historical factors affect protist distribution at least on large scales (Filker et al. 2016, Schiaffino et al. 2016, Fernández et al. 2017). Europe is a well-investigated (if not the best-investigated) area for biogeographic patterns (Hewitt 2000, Schmitt 2009). The European mountain ranges – the Alps in the center of Europe and further mountain ranges towards the west (Pyrenees, Massif Central), the south (Apennine), the east (High Tatras, Carpathian Mountains) and to the north (Scandinavian mountain ranges) – constitute a unique area for studying the effect of distribution barriers.

Based on the largest freshwater data set known to us, we demonstrate a significant contribution of geographical distance and historical factors to protist distribution patterns. We found spatial patterns of protist diversity with microbial community composition varying both with altitude in mountain ranges and with geographical distance. Seasonal taxon fluctuations influence protist community composition and may interfere with spatial analyses when samples originate from different seasons (Bock et al. 2014, Massana et al. 2015, Nolte et al. 2010). Considering that all samples in our campaign were taken in August 2012 and that seasonal effects and plankton succession dynamics is presumably lower in mid-summer lake plankton communities than in spring succession communities (Naselli-Flores et al. 2003), the observed divergence presumably reflects largely spatial differentiation of lake plankton communities rather than seasonal variation (Grossmann et al. 2016b). However, due to the large area included in the analysis a certain bias due to differential stages of the seasonal succession cannot be excluded.

Recently, de Vargas et al. 2015 have shown that the correlation between community similarity and geographic distance in marine ecosystems levels off for distances greater than 6,000 km. Our data confirm a similar relationship for limnic lowland systems which act, however, at shorter distances and is modified by local and regional geography. In particular, Bray-Curtis dissimilarity (supplementary material 3) increased with distance, but this correlation already levelled off for large distances and virtually vanished for sampling sites separated by geographic barriers. This indicates that protist freshwater communities in Europe are presumably largely driven not by geographic but by ecological factors.

For instance, we found a general trend of decreasing richness for high elevation lakes, i.e. high mountain lakes at elevations above approximately 1400 m showed a lower mean richness as compared to lakes at low elevations. This may be due to either the generally small size of mountain lakes as would be expected from species-area-relationships or indicate differential patterns of richness depending on elevation. As small lakes at low elevation can show a high richness we suspect that both, habitat size and elevation, or factors co-varying with elevation, may be repsonsible for the low richness.

Even though ecological factors are certainly of primary importance explaining for protist community composition, the contribution of geographic distance and historical factors (or geographic barriers) seems to have a stronger impact on protist freshwater communities than on marine communities (cf. de Vargas et al. 2015). Further, for high mountain ranges the effect of geographic distance is more pronounced as these regions are more isolated and may act as biogeographic islands at least to species with restricted dispersal capabilities. In this context it must be noted, that a considerable portion of protist distribution pattern is not explained by geograhic separation indicating that it is either restricted by other factors, i.e. presumably ecology, or not restricted by the geographic barriers, i.e. suggesting that the significance of geographic barriers may strongly differ between taxa. However, our data provide evidence that historical factors contribute to shaping protist distribution pattern (Fernández et al. 2016) beyond the presumably predominant effects of ecological factors.

According to the generally high beta diversity, our study revealed a high proportion of taxa exclusively found in only one or a few lakes. The very low contribtion of nestedness of the high elevation lakes and the somewhat higher contribution of nestedness to beta diversity in the northern lowland lakes may be interpreted as historical signals of the last glaciation and the post-glacial colonisation. However, as these differences are weak the interpretation as a historical signal remains hypothetical. How far postglacial resettlement, geographical barriers, other historical reasons, ecological parameters like nutrient supply or just biotic influences by other organismic groups dominate the exceptionality of the protist communities within lakes, remains to be shown. Further, it is likely that ecological and historical factors interact in shaping protist distribution pattern (Fernández et al. 2017). Irrespective of protist community differentiation due to ecological factors and geographic distance, we found a weak but nevertheless significant effect of mountain ranges as geographic barriers contributing to the separation of low elevation lakes (fig. 6). Future large scale analyses including cross-continental comparisons of such large datasets, in particular of lakes of a similar trophic status, will be a good test to confirm large scale patterns of geographic differentiation in protists and to differentiate regional from global patterns.