Main

Nutrients, materials and organisms transferred from terrestrial environments to freshwater networks play a crucial role in global biogeochemical cycles1,2,3. The shoreline of small lakes and ponds represents the largest surface area among all inland waters and thus serves as an important terrestrial–aquatic interface3,4,5. Uptake, transfer and release of nutrients occurring in the food web of small waterbodies thus have considerable consequences for the entire downstream aquatic network6,7,8,9. At the same time, small water volume and the tight linkage to the surrounding land render these waterbodies highly vulnerable to fluctuating hydrological conditions. Extreme rain events (ERE) shorten water retention times to values characteristic for streams, cause washout of planktonic organisms and bring extreme loads of nutrients and organic matter from flooded watersheds9,10,11,12,13. Such strong disturbances on planktonic communities can result in delaying, arresting or diverting seasonal succession from its usual pattern14,15.

Intensification in strength and frequency of EREs due to an ongoing and future climate change is predicted for many areas worldwide, including Central Europe16,17,18. Although terrestrial–aquatic boundaries have a crucial influence on the entire landscape as highly responsive hotspots of organic matter transformation, energy and nutrient flow19, flood-induced effects on nutrient budgets, composition and function of plankton communities as well as on food web structure remain largely unknown. While communities of small waterbodies cannot possess resistance to EREs due to direct washout effects, their resilience capacity20,21 is of high interest, raising an intriguing question to which extent these ecosystems are able to absorb EREs or if such disturbances force them beyond their resilience limits towards new alternative steady states22,23. To increase the predictability of consequences caused by EREs, extensive high temporal resolution studies along with an application of integrative and dynamic concepts extendable from population to meta-ecosystem levels are necessary.

In our study, we examined the planktonic community organization from viruses to zooplankton and fluxes of major elements in a small humic waterbody with strong aquatic–terrestrial coupling. Our high-resolution in situ sampling campaign during two EREs and subsequent quasi-stable conditions addressed two fundamental questions: (1) What are the patterns in responses of community compositions and functions to ERE disturbances? (2) What is the role of small headwater ecosystems in the transfer of nutrients and organisms across the terrestrial–freshwater interface at stable and extreme hydrological conditions?

We could show that although EREs caused a washout of resident organisms, the microbial community recovered within 2 weeks in four well-defined succession phases. While reassembly of phytoplankton and especially zooplankton took considerably longer, dissolved nutrients, which were enhanced hundredfold during the floods, returned to predisturbance levels even earlier than microbes. We thus conclude that the microbial community and the aquatic food web, in general, are not resistant but extraordinarily resilient and show a fast recovery after disturbances by EREs. Moreover, the microbial community of small ponds plays a crucial role as a sink of soluble nutrients that would otherwise be exported to the downstream network. This important function was rapidly re-established after the EREs, emphasizing high functional resilience of such ecosystems.

Results and discussion

General overview of sampling site and hydrological conditions

For our study, we have chosen Jiřická Pond24,25 (Supplementary Fig. 1) located in the headwaters of Malše River (Czech Republic). It is a small humic waterbody strongly influenced by forested landscape representing a freshwater ecosystem typical for Central and Northern Europe (Supplementary Note). The pond was studied during a highly dynamic hydrological period controlled by precipitation (5 May to 27 June 2014). Pond epilimnion and main tributary stream (500 m upstream from the pond) were sampled three times per week to unveil ERE-associated changes in the plankton community and to estimate influence of surrounding landscape on organismic dynamics in the pond. Sampling started 1 month after the snow-melt at low inflow rates (0.04–0.27 m3 s−1) resulting in a water retention time in the pond surface layer between 6.7 and 11.1 d. The precipitation rates did not exceed 11.3 mm d−1 before 16 May, when continuous rainfall for 3 d introduced 79.5 mm of water and caused the first flood event in the pond watershed (ERE1). Another heavy rainfall yielded precipitation of 112.3 mm and resulted in the second ERE (ERE2) equivalent to a 5-yr flood around 28 May. Both EREs severely decreased the water retention time of the pond surface layer to 0.46 d (ERE1) and 0.44 d (ERE2) and increased flow rates of the stream to 2.2 m3 s−1 (ERE1) and 2.5 m3 s−1 (ERE2) (Fig. 1a,b). Hydrological changes strongly affected physical and chemical parameters at both sampling sites (Supplementary Note); for example, seasonally rising water temperatures and stable pH values were interrupted by two pronounced ERE-related drops (Fig. 1c,d).

Fig. 1: Hydrological, physical, chemical and microbial parameters measured in the Jiřícká watershed.
figure 1

a,b, Precipitation (grey shade) in the watershed of Jiřícká pond, stream flow rate (a) and water retention time at 0.5 m depth in the pond (b). c,d, Water temperature and pH in stream (c) and pond (d). e,f, Abundances of prokaryotic cells and virus-like particles in stream (e) and pond (f). Prokaryotic abundances counted from two independent biological samples are presented as circles. The dashed lines indicate two EREs observed during the sampling campaign. Parameters measured in the stream are shown in orange (a,c,e) and those measured in the pond are shown in blue (b,d,f).

