Aquatic food webs in deep temperate lakes: Key species establish through their autecological versatility

Microbial planktonic communities are the basis of food webs in aquatic ecosystems since they contribute substantially to primary production and nutrient recycling. Network analyses of DNA metabarcoding data sets emerged as a powerful tool to untangle the complex ecological relationships among the key players in food webs. In this study, we evaluated co‐occurrence networks constructed from time‐series metabarcoding data sets (12 months, biweekly sampling) of protistan plankton communities in surface layers (epilimnion) and bottom waters (hypolimnion) of two temperate deep lakes, Lake Mondsee (Austria) and Lake Zurich (Switzerland). Lake Zurich plankton communities were less tightly connected, more fragmented and had a higher susceptibility to a species extinction scenario compared to Lake Mondsee communities. We interpret these results as a lower robustness of Lake Zurich protistan plankton to environmental stressors, especially stressors resulting from climate change. In all networks, the phylum Ciliophora contributed the highest number of nodes, among them several in key positions of the networks. Associations in ciliate‐specific subnetworks resembled autecological species‐specific traits that indicate adaptions to specific environmental conditions. We demonstrate the strength of co‐occurrence network analyses to deepen our understanding of plankton community dynamics in lakes and indicate biotic relationships, which resulted in new hypotheses that may guide future research in climate‐stressed ecosystems.