Source data

Full size image

Flood-caused microbial dynamics in the pond

Dynamics of bacteria, which represented an abundant group of plankton community even during and shortly after EREs, was studied with two approaches: (1) 16S ribosomal RNA gene amplicon sequencing to assess bacterial community diversity and (2) catalysed reporter deposition in situ fluorescence hybridization (CARD-FISH) for accurate enumeration of the 20 most abundant bacterial groups (Supplementary Table 1). Distance matrices calculated for data obtained with these methods displayed high agreement (R2 = 0.63, P < 0.0001, one-sided Mantel test) and revealed that bacterial community composition underwent pronounced shifts after both EREs but rapidly cycled back to the predisturbed stage during subsequent developments (Fig. 2). According to agglomerative clustering analysis (Extended Data Fig. 1a), four distinct phases could be distinguished in this succession pattern, which in their sequence and characteristics followed the course of an ‘adaptive cycle’ (Fig. 2a) that was repeated twice during our entire high-frequency sampling (Fig. 2b). This ‘adaptive cycle’, a concept initially developed to describe community organization and dynamics in complex ecosystems23, consists of four well-defined phases (Fig. 2a). Two of those, so-called ‘r’ (growth and exploitation) and ‘K’ (conservation) phases, were adopted from succession ecology26 and represent the slow fore-loop of the cycle characterized by biomass accumulation and an increase in connectedness within the ecosystem. The phases ‘Ω’ (collapse and release) and ‘α’ (reorganization) are parts of the fast back-loop of the cycle that results either in a reorganization of the system and escape towards an alternative steady state or in the repetition of the previous cycle. Distinguished twice-repeated phases were associated with ranges of water retention time and their separation could be additionally confirmed with similarity analysis (ANOSIM, R = 0.95, P = 0.0001). Beside hydrological characteristics, the phases were also consistent with the trends in several chemical parameters and food web organization and are described in detail next.

Fig. 2: Bacterial succession in the pond following the adaptive cycle.
figure 2

a, Representation of the adaptive cycle; α- and Ω-phases comprising the back-loop, r- and K-phases represent the fore-loop of the cycle. X indicates the path to an alternative development. Panel adapted with permission from ref. 23, Island Press. b, Multidimensional scaling on the basis of 16S rRNA gene amplicons from pond bacterial communities (Kruskal’s stress: 0.17; ANOSIM test for the four phases on one-sided 95% CI: R = 0.95, P = 0.0001). A line connects samples following chronological order. The diameters of the circles are proportional to the Pielou’s evenness of the samples (varying between 0.64 and 0.94). Pie charts depict the bacterial diversity at the class–phylum level. Letters and associated colours visualize the phases of observed succession. *Gammaproteobacteria are represented without Betaproteobacteriales, which are shown separately. ‘K’ and violet colour, conservation phase; ‘Ω’ and orange, collapse and release phase; ‘α’ and light green, reorganization phase; ‘r’ and aquamarine, exploitation phase. c–f, Radial plots depicting percentages of bacterial groups obtained by CARD-FISH (in percentage of hybridized bacteria). Radial axes represent percentages. Angular axes represent time. Colours of segments and symbols correspond to succession phases. The time series start with circles from violet (first K-phase) to aquamarine (first r-phase). Subsequent samples in triangles represent the second Ω-phase (orange), second α-phase (light green), second r-phase (aquamarine) and final K-phase (violet). This graph design enables easy comparisons between samples taken in the same phase. c, Relative abundances of general bacterial groups. d, Relative abundances of selected Actinobacteria groups. e, Relative abundances of Verrucomicrobia groups. f, Relative abundances of selected Betaproteobacteriales groups. Ca., Candidatus.

Source data

Full size image

Conservation phase (K) was observed at the beginning and the end of the sampling campaign, that is the first six and last six samplings (Fig. 2b and Extended Data Fig. 1a). During these periods, the pond was characterized by a water retention time longer than 3.9 d and high consumption rates of soluble inorganic nutrients (Fig. 3). Up to 60% of the dissolved reactive phosphorus (DRP) and up to 100% of nitrate and ammonia transported by the inlet were consumed during this period in the pond indicating that plankton growth was severely nitrogen limited (Extended Data Fig. 2). As expected, a relatively low bacterial diversity and community evenness (Supplementary Fig. 2) were characteristic for these resource-limited conditions23. According to the isolation source information of the closest relatives, most amplicon sequence variants (ASVs) detected during the K-phase (65.5%) were associated with aquatic environments and consisted of a few dominant taxa typical for humic lakes27 (Actinobacteria, Verrucomicrobia, Bacteroidetes and Betaproteobacteriales), which also dominated bacterial abundance derived by the CARD-FISH analyses (Figs. 2c and 4a). Highest proportions of K-phase specific ASVs (Extended Data Figs. 3 and 4) were affiliated with ‘Candidatus Nanopelagicales’ (Actinobacteria) and the family Puniceicoccaceae (Verrucomicrobia, Opitutae), which have small-sized cells and genomes28,29, and can be considered as typical K-strategists. According to the r/K-selection ecological concept, those are slow-growing organisms, for which a stable environment selects due to limited resources30,31. During the K-phase, the organization of the microbial food web, phytoplankton and higher trophic levels in the pond was most complex (Figs. 4 and 5), corresponding to the expected high connectedness in the system23. This resulted in the highest bacterial loss rate due to increased protistan bacterivory indicating a severe top-down control of the bacterial community (Extended Data Fig. 7) and provided conditions favourable for grazing-resistant microbes, for example ‘Candidatus Nanopelagicales’ (Fig. 2d), probably protected from protistan grazers due to their small cell size and Gram-positive cell wall structure28,32. Thus, both top-down control and K-selection contributed to ‘species sorting’ by habitat conditions33, which dominated microbial community assembly during the K-phase.

Fig. 3: Transfer rates of selected dissolved nutrients in Jiřícká pond.
figure 3

a,b, Daily import and export rates of DRP (a) and nitrate (b). Import to the pond is depicted in orange, export from the pond is depicted in blue. Dotted lines are separating three equal periods (23 d each) used for the calculation of exports and imports at different hydrological conditions. Black triangles indicate EREs.

Source data

Full size image
Fig. 4: Absolute abundances or biovolumes of selected plankton groups in the pond.
figure 4

a–d, Abundances of different bacterial groups obtained by CARD-FISH analyses. a, Abundances of general bacterial groups. b, Abundances of specific Betaproteobacteriales groups. c, Abundances of specific Actinobacteria groups. d, Abundances of specific Verrucomicrobia groups. e, Total chlorophyll a concentration and biomass of different phytoplankton taxa. f, Abundances of HNF and ciliates. Abundances counted from two independent biological samples are presented as circles. g, Abundances of rotifers and crustaceans. Panel adapted with permission from ref. 24, PAGEPress, under a Creative Commons licence CC BY-NC 4.0. Violet, ‘K’ conservation phase; orange, ‘Ω’ collapse and release phase; light green, ‘α’ reorganization phase; and aquamarine, ‘r’ exploitation phase.

Source data

Full size image
Fig. 5: Food web structure of the Jiřická Pond plankton community during flood-caused cyclic succession.
figure 5

Different phases are indicated in an external circle with colours: violet, ‘K’ conservation phase; orange, ‘Ω’ collapse and release phase; light green, ‘α’ reorganization phase; and aquamarine, ‘r’ exploitation phase. Relative abundances (where ‘1’ corresponds to maximum abundances) of different plankton groups are depicted in the four inner circles (from outside to inside): Zooplankton is represented by rotifers (dominated by Polyarthra sp.) and crustaceans (dominated by Eubosmina sp.). Protists are represented by HNF and ciliates. Food preferences of zooplankton and protists at different phases are indicated by bubbles. Phytoplankton is represented by summarized dynamics of all algal species. Bacterioplankton is divided into three categories: emigrants, r-strategists and all-rounders, and K-strategists. Nutrient gradient based on nitrate concentrations (0.00–0.14 mg N l−1) is displayed by different background colours; stream flow rates are indicated by different sized arrows.

Full size image

Mass effects-driven collapse and release phase (Ω) consists of two samples collected during the EREs (Extended Data Fig. 1a). They were characterized by short water retention time (<0.46 d) and up to twofold increased concentrations of dissolved organic carbon (DOC) and dissolved nutrients (Fig. 3 and Extended Data Figs. 2 and 8). The phase displayed low prokaryotic densities but yielded the highest diversity and evenness of the bacterial community (Supplementary Fig. 2), resembling flood-associated microbial dynamics in subsurface ponds34,35 or post-typhoon recovery pattern36. The absolute and relative abundances of bacterial groups dominating in the previous phase collapsed, except for Verrucomicrobia and Alphaproteobacteria (Figs. 2c and 4a,d). The proportion of ASVs, which had closest relatives reported from aquatic habitats, decreased to less than half (18%), with only 3% being associated with lake bacteria. A high proportion of soil-related ASVs (46%) in the pond indicated a tight terrestrial–aquatic coupling and an important role of ‘mass effects’, a term used in metacommunity theory for emigration processes33. Emigrant phylotypes (from the terrestrial surrounding) within Verrucomicrobia, Alphaproteobacteria, Acidobacteria, Solirubrobacteraceae, Omnitrophicaeota and Patescibacteria displayed short-lived peaks but did not persist in the epilimnion during the subsequent stabilizing hydrological conditions, revealing high ‘colonization resistance’37 of the pond (Extended Data Fig. 5a).