| INTRODUC TI ON
Aquatic food webs are the basis for ecosystem functioning (Azam, 1998;Sommer et al., 2012). Relationships among involved organisms form extraordinarily complex networks in which protists hold a key position and together with prokaryotes (including phototrophic forms such as cyanobacteria) support aquatic ecosystem resilience and stability. For example, phytoplankton and their exudates are an important food and nutrient source for heterotrophic consumers and bacteria, respectively (e.g. Saleem et al., 2013;Šimek et al., 2014;Zingel et al., 2007). Additionally, algal and ciliate communities are valuable indicators for the trophic status and health of lakes and running waters (e.g. Foissner et al., 1991Foissner et al., , 1999Smol & Stoermer, 2010). Ciliates and heterotrophic flagellates are efficient grazers controlling bacterial and phytoplankton biomass and thus strikingly influence major biogeochemical cycles in aquatic habitats (e.g. Saleem et al., 2012;Sanders et al., 1989;Šimek et al., 2019;Weisse & Müller, 1998;Zingel et al., 2006). At the same time, they are a significant carbon source for metazoans such as rotifers, which again are the nutritional basis for larger crustaceans and fish larvae (e.g. Zingel et al., 2012Zingel et al., , 2019. Understanding such complex interactions at the base of aquatic food webs is therefore a cornerstone for elucidating ecosystem functioning and functions. As climate change is rapidly changing aquatic ecosystem functioning and functions (Foley et al., 2012;Posch et al., 2012), it is of special importance and interest to enhance our understanding of the basic ecology of complex biotic interactions in different lake ecosystems.
In the first place, the inference of the complex interplay between biotic and abiotic interactions which define ecosystem functioning and functions requires a deep knowledge of the communities and factors involved. The traditional method for the identification of protistan plankton in ecological studies is microscopic identification. However, for several reasons discussed in detail previously (e.g. Caron et al., 2004), this approach has several shortcomings, especially when analyses of a high number of samples are required. Such a deep exploration of protistan plankton community structures at a high spatial and temporal resolution is for example mandatory when including climate change effects in the Plankton Ecology Group (PEG) model , which provides the theoretical framework to explain seasonal succession patterns of planktonic lake communities.
DNA metabarcoding and high-throughput sequencing (HTS) technologies contributed substantially to the toolbox of plankton ecology, compensating for the shortcomings of traditional microscopy-based analyses (Banerji et al., 2018;Mikhailov et al., 2019;Stern et al., 2018). A molecular approach allows for analyses of a high number of samples and deeper insights into complex plankton community structures. However, only the integration of specific analytical methods, such as from graph theory, guaranteed the exploitation of massive DNA metabarcoding data sets beyond descriptive diversity (Röttjers & Faust, 2018;Watson et al., 2019).
Here, we constructed and analysed networks of protistan plankton communities in the surface (epilimnion) and deep (hypolimnion) water layers of two deep temperate lakes in the Alpine European Region, namely Lake Mondsee (Austria) and Lake Zurich (Switzerland). Both lakes were strongly influenced by anthropogenic induced eutrophication during the 20th century, but were successfully restored owing to the reduction of nutrient immissions (Dokulil & Teubner, 2005). By now, each lake again reached nutrient levels similarly low to the time before pollution caused striking changes in both ecosystems. Thus, we compared two lakes that presently reflect a relatively natural trophic status, that is oligo-mesotrophic. Nevertheless, the two lakes are contrasting with regard to their water turnover (mixis) dynamics. Water turnover is extremely important for transporting nutrients from deep layers back into surface layers and, thus, has a strong influence on the seasonality of organisms and on entire food web structures. While Lake Mondsee is dimictic with two turnover events, one in autumn and one in spring, Lake Zurich is monomictic with one main mixis event in spring (Ficker et al., 2017;Yankova et al., 2017). However, global warming has led to thermal stratification patterns of water bodies in lakes of the alpine region that severely altered these turnover dynamics.
Numerous studies already documented the impact of climate change on Lake Zurich and the resulting dramatic lake effects, including economic consequences (Livingstone, 2003;North et al., 2014;Posch et al., 2012;Yankova et al., 2016Yankova et al., , 2017. This impact is highlighted by the thermal stratification of the water column, which has led to a stop of complete water turnover in spring (vernal holomixis) for several years (Yankova et al., 2017). As a consequence, nutrients accumulated in the hypolimnion and were not transported back into the epilimnion. Lake Mondsee, in contrast, does not yet show such significant ecological effects resulting from climate change although the average surface-temperature increased by almost 2°C during the last decades and leads to a prolongation of the stratification period (Ficker et al., 2017).
However, complete oxygen depletion in the hypolimnion is exceptional and complete water turnover (holomixis) reaching the deepest zones in Lake Mondsee is still ensured.
In both lakes, ciliates belong to ecologically and functionally important components of the plankton communities that are well-characterized, owing among others to their role as top predators among protists (e.g. Müller et al., 2002;Pitsch et al., 2019;Posch et al., 2015;Salbrechter & Arndt, 1994). Even though both lakes shared numerous ciliate species, we hypothesized that the ecological role and importance of individual species may differ because of the contrasting environmental characteristics (e.g. lack of holomixis and resulting nutrient availability) of both lakes. Furthermore, we assume a higher vulnerability of ecosystem structures in Lake Zurich owing to its described sensitivity to climate-induced stressors (Yankova et al., 2017).

| Sampling sites and sampling procedures
Water samples were collected biweekly from the epilimnion and the deep hypolimnion of Lake Mondsee (LM), Austria, and Lake Zurich (LZ), Switzerland, from June 2016 through May 2017. In both lakes, the epilimnion (LM epi, LZ epi) was sampled at a depth of 5 m, corresponding to the chlorophyll a (chl a) maximum (= most productive depth in both lakes). The deep hypolimnion referred to the depth close to the lake bottom but not yet anoxic (LM hypo, LZ hypo). Due to the long-term and ongoing monitoring programmes (more than 40 years) in both lakes, temporal and spatial patterns of oxygen concentrations were known and resulted in the following sampling depths. In Lake Mondsee, the deep hypolimnion referred to a depth of 40 m (max. depth of Lake Mondsee: 68 m; Dokulil et al., 2006). In Lake Zurich, the hypolimnion was sampled at a depth of 120 m (max. depth of Lake Zurich: 136 m; Bossard et al., 2001). At each sampling date, 5 L of raw water was taken from the required depth with a water sampler (Uwitec). From this sample, two times 2 L (sample duplicates, 4 L in total) of water were prefiltered through a 150-μm net gauze to remove larger zooplankton and then filtered onto a 0.65μm membrane filter (Durapore, Merck Millipore) with a peristaltic pump at a speed of 50 ml/min. Each filter was stored in a cryovial containing 1.5 ml RNAlater (Qiagen) and stored at −80°C until further processing. For the determination of bacterial abundance, we added formaldehyde (2% final concentration) as fixative to a 40 ml subsample of the raw water sample. Fixed samples were stained with SYBR-Green I (Sigma-Aldrich) in the dark and analysed via flow cytometry (Cytoflex S, Beckman Coulter). For the determination of the total algal biovolume, a 100 ml subsample of the raw water sample was preserved with neutral Lugol's solution (0.25% final concentration) and processed according to the method of Utermöhl (1958).
Environmental parameters as well as bacterial and phytoplankton cell counts throughout the sampling year are shown as supplementary information (SI 1 and SI 2).

| DNA extraction, PCR and illumina sequencing
For DNA extraction, each filter was transferred into a lysing matrix tube (Lysing Matrix E, MP Biomedicals) and 600 µl RLT buffer (Qiagen) and 6 µl β-mercaptoethanol were added. The remaining liquid in each cryovial was centrifuged and discarded. To the residual pellet, 200 µl RLT buffer and 2 µl β-mercaptoethanol were added and mixed. This mix was added to the lysing matrix tube. Each matrix tube was subjected to bead-beating for 45 s by 30 Hz (MM 200,Retsch). Total DNA was extracted using the AllPrep DNA/RNA Mini Kit (Qiagen).
The hypervariable V9 region of the 18S rDNA was amplified following a standard protocol (Stoeck et al., 2010). As forward primer, we used 1391F (5′-GTACACACCGCCCGTC-3′; Lane, 1991) and as reverse primer EukB (5′-TGATCCTTCTGCAGGTTCACCTAC-3′; Medlin et al., 1988). The PCR protocol for V9 amplification employed an initial denaturation step at 98°C for 30 s, followed by 30 cycles of 10 s at 98°C, 20 s at 61°C, 25 s at 72°C and a final five-minute extension at 72°C. The reactions were run in volumes of 50 µl using 0.5 µl Phusion polymerase (New England Biolabs (NEB)), 10 µl 5× Phusion GC buffer (NEB), 1 µl 10 mM dNTPs, 0.5 µl template DNA, 32.5 µl pure water, and 0.5 µl of each forward and reverse primer. From each DNA extract, triplicate PCR reactions were conducted to minimize PCR bias. Prior to purification (MinElute Kit, Qiagen), PCR sample replicates were pooled. Positive (eukaryotic DNA isolated from a culture of Paramecium caudatum) and negative (distilled water) controls were successfully applied during each step of the process to ensure highest quality results and exclude possible contaminations.
From the resulting PCR products, sequencing libraries were constructed using the NEB Next Ultra DNA Library Prep Kit for Illumina (NEB). Library quality was assessed with a Bioanalyzer

| Processing of illumina amplicon data
The HTS data quality was evaluated as following: in the initial step, excessive primer overhangs were clipped via cutadapt version 1.18 (Martin, 2011). Reads were then quality-filtered using the split.libraries.py command in QIIME version 1.8.0 (Caporaso et al., 2010) considering the following criteria: (1) containing exclusively unambiguous nucleotides, (2) containing the complete and correct forward or reverse primer (allowing four mismatches max.) and (3) having a minimum length of 90 bp after primer removal. In the final step, all reads underwent de novo chimera analysis in uchime version 5.2.236 (Edgar et al., 2011). All high-quality reads were eventually de-replicated into amplicons and clustered into OTUs (operational taxonomic units) in swarm version 2.2.2 (Mahé et al., 2015) using -d = 1 and the fastidious option -f (for further information on the use of Swarm, please see the developer's instructions at https://github. com/frede ric-mahe/swarm/ wiki/Fred%27s-metab arcod ing-pipeline). OTUs consisting of only one read (singletons) were removed from the Swarm clustering results. The remaining first-level clustering results of Swarm were subsequently used as input for a secondlevel clustering via sequence similarity networks at a 97% similarity threshold . The second-level clustering strategy aggregates highly similar first-level OTUs into one second-level OTU, thus preventing an over-splitting of genetic diversity. The OTU table was employed to a blast query against the genbank database (version 230.0) and against a custom database of ciliate sequences which we obtained from Sanger sequencing of species isolated from both lakes to assign taxonomies. This custom database comprised reference sequences of 15 different ciliates which were isolated from either Lake Mondsee or Lake Zurich. Isolation and identification of the ciliates as well as subsequent sequencing are described in Pitsch et al. (2019). More details on the custom ciliate reference database (e.g. GenBank accession numbers) are provided as supplementary information (SI 3).
However, not all OTUs could be assigned to a reference sequence for which an informative taxonomic classification was available. For downstream analyses, we recorded blast results of 'uncultured eukaryotes' as 'environmental hits', while considering every taxonomic information provided by the database for other 'uncultured' blast results (e.g. 'alveolate' in 'uncultured alveolate').

| Compositional variation analyses
All statistical analyses were conducted in r version 3.5.3 using the program packages vegan version 2.5-5 (Oksanen et al., 2019) and ggplot2 version 3.2.0 (Wickham, 2016) for graphical visualization. Rarefied (to smallest sample size) OTU richness was calculated for each sample as a measure of alpha diversity. The abundance-based OTU matrix comprising all samples was used as an input for nonmetric multidimensional scaling (NMDS; Kenkel & Orloci, 1986). In the initial step of the NMDS, the abundance-based data set was standardized by a Wisconsin double standardization and an additional square root transformation, both of which are recommended for improving NMDS results (Oksanen et al., 2019). After standardization, Bray-Curtis dissimilarity values across all pair of samples were calculated to generate distance matrixes for the NMDS analysis. Three environmental variables, water temperature, oxygen concentration and total chl a concentration, were fit to the ordination using the envfit function of the vegan package. The fitting (R 2 ) of each variable is based on its maximal correlation with the ordination configuration of the community data and detects the direction in the ordination space towards which the variable is most rapidly changing (Oksanen et al., 2019); the significance of each fit was assessed with 1000 permutations of the analysis.

| Co-occurrence network construction
There exist numerous strategies for network analyses of HTS data; however, not all of these are equally suited to detect co-occurrences between OTUs which are most likely to resemble interactions between organisms (Blanchet et al., 2020;Carr et al., 2019). One simple but important consideration is to infer co-occurrence networks not from presence-absence data, but from abundance data (Blanchet et al., 2020). Following this advice, the input data for our network analyses were abundance-based OTU-to-sample matrices. However, there are considerable difficulties owing to differences among correlation techniques which infer co-occurrences from abundance data (Weiss et al., 2016). Among the most severe problems is the detection of excessive false-positive signals, which are related to a large number of low-abundant OTUs and the high sparsity of microbial abundance matrices (Carr et al., 2019;Connor et al., 2017).
To overcome these issues and distinguish false positives from real signals, recent approaches have suggested to insert random noise into networks and to incorporate null models for identifying those signals which do not arise from neutral processes, but from complex biological associations (Connor et al., 2017). Although the use of null models requires permutations of the same network and is therefore computationally demanding, it is currently regarded as the most reliable strategy for inferring co-occurrence networks from microbial OTU abundance matrices (Carr et al., 2019). Taking these methodological considerations into account, we conducted our network analyses with NetworkNullHPC (https://github.com/lente ndu/Netwo rkNul lHPC), which is an implementation of the null model approach for inferring statistically significant co-occurrences in metabarcoding data sets (Connor et al., 2017). In this process, rare OTUs that occurred in <10% of all samples were removed from network analyses as they are unlikely to produce significant correlations. Read counts of remaining OTUs were normalized, and random noise was added to the OTU matrix for breaking ties in the Spearman rank correlation results. To derive the consensus network and conduct the null model approach on a statistically reliable number of network permutations, NetworkNullHPC conducts by default 1000 permutations of the network. The null model approach also avoids the use of arbitrary set thresholds for Spearman rank correlations applied across independent analyses, but uses the permutated networks to define a specific threshold for each independently calculated network. Only edges exceeding these Spearman correlation significance thresholds were maintained.
After obtaining the consensus network, we applied one more step of quality filtering which was inspired by recommendations in the review of Röttjers and Faust (2018), who proposed to agglomerate taxonomically highly similar nodes. This additional step offers the advantage to reduce excessive numbers of edges and nodes which may relate to the co-occurrence of OTUs from the phylogenetically same organism (autocorrelations). We implemented this workflow into our approach by merging nodes which shared a sequence similarity of at least 98% to the same reference sequence, while keeping all abundance and co-occurrence information of the individual nodes that were merged. In the same step, nodes which shared less than 98% sequence similarity to any reference sequence were removed from the networks.
The resulting edge lists of the four co-occurrence networks (Lake Mondsee epilimnion, Lake Mondsee hypolimnion, Lake Zurich epilimnion, Lake Zurich hypolimnion) are provided as supplementary information (SI 4 -SI 7).
Additional subnetworks were constructed for selected nodes of the four main co-occurrence networks. To create these subnetworks, we extracted the selected node along with its complete direct neighbourhood (i.e. all other nodes to which the selected node was directly connected by an edge) from one of the four main networks. The subnetworks allow for comparing the position and co-occurrence neighbours of specific nodes between different lakes as well as between different depths of the same lake.

| Statistical evaluation of network properties
Visualization of the networks was carried out in gephi version 0.9.2 (Bastian et al., 2009) using the Yifan Hu layout. Network metrics (average shortest path length, average degree, connectivity, diameter, modularity, transitivity) were assessed in R, using the package igraph version 1.2.2 (Csardi & Nepusz, 2006). A summary of all applied network metrics is shown in Table 1; further information about network metrics is briefly provided in the following.
We used betweenness centrality analyses to identify key nodes in the networks. Betweenness measures a node's capability of connecting two other nodes in the network (i.e. nodes with a high betweenness will be on a large number of shortest paths between pairs of nodes). Such three-way interactions are considered most promising when assessing the structure and function of microbial ecosystems (Connor et al., 2017). We inferred key nodes by applying a 95% confidence interval on the betweenness scores of all nodes within one network, as suggested by Röttjers and Faust (2018 TA B L E 1 Network metrics and their adaption to ecological questions a betweenness score exceeding the confidence interval were defined as key nodes of the respective network. The robustness of the network was assessed with the r package netswan version 0.1 (Lhomme, 2015). We simulated species extinction by performing a cascading attack scenario (Albert et al., 2000), in which the loss of connectivity was recorded when nodes were removed in decreasing order of their betweenness (as a proxy of their importance for the network structure), but with recalculating the betweenness after the removal of each node. A continuous decrease of connectivity will lead to the disintegration of a network, but this disintegration will be faster or slower depending on the specific network structure. If the associations that were mediated by the removed node can be maintained via other nodes, the network is robust towards the removal of that node. Thus, the loss of connectivity will be small and the network will not disintegrate into isolated subgraphs. Connectivity is used to determine the robustness of a network towards attacks and has been introduced into microbial ecology to assess the robustness of a community towards disturbances (Röttjers & Faust, 2018;Shi et al., 2016). Likewise, connectivity and robustness are measures which allow for assessing the complexity of a network (Albert et al., 2000).

| Compositional variation of protistan plankton communities
Ordination analyses of the Bray-Curtis (Figure 1) index resulted in two distinct clusters for Lake Mondsee and Lake Zurich samples, with no overlap of the confidence intervals for these clusters (stress: 0.075). These clusters were separated along axis 1, while axis 2 correlated significantly with temperature and oxygen con- netic assemblages in the shallower Lake Mondsee were not as distinct as in Lake Zurich.

| Co-occurrence network properties
Network analyses resulted in 160 nodes and 1122 edges for the epilimnetic (Figure 2a) and in 132 nodes and 843 edges for the hypolimnetic co-occurrence network of Lake Mondsee (Figure 2b), in 131 nodes and 522 edges for the epilimnetic (Figure 2c) and in 112 nodes and 224 edges for the hypolimnetic co-occurrence network of Lake Zurich (Table 2; Figure 2d). An overview on the proportions of nodes and edges shared among the networks is presented in Figure 3. Edges between pairs of OTUs that can be concordantly found in independent networks are likely to resemble actual biotic associations. The networks of Lake Mondsee shared 61% of their nodes and 22% of their edges. The networks of Lake Zurich shared 51% of their nodes and only 4% of their edges. Thus, even though a high proportion of nodes (OTUs) were shared between surface waters and deep-water layers, the associations (edges) of these OTUs to other OTUs were distinctively different in the different habitats.
A similar pattern emerged from habitat-specific comparisons across different lakes. LM hypo and LZ hypo shared only 28% of their nodes and 1% of their edges. These numbers illustrate that the portion of shared nodes between the networks was for all pairwise combinations of habitats and lakes considerably higher than the portion of shared OTUs between the respective pairs of input data sets. This observation can be retraced to low-abundant OTUs, which make up the majority of HTS data sets and often occur in single samples, but are rarely involved in statistically significant co-occurrences with other OTUs and therefore excluded from the networks. Their exclusion underlines the stringency of the network approach applied in this study, which rigorously removed co-occurrences with weak statistical support that could not be distinguished from background noise of the data set. Further network properties are detailed in Table 2 and can be summarized as follows: (i) networks of Lake Mondsee were more complex compared to Lake Zurich, (ii) networks of the epilimnion in both lakes were more complex compared to their hypolimnetic counterparts, and (iii) LM epi had the most complex and largest protistan plankton co-occurrence network, while LZ hypo showed the least complex one.

| Species extinction scenarios
With increasing removal of key nodes with high betweenness, the loss of connectivity (Figure 4)  TA B L E 2 Properties of protistan plankton networks obtained from the epilimnion (epi) and hypolimnion (hypo) of Lake Mondsee (LM) and Lake Zurich (LZ)

F I G U R E 3
Pairwise comparisons of shared and exclusive network nodes (OTUs) and edges (significant co-occurrences) between the cooccurrence networks constructed for protistan plankton of Lake Mondsee epilimnion (LM epi), Lake Mondsee hypolimnion (LM hypo), Lake Zurich epilimnion (LZ epi) and Lake Zurich hypolimnion (LZ hypo) (second most: Chrysophyceae, 5 nodes). In addition, Ciliophora also contributed the most nodes shared between all networks. These 5 nodes were taxonomically assigned to (ordered in descending read abundance): Vorticella sp., Tintinnidium fluviatile, Coleps sp. (reference sequence isolated from Lake Zurich), Urotricha sp. (reference sequence isolated from Lake Zurich) and Cyclidium marinum (probably a misleading entry in the GenBank database; we suppose a taxonomic placement among the morphologically hardly distinguishable Scuticociliata). The taxonomic affiliation of each node in every network is listed in the supplemental information (SI 9-SI 12) as well as a heatmap that visualizes the occurrence of each node affiliated to ciliates in the four networks (SI 13).
The importance of Alveolata nodes was further underlined by betweenness centrality analyses, though Stramenopiles contributed more key nodes to the Lake Mondsee hypolimnion network ( Figure 5). Unfortunately, the largest fraction of key nodes could not be clearly taxonomically assigned, which highlights the need for more morphological analyses even in ecosystems with long-term observational data. Among the Alveolata, the group which contributed most key nodes to each network were the Ciliophora (LM epi: 2 key nodes; LM hypo: 2 key nodes; LZ epi: 4 key nodes; LZ hypo: 1 key node). However, Ciliophora did not dominate the key nodes as clearly as may have been expected from the taxonomic composition of the networks. Bacillariophyceae (diatoms) emerged as key nodes in both networks of Lake Mondsee but not in Lake Zurich (LM epi: 5 key nodes; LM hypo: 8 key nodes; LZ epi: 0 key nodes; LZ hypo: 0 key nodes), while Chrysophyceae contributed similar amounts of key nodes to the networks as Ciliophora (LM epi: 1 key node; LM hypo: 4 key nodes; LZ epi: 2 key nodes; LZ hypo: 3 key nodes).

| Ciliate-specific subnetworks
An important and often neglected point in the context of network analyses is the attempt to link significant correlations between nodes in the networks to actual interactions between living species in nature. To contribute to this crucial point and to further analyse the similarities and differences among the two lakes and depth layers, we constructed subnetworks of four key ciliates as representatives: (i) one OTU, which was identical to the sequence of Urotricha castalia in our custom reference database ( Figure 6).
(ii) One OTU, which was 99.2% similar to the sequence of an undescribed ciliate in the NCBI nucleotide database (Accession no. KX465193, closest taxonomic blast hit was the haptorid Cultellothrix lionotiformis with 97.1% sequence similarity; Figure 7).
(iii) One OTU which was identical to the sequence of Paradileptus F I G U R E 4 Species extinction scenario conducted via a cascading attack. For each of the four co-occurrence networks (Lake Mondsee epilimnion: LM epi; Lake Mondsee hypolimnion: LM hypo; Lake Zurich epilimnion: LZ epi; Lake Zurich hypolimnion: LZ hypo), key nodes with the highest betweenness were stepwise removed while the loss of connectivity was recorded. This loss of connectivity documents the collapse of the network structure and is an indication of the robustness of each individual network elephantinus in our custom database ( Figure 8). (iv) One OTU which was identical to the sequence of Histiobalantium bodamicum in our custom database (Figure 9). Ciliates were chosen for these subnetworks, since decades of research in both lakes have led to profound knowledge about the behaviour and autecology of this group of organisms. Furthermore, ciliates also emerged as quantitatively (number of nodes) and qualitatively (key nodes with high betweenness centrality) important taxa within each network.
Three of the four focal ciliate species in the subnetworks were chosen because they were previously studied in at least one of the two lakes (Müller et al., 2002;Pitsch et al., 2019;Salbrechter & Arndt, 1994) and some could be cultured again for future experiments. The fourth species was chosen as an example for an unknown species, to highlight the potential of network analyses for facilitating species identification.
Urotricha castalia subnetworks: Urotricha castalia had notably more co-occurring taxa in LM epi, LM hypo and LZ epi, compared to LZ hypo ( Figure 6). However, surprisingly few connections were shared between the four U. castalia subnetworks. These included, Ciliate KX465193 subnetworks: The unassigned ciliate with the Accession no. KX465193 in GenBank had many co-occurring taxa in both Lake Mondsee water layers, even though it was characterized by a low sequence abundance (Figure 7). In LM epi (Figure 7a), KX465193 was predominantly associated with diatoms (nine species); one chrysophyte; followed by alveolates (six ciliates, two dinoflagellates and one apicomplexan); and green algae (six species). Whereas in the deeper water of Lake Mondsee, most of its connections were to other alveolates (four ciliates, one dinoflagellate and one undescribed alveolate); two telonemid flagellates; two green algae; one chrysophyte; one diatom; and the crustacean Mesocyclops. In Lake Zurich, KX465193 had a high relative sequence abundance, but its subnetwork comprised only very few connections: one in the epilimnion to the ciliate Cinetochilum margaritaceum; and five connections in the hypolimnion, three of which were to unknown eukaryotes, and two to other ciliates.
Paradileptus elephantinus subnetworks: The ciliate P. elephantinus only occurred in networks of the epi-and hypolimnion of Lake Mondsee. In the epilimnion (Figure 8a), this ciliate was associated with a similar set of taxa as the above-described Urotricha castalia in the corresponding network ( Figure 6). The subnetwork of P. elephantinus was notably less complex in the hypolimnion (Figure 8b) than in the epilimnion. In the deeper water of Lake Mondsee, P. elephantinus was only associated with four different stramenopiles (two Dinobryon species, one uncultured Chromulinaceae, one Mallomonas species).
Histiobalantium bodamicum subnetworks: The ciliate H. bodamicum only occurred in networks of the epi-and hypolimnion of Lake Zurich. In the epilimnion (Figure 9a), H. bodamicum was associated with three cryptophytes, two alveolates, one chlorophyte, one stramenopile, one fungus and one unknown eukaryote. In the F I G U R E 5 Pie charts displaying the taxonomic assignment of key nodes that exhibit significantly high betweenness centrality (using 95% confidence intervals) within protistan plankton networks of Lake Mondsee epilimnion (a), Lake Mondsee hypolimnion (b), Lake Zurich epilimnion (c) and Lake Zurich hypolimnion (d) hypolimnion (Figure 9b), the ciliate was associated with one unknown ciliate, one unknown eukaryote and to the choanoflagellate Salpingoeca fusiformis.

| Network analyses in microbial ecology
With the increasing amount of complex and massive data sets compiled by environmental sequencing studies, network analyses have become increasingly important in microbial ecology for interpretating these data (e.g. Karimi et al., 2017;Lima-Mendez et al., 2015).
Although many studies in different environments have impressively demonstrated how network analyses may be successfully applied to answer central ecological questions (e.g. Gilbert et al., 2012;Mikhailov et al., 2019;Shi et al., 2016), there still remain methodological limitations. Most importantly, all conclusions drawn from co-occurrence networks are based on the assumption that organisms whose abundances statistically significantly correlate across a given number of samples are likely to be ecologically related (Connor et al., 2017). However, the detection of a significant correlation alone does not inform about the nature of the relation between the organisms. There are several explanations for correlations among organisms. These range from shared preferences for the same environmental conditions to actual interactions between the organisms.
It is difficult to infer these reasons directly from networks. However, networks provide the basic information that allow to infer hypotheses across ecosystems (i.e. by comparing network structures among communities of different habitats) and also for individual species (i.e. by inferring autecological associations), which then can be further tested in experiments.
In the framework of this study, we conducted network analyses that implemented the most recent methodological improvements (Connor et al., 2017;Röttjers & Faust, 2018), which were even positively acknowledged by studies that urge caution in the interpretation of ecological network analyses (Blanchet et al., 2020;Carr et al., 2019). Our results (as well as previous studies, e.g., Lima-Mendez et al., 2015;Shi et al., 2016;Steele et al., 2011) highlight that network analyses may allow for novel insights into microbial F I G U R E 6 Subnetworks of Urotricha castalia extracted from the protistan plankton networks of Lake Mondsee epilimnion (a), Lake Mondsee hypolimnion (b), Lake Zurich epilimnion (c) and Lake Zurich hypolimnion (d). Node sizes mirror relative sequence abundances of each node (OTU) within the individual subnetwork. Taxonomic identities of nodes (OTUs) are provided as Supplementary Information SI 14 F I G U R E 7 Subnetworks of the uncultured ciliate KX465193 extracted from the protistan plankton networks of Lake Mondsee epilimnion (a), Lake Mondsee hypolimnion (b), Lake Zurich epilimnion (c) and Lake Zurich hypolimnion (d). Node sizes mirror relative sequence abundances of each node (OTU) within the individual subnetwork. Taxonomic identities of nodes (OTUs) are provided as Supplementary

| Water turnover dynamics affect protistan plankton communities
Previous studies have shown that ongoing global warming has distinctively affected the environmental conditions in Lake Zurich (North et al., 2014;Posch et al., 2012;Yankova et al., 2017). This trend climaxed in an abrupt stop of deep-water turnover in the lake, which turned from a holomictic (with complete turnover) into a periodically meromictic ecosystem (with only upper parts of the water body affected by turnover, see Yankova et al., 2017). This drastic environmental change led to consequences for the protistan plankton community, such as the absence of traditional phytoplankton spring blooms predicted by the PEG model . To cope with these changes, organisms have to adapt In addition, the above-described different mixis patterns of the two lakes may contribute to a higher dissimilarity of epi-and hypolimnetic plankton communities in Lake Zurich (maximal depth: 136 m) compared to Lake Mondsee (maximal depth: 68 m). In monomictic Lake Zurich, the deep hypolimnion gets infiltrated by surface water only once a year in spring (or not at all, when thermal stratification impedes mixis), whereas water turnover affects deep layers in Lake Mondsee twice a year. Thus, the different frequency of water exchange may lead to a higher dissimilarity of microbial communities in Lake Zurich.

| Plankton richness, ecosystem robustness and functional diversity
The differences in protistan plankton composition are further reflected by the co-occurrence networks which were constructed from these data sets. It may be expected that differently composed communities result in differently composed co-occurrence networks. But it is intriguing that despite similar levels of OTU richness in the input data sets, the resulting network structures F I G U R E 9 Subnetworks of Histiobalantium bodamicum extracted from the protistan plankton networks of the epi-(a) and hypolimnion (b) of Lake Zurich. Node sizes mirror relative sequence abundances of each node (OTU) within the individual subnetwork. Taxonomic identities of nodes (OTUs) are provided as Supplementary Information SI 17 were considerably different between the two lakes. All network metrics pointed towards a less tightly connected and more fragmented network structure for Lake Zurich communities in which associations between larger groups of organisms are impeded.
By contrast, associations between larger groups of co-occurring organisms are easier established in networks of Lake Mondsee protistan plankton communities. This finding underlines the potential of network analyses for uncovering differences between communities that go beyond conclusions which can be drawn from diversity analyses.
On top of that, in a scenario which modelled the successive extinction of key taxa after perturbation, both protistan plankton networks of Lake Zurich emerged as notably less robust and thus, more vulnerable than the Lake Mondsee networks (Figure 4). Several studies have shown that network metrics such as connectivity provide valuable insights into the structure of microbial communities and may predict responses of the communities towards disturbances (see, e.g. Dunne et al., 2002;Proulx et al., 2005, and references therein; Röttjers & Faust, 2018;Shi et al., 2016). The tenor of these studies is that a low robustness modelled from extinction scenarios of the community network can be interpreted as a low robustness of the studied ecosystem.
The complex consequences caused by a loss of biodiversity in general and a decline of species richness in specific are discussed among others as major reasons for lower ecosystem stability (Tilman et al., 2006). But despite the general notion that climate stress results in a decline of multicellular and microbial species richness (e.g. Banerjee et al., 2020;Bellard et al., 2012), the average protistan plankton (OTU) richness was highly similar in both lakes under study. Thus, we agree with previous studies that a high species richness alone is insufficient for making predictions on ecosystem stability (May, 1973;Tilman, 1996). Network analyses may be more suited for this task, since they integrate species richness but focus on potential interactions between the species for assessing the complexity of the studied ecosystem. Our specific strategy of performing targeted attacks on the community (removing key nodes with high betweenness centrality), and inferring the robustness of the community networks from the loss of connectivity corroborates earlier work by Dunne et al., 2002. Similar to our own approach, the authors concluded that connectivity may serve as an indicator of ecosystem complexity and robustness. Furthermore, removing a species from a network by a targeted attack often leads to secondary extinctions of other species (because these were only connected to the attacked species) or the disintegration of the network (because the attacked species connected two otherwise independent modules of the network).
The total effect of the species extinction thus indicates the species' functional position within the ecosystem's network (Dunne et al., 2002). As reviewed by Weisse (2017), there is strong evidence that communities with greater functional diversity are more stable, while greater ecosystem stability implies lower susceptibility to perturbation. Taking the results of the species extinction scenarios into account, we therefore assume the lower robustness of the co-occurrence networks towards disturbances may indicate a lower functional diversity in protistan plankton communities in Lake Zurich than in Lake Mondsee.

| Predictions by network analyses may complement our perception of species autecology
As expected from previous microscopical studies (e.g. Lischke et al., 2016;Posch et al., 2015;Sherr & Sherr, 1988;Weisse, 2006;Weisse & Sonntag, 2016), ciliates emerged as major contributors to the co-occurrence networks of both lakes, potentially acting as predators of bacteria, algae, flagellates, other ciliates and even some metazoans. In addition, ciliates represent an important food source for zooplankton, thus, ultimately, supporting the productivity of lakes (e.g. fish stock). Mixotrophic ciliates and functionally autotrophic species may even significantly contribute to primary production in lakes.
Several ciliates were found in networks of only one of the two studied lakes. Such differences may reflect the preference of species with specific traits best adapted to local environmental conditions . Examples for ciliates which occurred in only one of the two lakes are Paradileptus elephantinus (exclusively detected in Lake Mondsee networks) and Histiobalantium bodamicum (exclusively detected in Lake Zurich networks). Paradileptus elephantinus is a conspicuous dileptid ciliate which is widespread but usually not abundant in the plankton of lakes and large rivers (for a compilation of autecological data, see Vďačný & Foissner, 2012).
Due to its large size of up to 600 µm, we consider this ciliate as a planktonic flagship species that is a voracious predator hunting on bacteria, protists and rotifers. Thus, by interacting with several other organisms and having a strong top-down effect on diversity, P. elephantinus meets one definition of a keystone species (Davic, 2003).
However, the associations with numerous diatoms and chlorophytes strongly link the occurrence of this species to the availability of such algae as possible prey. In Lake Zurich, diatoms and chlorophytes are largely outcompeted: total phytoplankton is in general dominated by the cyanobacterium Planktothrix rubescens, which accounts for ~80% of total biomass. This situation decreases the importance of ciliates such as P. elephantinus and selects for species with traits better adapted to local conditions in Lake Zurich. Such a species is Histiobalantium bodamicum, which rejects diatoms and prefers small flagellated cryptophytes as prey (Müller & Schlegel, 1999), to which it was also connected in its subnetworks (Figure 9). This concurs with abundance peaks of H. bodamicum in Lake Zurich during spring with up to 35 Ind. ml −1 accounting for 55% of the total ciliate morphotype abundance, when also cryptophytes were highly abundant (Posch et al., 2015).
In spite of different lake effects, several ciliate species co-occurred in Lake Mondsee and Lake Zurich and contributed considerably to the protistan plankton networks of both lakes. At first thought, this seems to contradict the selection of species according to their response traits to environmental factors as discussed above.
However, these ciliate species are characterized by a high versatility in functional traits and, thus, a high capability of adaption to different ecological conditions and to different sets of locally occurring interaction partners. This conclusion can be drawn from the distinct subnetworks of the species in the two lakes. Representative examples discussed here in detail are Urotricha castalia and an unknown ciliate. Urotricha castalia is a common ciliate in lake plankton causing the breakdown of classic eukaryotic vernal phytoplankton spring blooms, together with oligotrich ciliates and small metazooplankton such as rotifers (Sommer et al., 1986). In temperate lakes, algal mass developments during spring are in many cases formed by small-sized, fast-growing centric diatoms (stramenopiles) and several genera within the cryptophytes. These algae seem to be not only in the adequate food size range of various typical spring bloom ciliates, but are probably also of high nutritional value, as primary producers have low N:P nutrient requirements, thus being a phosphorus source for ciliates. Experimental studies showed that U. castalia also grew exceptionally well on small-sized cryptophytes as a food source (Weisse et al., 2001). Its versatile foraging trait allows for different responses of U. castalia, which are triggered by the prevailing resources in both lakes: diatoms and green algae in Lake Mondsee and cryptophytes in Lake Zurich ( Figure 6). Its low sequence abundance suggests that U.
castalia may be top-down controlled by predation. This is supported by the ciliate's associations to several metazoans in the subnetworks.
The unknown ciliate, which has a key position in the co-occurrence networks of both lakes, has a high sequence similarity to an undescribed ciliate with GenBank Accession no. KX465193 (Luo et al., 2017) detected in Lake Oloidien in Kenya. Previous studies on ciliates in Lake Zurich and Lake Mondsee confirmed the occurrence of several small-sized haptorid species of the genera Askenasia, Mesodinium or Didinium to which the unknown species may belong (e.g. Pitsch et al., 2019;Salbrechter & Arndt, 1994). Its association with ciliates of the bacterivorous genera Cinetochilum and Cyclidium and its significant connection to autotrophic and heterotrophic flagellates indicate that this unknown ciliate may also have a versatile lifestyle feeding on different types of small-sized prey organisms. Its high (sequence) abundance in spite of its assumed small size speaks for a low sensitivity of this species towards predation. Haptorids, in general, are well-adapted to predation pressure by performing 'escape jumps' or the use of extrusomes as a defence mechanism (for details see Foissner et al., 1999). This example as well as other un- identified key nodes in all four habitat-specific networks (Figure 2) demonstrates the necessity to intensify analyses of morphological protistan inventories in lakes. In addition, these results demonstrate the strength of network analyses based on molecular data, which allow the identification of important but undescribed key taxa in the pelagic food web and thus guide the design of further research studies with the aim to morphologically identify such taxa and investigate their autecology in experimental studies.
Our results demonstrate that it is crucial to identify the full set of potential functional traits of individual taxa under different environmental conditions if we want to understand community function at an ecosystem level. The general practical framework for functional diversity analysis in an ecosystem context is based on a functionaltrait matrix and a relative abundance (or biomass) matrix of each species in a community (Villéger et al., 2008). Functional diversity can then be linked to matrices of environmental variables and ecosystem properties such as nutrient cycling, productivity, robustness or sensitivity to invasions. The functional traits of several ciliates (and other protists) might be much more complex than previously thought and might vary as a response to environmental changes (see also discussion in Weisse, 2017). While traits of microbial organisms are best analysed by microscopical (and possibly molecular) approaches, the design of such analyses will profit from associations predicted in networks, which appear as a powerful tool to phrase hypotheses about individual functional trait versatility under different environmental settings.
Once confirmed, the expressed functional traits may serve as indicators for present ecosystem status and expected changes in the future.

DATA AVA I L A B I L I T Y S TAT E M E N T
Illumina amplicon data NCBI SRA: SRR11351406-SRR11351499.