Virus-like particles in both pond and stream had local maxima during EREs (Fig. 1e,f), which were most likely connected to the introduction of virus-like particles from soil and/or their resuspension from the stream sediment. In contrast, larger plankton groups (protists, phytoplankton and zooplankton) were almost entirely washed-out from the pond (Fig. 4e–g).

Reorganization phase (α) comprises three samples collected between the third and fifth day after the ERE peaks. This phase was characterized by a still low water retention time in the pond (<2.0 d) and enhanced DOC and nutrients concentrations (Fig. 1 and Extended Data Fig. 8). Instead of soil-associated emigrant ASVs detected in Ω-phase, we observed high read numbers of ASVs associated with the stream (Extended Data Fig. 5b), for example representatives of Duganella, Rhodoferax, Flavobacterium and Pseudomonas lineages. At the same time, we detected several pond-specific groups, for example representatives of Polynucleobacter cluster C38 (PnecC), ‘Candidatus Methylosemipumilus turicensis’39, as well as a pond-specific ASV of Limnohabitans lineage LimA (Fig. 2f and Extended Data Fig. 6a,b). The observed occurrence patterns exemplify an evident reorganization in assembly mechanisms of the pond bacterial community, namely a decrease of mass effects and the appearance of species sorting selecting during this phase for so-called r-strategists, fast-growing organisms able to rapidly exploit resources and to colonize unstable environments30,31,40. Epilimnetic populations of protists, phytoplankton and zooplankton still remained at negligible levels.

Exploitation phase (r) contains seven samples (two after ERE1 and five after ERE2). In the course of this stage, chemical and hydrological characteristics, as well as consumption rates of soluble inorganic nutrients changed to levels comparable to stable conditions due to increasing exploitation. Linkages between bacterial communities of pond and stream were reduced (Extended Data Figs. 5 and 6); as with increasing water retention time, stream relevant ASVs were replaced by pond-specific ones. Betaproteobacteriales displayed maxima in CARD-FISH and sequence read proportions accounting for up to 57.5% of all bacteria (Fig. 2b,c). Among them, several ASVs were mainly restricted to the exploitation phase: the majority of Limnohabitans lineage LimA, previously reported as EREs associated copiotrophs41, seemed to represent typical r-strategists; while an ASV of PnecA (Polynucleobacter rarus) could be favoured by a decrease in pH (ref. 42). On the other hand, PnecC and ‘Candidatus Methylosemipumilus turicensis’ showed persistent read and cell count proportions also in the α- and K-phases (Fig. 2f and Extended Data Fig. 6b). Both taxa are small- to medium-sized bacteria with passive lifestyles based on the use of photodegradation products or C1 compounds39,43. This seemingly successful strategy renders them to abundant all-rounders in the described ecosystem. During the r-phase, also fast-growing algae and heterotrophic nanoflagellates (HNF) attained their maxima, while ciliates, rotifers and crustaceans displayed their peaks later24 (Figs. 4e,g and 5). Uncontrolled growth of HNF together with nutrient limitation resulted in the collapse of bacterial r-strategists, which thereafter were replaced by grazing-resistant K-strategists, for example ‘Candidatus Nanopelagicales’ and Bacteroidetes28,40,44. This transit from r- to K-phase completed the observed succession cycle.

The ERE-induced succession of the entire planktonic community

The bacterial community in the pond revealed an unexpected resilience in response to EREs and returned to the initial steady state through the above-described well-determined pathway within 16 d after ERE2. The 50 most common pond ASVs displayed highly repetitive appearance patterns in the two loops and were mainly restricted to one or two distinct phases. Though, we missed the recovery of three ASVs well represented during the first conservation phase (one representative of Verrucomicrobia and two of Bacteroidetes, Extended Data Fig. 3), not even a single emigrant ASV could establish successfully in the pond after EREs. The succession of pond eukaryotes generally followed patterns outlined in the PEG (Plankton Ecology Group) model15,45 tightly related to the organisms’ size and growth strategy. Fast-growing bacterial r-strategists were followed by bacterivorous HNF, having comparable growth potential46 as these prey bacteria. In parallel, we observed a rapid increase of small green flagellates (Chlamydomonas spp.; Extended Data Fig. 7c). With a short delay of 7 d, the growth of these edible algae, as well as HNF, became top-down controlled by the subsequently recovered ciliate community (Figs. 4 and 5 and Extended Data Fig. 7), represented mainly by predatory prostome ciliates47. After the decline of edible HNF and algal species, more omnivorous ciliates were detected in the pond. This led to an increase in ciliate bacterivory during the second conservation phase (Extended Data Fig. 7d) and higher loss rates of Verrucomicrobia and Bacteroidetes (Fig. 4a,d), which appear to be less susceptible to HNF grazing but highly vulnerable to omnivorous–bacterivorous ciliates (for example, Halteria sp. and Rimostrombidium spp.). Similar to bacteria, eukaryotic dynamics were driven by resource-limitations and disturbances but also by top-down control via food web structure, dominated at different successional stages by HNF, ciliates, crustaceans or copepods (Figs. 4 and 5 and Supplementary Fig. 3). As expected for dynamic waterbodies, zooplankton in the pond was greatly dominated by rotifers (average ratio to crustaceans 126:1 (refs. 24,48); Fig. 4 and Supplementary Fig. 3). After EREs, reassembled populations of cladocerans and rotifers displayed a higher diversity but did not reach the predisturbance densities because of the establishment of a predatory copepod population (Mesocyclops leuckarti24; Supplementary Note). Recovery of this species in early summer represents a well-known seasonal development49. Temporal patterns were also observed in the phytoplankton dynamics; for example, during the second conservation phase, Chlamydomonas spp. persisted in the phytoplankton with Cryptophytes and Chrysophytes, typical representatives of spring assemblages32 (Fig. 4e and Supplementary Fig. 4). Seasonal dynamics of phytoplankton and zooplankton seemed to be accompanied by stochastic processes in the post-EREs reassembly (back-loop of the adaptive cycle) that resulted in alternations in their composition (Supplementary Figs. 3 and 4). Similar pattern, predictable recovery of bacterial community and striking differences in the recovery of phytoplankton, was observed in lake communities after typhoon disturbances36. Surprisingly, the reorganization of higher trophic levels from ciliates to cladocerans showed a limited influence on bacterial community composition observed in both K-phases.

Segregation between stream and pond

Most physical and chemical characteristics in both stream and pond displayed comparable dynamics tightly linked to EREs (Extended Data Fig. 2 and 8). Significant differences between the stream and pond were observed only for temperature, which was higher in the pond (Mann–Whitney rank test: P < 0.001) and soluble fractions of nitrogen and phosphorus, which were higher in the stream (Mann–Whitney rank test: P < 0.001 for dissolved phosphorus (DP) and nitrate; P = 0.022 for ammonia; two-tailed t-test: t(46) = −0.8, P = 3 × 10−9 for DRP). Despite similar chemistry, constantly low water retention time did not allow the establishment of a stable plankton community in the stream. In result, we could not identify typical stream eukaryotic plankton populations due to their low densities or recognize clear succession patterns caused by EREs even in bacterial populations (Extended Data Fig. 10). High mass effects resulted in 26% of soil-associated bacterial ASVs in the stream even during a low hydrological regime. Despite high stochasticity in the community assembly, high nutrient contents and low grazing pressure provided favourable conditions for r-strategists in the stream. Such ASVs were detectable in the stream at all non-EREs conditions but in the pond exclusively during the reorganization period that increased the similarity between both sides (Extended Data Figs. 1c, 5b and 9). However, a water retention time of 2 d was sufficient to separate inlet and pond microbial communities and to prevent an overlap of the most abundant phylotypes (Extended Data Figs. 1 and 9). This notion demonstrates the compositional discontinuity between the hydrologically connected stream and pond ecosystems and limited influence of mass effects on pond community assembly during stable hydrological conditions, that is the K-phase.

Transfer of nutrients through the pond

Low initial flow rates (0.07 ± 0.03 m3 s−1; average ± s.d.) at the beginning of the sampling campaign (23 April to 16 May 2014) were associated with very low loads of dissolved nutrients from the watershed (Extended Data Fig. 2). The average daily export of nitrate and DRP from the pond amounted to 0.05 ± 0.07 kg and 0.04 ± 0.03 kg, respectively. Altogether, during this period, 1.05 kg of nitrate N (8% of the import) and 0.96 kg DRP (48% of the import) were exported from the system. Loads of soluble nutrients increased dramatically with rain events. Daily export of nitrate and DRP reached up to 22.6 kg and 2.6 kg, respectively. In 23 d, covering both EREs and community recovery periods, 103.7 kg of nitrate N and 14.6 kg of DRP were released to the downstream network. The subsequent equivalent period was characterized by slightly elevated flow rates and loads from the watershed (0.13 ± 0.03 m3 s−1; 31.8 kg nitrate N; 4.1 kg DRP) but also by an already recovered microbial community and hence low nutrient export (0.7 kg nitrate N; 1.5 kg DRP). This indicates that ponds substantially control the nutrient release from the watershed, as nitrate export remained in the same range during pre- and post-EREs’ periods, while nitrogen consumption in the pond differed more than three times following EREs (Extended Data Fig. 2).

Conclusions

In this study, we present the detailed mechanism underlying functional and compositional resilience of the small pond to EREs and connect the flow of nutrients, bacterioplankton dynamics and reorganization of higher trophic levels to the entire ecosystem development. The backbone of community recovery is represented by the twice-repeated adaptive cycle in bacterioplankton composition. Four well-defined phases of this cycle reflected the transition from the dominance of mass effects during and shortly after the EREs to the prevalence of species sorting in community assembly mechanisms during conservation phases. Within this successional development, we report a remarkable fast return of nutrient export rates to predisturbance levels in just 9 d, rapid reestablishment of the prokaryotic and eukaryotic microbial community within 16 d and slightly longer reassembly of phytoplankton and zooplankton, which was to a greater degree influenced by the seasonal succession of different taxa. Thus, our high-resolution sampling campaign revealed that small waterbodies have the potential of surprisingly fast and almost complete compositional and functional recoveries after ERE disturbances. However, as the frequency of EREs may increase in the future, it remains questionable to which extent this resilience can be challenged and how possible alternations in the composition could affect ecosystem function and provided services. Knowledge on functional and compositional resilience of small waterbodies, representing the most abundant freshwater habitats of Earth4, is of high importance for sustainable management and for predicting the impact of future climate change on their role in transformation and control of organic matter and nutrient flow through the landscape.

Methods

Study site

Jiřická Pond (Pohořský Pond) is located in Novohradské Mountains at an elevation of 892 m above sea level (48.6189° N, 14.6731° E) and was formed by damming of the Pohořský stream in the eighteenth century. Jiřická is a small (0.04 km2), shallow (4 m maximum depth) humic pond; the watershed (12.5 km2) is composed to a high extent of forest but also of pastures and raised peat bogs24. The intensive sampling of the pond and the main tributary of Jiřická (stream) took place between 5 May and 27 June 2014, sampling three times per week (Monday, Wednesday and Friday; total samples 24). Samples for chemical and biological analyses were taken with a Friedinger sampler at 0.5 m depth at the deepest point of the pond. The stream was sampled 500 m upstream from the pond.

Meteorological and hydrological characteristics

Meteorological parameters, including the level of precipitation, were monitored with the meteorological station M4016 (Fiedler-Mágr, Czech Republic) installed near the pond during the study period. The water level of the inlet was measured three times per week and the flow rates were calculated on the basis of the relation between water level and discharge calibrated at different flow regimes using handheld Acoustic Doppler Velocimeter (FlowTracker). Water age at 0.5 m was calculated using the CE-QUAL-W2 model of reservoir hydrodynamics50.

Enumeration of planktonic organisms

Samples for microbial analyses (2 l) were split into subsamples that were either immediately fixed on site or transported to the laboratory within 2 h. For enumeration of virus-like particles, 1 ml of water samples from both sites was fixed with glutaraldehyde (1% final concentration) for 10 min, shock-frozen in liquid nitrogen and transported in liquid nitrogen to the laboratory, where samples were stored at −80 °C. Enumerations of virus-like particles were done using an inFlux V-GS cell sorter (Becton Dickinson) as previously described51. For bacterial enumeration, water subsamples were fixed with formaldehyde on site (2% final concentration). In the laboratory, duplicates of 1–2 ml were stained with DAPI and filtered onto black 0.2-μm pore-size polycarbonate membranes. The bacterial abundance was determined via epifluorescence microscopy52. Numbers of HNF and ciliates were estimated in a similar manner from Lugol–formaldehyde–thiosulfate fixed subsamples53. For this purpose, 5–30 ml were stained with DAPI and filtered onto a black 1.0-μm pore-size membrane. Protozoan bacterivory was estimated using fluorescently labelled bacteria (FLB)54 prepared from a mixture of cultures of small rods of Limnohabitans lineage LimC46,55 and two undescribed strains of Polynucleobacter lineage PnecC. The sizes of FLB reassembled the typical size class distribution of the bacterioplankton in the pond. HNF and ciliate abundances and FLB uptake rates were determined in short-term FLB direct-uptake experiments, with ciliates being determined to genus level56. To estimate total protozoan grazing, we multiplied the average uptake rates of HNF and ciliates with their in situ abundance. Phytoplankton samples were preserved with acid Lugol solution and stored in a dark box at room temperature. Species were enumerated with the Utermöhl method57 using the microscope Olympus IMT2. The mean algal cell dimensions for biovolume calculation were obtained using the approximation of cell morphology to regular geometric shapes58.

Zooplankton was sampled once a week. Crustaceans were sampled by vertical hauls using an Apstein plankton net (200-μm mesh) from the entire water column59 (4 m). Rotifers were sampled from the whole water column using a plastic tube of the appropriate length. A total of 55 l of water was subsequently filtered using a 35-μm net. The zooplankton samples were preserved in 4% formaldehyde and species abundances were determined microscopically.

Nucleic acids sampling, isolation and processing

Microbial biomass was collected on Sterivex filter units (polyethersulfon membranes, 0.22-µm pore size, Millipore) by filtering approximately 600 ml with a sterile syringe (n = 48). Additional samples for RNA isolation were prepared in the same manner every Wednesday (n = 14). The units were sealed, flash-frozen in liquid nitrogen and transported in liquid nitrogen to the laboratory and stored at −80 °C until further processing. Nucleic acids (DNA and RNA) were isolated using phenol–chloroform–isoamyl alcohol extraction method according to a previously described protocol60. After washing with 70% ethanol, nucleic acids were briefly dried, dissolved in 20 µl of nuclease-free water and stored at −80 °C. For complementary DNA synthesis, DNA was degraded in the DNA/RNA extract by twice-repeated incubation with Turbo DNase at 37 °C for 45 min, followed by addition of DNase inactivation reagent according to the Turbo DNA-free Kit protocol (Ambion). Thereafter, 20–50 ng of RNA were mixed with 1.5 µl of primer mix (Random Primers 6, New England BioLabs) in RNA free water (final volume 10.5 µl) and incubated for 5 min at 70 °C. Reverse transcription was performed with the Array Script enzyme (Ambion) according to the manufacturer’s instructions at 42 °C for 2 h and terminated by enzyme deactivation at 95 °C for 5 min.

Amplicon sequencing and analysis

Primers F27 and R519 were used to amplify the V1–V3 region of the 16S rRNA gene. A single-step PCR using HotStarTaq Plus Master Mix Kit (Qiagen) was conducted as follows: 94 °C for 3 min, followed by 28 cycles of 94 °C for 30 s, 53 °C for 40 s and 72 °C for 1 min; after which a final elongation step at 72 °C for 5 min was performed. After PCR, all amplicon products were mixed in equal concentrations and purified using Ampure beads (Agencourt Bioscience). The amplicons were sequenced using the Roche 454 GS FLX Titanium platform at MR DNA laboratory (Shallowater) with an average sequencing depth of 50,000 raw reads per sample.

Demultiplexed reads were processed using DADA2 package v.1.12 in R61,62. The primer free sequences were trimmed to 350 base pairs length according to the quality report. The values recommended for 454 data were used for the denoising, namely homopolymer gap penalty ‘−1’ and increased to 32 band size. Chimeric ASVs were detected and removed with de novo algorithm. The remaining unique ASVs were taxonomically assigned using naive Bayesian classifier method and trained Silva SSU database release 132 (refs. 63,64). Variants affiliated to mitochondria and plastids were excluded. Two samples with read numbers <10,000 were removed from analysis (stream DNA from 25 June 2014 and stream cDNA from 11 June 2014). The final number of reads in samples ranged from 19,390 to 87,185; the dataset was rarefied to 19,390 reads per sample before diversity and statistical analyses.

From the pond samples, the 50 most common ASVs were selected according to the percentages of reads per sample for a graphical representation of read abundance dynamics. Since these phylotypes were not detected in the EREs samples, the same procedure was done additionally for most frequent ASVs in the EREs samples (>0.3% of read number in the sample; n = 14). The next relatives of these sequences and of the 50 most common ASVs from the stream samples were assigned with BLAST 2.2.21 against Silva SSU database release 132 (refs. 63,65) and were used for constructing a bootstrapped Randomized Axelerated Maximum Likelihood tree (GTR-GAMMA model, 100 bootstraps66). The isolation source information for the next relatives with identity scores >97% was collected from NCBI using reutils v.0.2.3 package67 in R62 and manually classified to water- versus soil-related environments and others.

Probe design for CARD-FISH

Clone libraries of 23S rRNA gene were constructed for several dates to recruit the sequences affiliated with bacterial taxa of interest (Betaproteobacteriales and Actinobacteria). Clones were sequenced with five primers targeting 23S rRNA genes (Supplementary Note and Supplementary Table 2), assembled with Geneious R9.1 and aligned with the SINA web-aligner68. Afterwards, sequences were imported to Silva LSU database release 132 using ARB software (v.arb-6.0.6)69 for tree reconstruction and CARD-FISH probe design. Alignments of 68 almost full-length sequences (>2,000 nucleotides) and their closest relatives were manually improved and two Randomized Axelerated Maximum Likelihood trees (GTR-GAMMA model, 100 bootstraps)64 were constructed for Limnohabitans spp. (Supplementary Fig. 5) and Actinobacteria affiliated with ‘Candidatus Nanopelagicales’ and the luna2 lineage (Supplementary Fig. 6), respectively. Probe design for the LimB cluster of Limnohabitans spp. and the luna2 cluster of Actinobacteria was done with the tools probe_design and probe_check. Three competitor probes were designed for probe LimB-23S-920 and two competitor probes and three helper oligonucleotides were designed for probe luna2-23S-1239 to further increase mismatch discrimination and rRNA accessibility70. Additionally to 23S rRNA targeted probes, we also designed and modified a set of probes targeting the 16S rRNA genes of Verrucomicrobia as an online check with SILVA64 revealed that the general probe ver47 (ref. 71) targets only a part of freshwater Opitutales (79%; Supplementary Fig. 7). We modified this probe in a way that it discriminated against this order (probe ver46: 0.1% coverage of Opitutales) and designed an additional probe that targets 97.1% of all available sequences of Opitutales (probe opitu-346; Supplementary Fig. 7). All probes were checked in silico for specificity and coverage using the TestProbe function of SILVA64. Probes were further tested in silico with the mathfish online tool72 and in the laboratory with different formamide concentrations in the hybridization buffer until stringent conditions were achieved (Supplementary Table 1).

CARD-FISH analysis

For CARD-FISH analyses, formaldehyde-fixed samples from the pond epilimnion (n = 24) were filtered within 10 h after fixation onto white polycarbonate membrane with 0.2-µm pore size (47 mm diameter, Millipore). Filters were stored at −20 °C before subsequent analysis. The filter sections were used for specific hybridization of 24 bacterial groups with horseradish peroxidase labelled oligonucleotide-probes (Supplementary Table 1) and subsequent amplification with fluorescein labelled tyramides as previously described73. Proportions of CARD-FISH positive signals were evaluated with an epifluorescence microscope (Olympus BX-53F) for at least 500 DAPI-stained cells per probe.

Physical and chemical characteristics

Water temperature was monitored during the entire sampling period at the pond and stream sampling sites at hourly intervals with temperature data loggers (Tidbit) and averaged for daily values. Samples for other physical and chemical analyses were collected from pond and stream sites at every sampling in prewashed and rinsed polyethylene terephthalate sampling bottles and were delivered to the laboratory in thermos-boxes within 2 h. Chlorophyll a concentration was determined spectrophotometrically after extraction with acetone74. Acid-neutralizing capacity (determined by Gran titration) and pH were measured by Radiometer TIM 865 station (Hach). Conductivity was evaluated using inoLab Cond 720 (WTW). Aliquots of samples were filtered through 0.45-µm nylon membrane filters (Fisher Scientific) for ion chromatography and silica determination. Dissolved reactive silica was determined by the molybdate method75. Concentrations of cations (Na+, K+, Ca2+, Mg2+, NH4+) and anions (Cl–, SO42–, F–, NO2–, NO3–) were determined by ion chromatography (Dionex IC25). All following chemical analyses, except for total phosphorus, were performed after filtration of water samples through 0.4-µm glass fibre filter (GF-5, Macherey-Nagel). DRP was determined by the molybdate method76 and total TP and DP with a modification of this method with sample preconcentration and perchloric acid digestion77 yielding detection limits of 1 and 0.5 µg l−1, respectively. DOC was analysed as a non-purgeable organic carbon by catalytic combustion at 680 °C (Elementar), the sample was acidified with HCl (1 M) to pH <4 and air purged for 3 min before the analysis. The detection limit was ~0.05 mg l−1. Dissolved nitrogen concentrations were obtained using a vario TOC cube (Elementar), nitrogen was determined as NO by infrared (NDIR) detector below the ppb level. Additionally, molecular weight distribution and detection of low molecular weight organic acids were performed according to protocol described in the Supplementary Note.

Data analyses

Richness, abundance-based coverage estimator (ACE), diversity estimator Chao1, Shannon index and Pielou’s evenness were calculated for rarefied data of 16S rRNA gene amplicons with the R package vegan78. To test the positive correlation between 16S rDNA amplicon data and CARD-FISH based relative abundances, Bray–Curtis dissimilarity matrices were computed and a one-tailed Mantel test (Pearson) was performed on the 95% confidence interval using XLSTAT 14 (Addinsoft). To study the community dynamics patterns caused by EREs, the agglomerative hierarchical clustering (AHC) based on Bray–Curtis dissimilarity was calculated for 16S rDNA data from pond and stream samples together or separately using the complete linkage algorithm and automated truncation settings (on the basis of entropy) in XLSTAT 14 (Addinsoft). Grouping of pond samples on the basis of the clustering analysis was associated with water retention time ranges and was additionally tested with ANOSIM statistics using 9,999 permutations in the R package vegan78. Multidimensional scaling plots for the visualization of distances of pond and stream samples were done on the basis of Bray–Curtis dissimilarity matrices using XLSTAT 14 (Addinsoft). To compare physical and chemical parameters in the stream and pond, t-tests or alternatively Mann–Whitney rank tests (for non-normally distributed data, for example TP, DP, ammonia, nitrate, nitrite, pH and temperature) were performed on the 95% confidence interval with SigmaPlot 12.5 after applying the normality Shapiro–Wilk tests. The nutrient balances represent the differences between values obtained from a transfer model (theoretical values calculated for the nutrient transport to epilimnion on the basis of nutrient concentrations measured in the stream and hydrological parameters) and values measured in the pond epilimnion. The transfer of different compounds from the inlet to the epilimnion of Jiřická Pond was modelled using the hydrodynamic and water quality model CE-QUAL-W2 (ref. 48). For modelling purposes, the waterbody was separated into 12 segments and all calculations were done for 0.5 m depth layers48. The import and export of nutrients were calculated by multiplying the average daily flow data by the nutrient concentrations.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.