Changes in Ciliate Communities Reveal Modi cation of Lake Functioning Over the Last Century

Barouillet Cécilia (  cecilia.barouillet@gmail.com ) INRAE, Université Savoie Mont Blanc, CARRTEL Valentin Vasselon INRAE, Université Savoie Mont Blanc, CARRTEL François Keck INRAE, Université Savoie Mont Blanc, CARRTEL Laurent Millet CNRS, Chrono Environnement David Etienne Université Savoie Mont Blanc, INRAE, CARRTEL Didier Galop GEODE UMR 5602 CNRS, Université de Toulouse Damien Rius CNRS, Chrono Environnement Isabelle Domaizon INRAE, Université Savoie Mont Blanc, CARRTEL


Introduction
Lakes are largely recognized as integrators and sentinels of environmental changes 1 . Pressures from anthropogenic activities have largely increased in magnitude since the mid-twentieth century, a period that has also been referred to as the Great Acceleration 2,3 . More speci cally, climate change and the human alteration of the landscapes can have a profound impact on the physical and chemical characteristics of lakes 4,5 , thereby in uencing the communities living in and depending on these ecosystems 6 . While challenging, assessing the biological response to environmental changes over large geographical scales can provide important insight into the vulnerability of lakes to anthropogenic and climate forcing. The top-bottom paleolimnological approach allows such assessment for aquatic communities through the comparison of sedimentary archives of recent (i.e. top samples) to past communities (i.e. bottom sample) 7 . This comparative approach is quite powerful as it can be applied on large geographical scales and provides important insight into reference conditions (i.e. prior to major anthropogenic in uences), thereby allowing an assessment of the amplitude of change 8,9 . Although the diagenesis and archiving of DNA in the sediments through time is often discussed in literature 10 , previous studies have demonstrated that the genetic information of the micro-organisms living in the water column archives in the sediments [11][12][13] . The development of molecular biology techniques to target DNA molecules preserved in lake sediments (sed-DNA) largely expanded the eld of paleolimnology over the last decade 14 . For instance, applying molecular biology tools in paleolimnological investigations allowed to reconstruct historical freshwater biodiversity 15 , as well as uncovered past ecological changes of a wide range of organisms, including specimens for which visible remains are not preserved 16, 17 . This includes the reconstruction of the long-term dynamic of micro-organisms 18 which represent the largest source of biodiversity and ecological functions, but have long been omitted from debates concerning the response of the global biodiversity, in part due to the lack of data 19 . As such, sed-DNA offers new opportunities to reconstruct a more holistic view of the long-term response of lakes to environmental changes 18,20 . For instance, through the implementation of new paleoindicators, sed-DNA can allow to reconstruct the long-term dynamic of overlooked communities such as the communities of the microbial loop that play a key role in the functional ecology of lakes.
Ciliates are unicellular microorganisms occupying diverse ecological niches and a wide biogeographic range 21 . These protists display a large functional diversity, acting as predators of bacteria and other protists, including algae, as well as small metazoans 22 , while mixotrophic ciliates can signi cantly enhance primary production 23,24 . Altogether, they play a key role in aquatic pelagic and benthic food webs, especially in the transfer of energy and nutrient cycling from the microbial loop to the metazooplankton 21,25 . Long-term reconstruction of the microbial eukaryote diversity of large peri-alpine lakes reported a strong relative contribution of the ciliate communities to the total abundance of microeukaryotic communities and demonstrated that ciliates are sensitive to changes in phosphorus concentration 13 . Although the ciliates are good indicators of environmental changes 25,26 , they are seldom used in routine limnological surveys. Indeed, the analysis of ciliate communities through microscopicbased approach and morphospecies identi cations is challenging 27 , mostly because the monophyletic group of ciliates includes very different organisms that display diverse and complex life cycles 28 . Recent advances in molecular biology and the development of more robust taxonomic libraries continue to strengthen eDNA and sed-DNA techniques 28,29 , thereby allowing to include these understudied groups and useful indicators for a more holistic ecological diagnosis.
In the present study, we combine metabarcoding methods and the top-bottom paleolimnological approach to reconstruct recent and past ciliate communities of 48 temperate lakes ( Fig. 1), for which drastic changes in the composition of the micro-eukaryotic communities were previously observed as a response to the Great Acceleration 18 , with a signi cant increase in phototrophic and mixotrophic communities consistent with the global enhancement of primary productivity. Focusing on these same lakes, we aimed at exploring how the heterotrophic and mixotrophic group of ciliates responded to these recent changes by (1) assessing the amplitude of changes in the ciliate communities across a wide variety of lakes and (2) evaluating the potential of using the ciliate as indicators of functional and biological changes.

Results
Ciliates metabarcoding analyses. The sequencing resulted in a total number of 2,746,319 DNA reads with an average of 28,650 reads per sample. After the ltering steps, 1,745,532 reads were retained and clustered into 2446 OTUs. Detailed information about the effect of bioinformatics treatments on the DNA reads and a summary of the resulting total number of OTU and number of reads taxonomically assigned and associated to a functional trait are accessible in the supplementary information (Supplemental Tables S1, S2).
Ciliates community diversity. The NMDS based on the Bray-Curtis distance showed that the overall dispersion of the recent and past samples differed, with the distance between the samples and their geometric median being much more variable and greater for the bottom samples (Fig. 2).
The average distance between the samples and their centroid was signi cantly different between the past 0.79 and recent samples 0.35 (Wilcoxon test, p<0.01). However, there was no large displacement of the centroid relative to each other as illustrated by the overlapping con dence ellipses of the two groups (Fig.  2a). This pattern was even more pronounced on the NMDS calculated from the OTU table (Supplemental Fig. S1). The hierarchical analysis identi ed four clusters according to their community composition, separating most past samples from the bottom samples (Fig. 3).
More speci cally, the clusters 1, 2 and 3 were only composed of past samples except for the cluster 2 which contained the present and past samples of two lakes that did not differ. The samples of cluster 1 were dominated by the Oligohymenophorea (Class) Scuticociliatia (Subclass), samples of cluster 2 by Colpodea (mainly Cyrtolophosidida), while samples of cluster 3 were dominated by Spirotrichea Hypotrichia (Fig. 3). The fourth cluster regrouped the samples for which the ciliate communities were dominated by the Class Litostomatea, subclass Haptoria. Although the fourth cluster regrouped past and recent samples, sub-clusters indicated a clear separation of the recent from the past samples. The SIMPER analysis applied on the subclass level indicated that the increase of the Haptoria, the Oligotrichia and Armophorea contributed by 26%, 11% and 2% respectively to the Bray-Curtis dissimilarity between the recent and past samples, while the decline of Colpodea, Hypotrichia, Scuticociliatia and Prostomatea contributed by 19%, 18%, 6% and 3% to the dissimilarity. The differential abundance analysis (DESeq2) provided information about the intensity of the changes between the recent and past samples 30 and applied at the Genus level indicated a signi cant increase of ve genus and families of the subclass Haptoria, two genus and families of the subclass Oligotrichia and the genus Metopus of the class Armophorea, as well as a signi cant increase of the Peritrichia and Suctoria (Supplemental Fig. S2). The DESeq2 analysis also indicated a signi cant decline of Urotricha a genus of the class Prostomatea, two Scuticociliatia including Cyclidium, ve genus and families associated to the class Colpodea, and ve genus and families associated to the class Hypotrichia, as well as the signi cant decline of Peniculia. The regression tree analysis identi ed a signi cant split in the dataset at an elevation of 1400m, with a Bray-Curtis dissimilarity coe cient averaging 0.37 for lakes situated above 1400m in elevation while lakes situated under 1400m in elevation displayed an average of 0.57 (Fig. 4).
Ciliates functional traits and lakes ecosystem functioning. The DESeq2 indicated a signi cant increase of the facultative or strict anaerobe benthic ciliates while the benthic ciliates that require oxygen signi cantly declined in the recent period (Fig. 5). This increase in facultative or anaerobe benthic ciliates was mostly in uenced by the signi cant increase in the abundance of the genus Metopus (order: Armophorida, class: Armophorea) (Supplemental Fig. S2). Moreover, there was a signi cant increase in sessile ciliates (i.e. attached to a substrate), as well as a signi cant increase in the pelagic ciliates that tends to be found in the epilimnion or distributed in the entire water column during mixing periods 31,32 . In contrast, ciliates associated with the hypolimnion of the pelagic zone signi cantly declined. In addition, the same analysis applied on the foraging trait categories indicated a signi cant increase in mixotrophs and non-selective lters feeders (which includes some detritivores), while bacterivores declined signi cantly in recent samples (Fig. 5).

Discussion
Over the last decades, lakes have been exposed to environmental changes (anthropogenic stressors and natural changes) with important implications on the biological communities. Recent paleolimnlogical studies have provided new insights on the long-term responses of microbial assemblages by using DNAbased methods 15 . Although previous paleolimnological investigations of overall microeukaryotes communities reported a strong response of the ciliates to nutrients inputs 13,33 , speci c long-term changes in the ciliates communities structure and functional ecology had not yet been investigated. The present study represents the rst paleolimnological reconstruction of ciliate communities and demonstrate the potential of using heterotrophic protists as indicators of environmental change at the community and functional levels.
Analyses of the ciliate communities indicate an overall decline in the β-diversity in recent times following the same trends as the overall micro-eukaryotes diversity 18 . Through the use of primers speci cally targeting ciliates, our study brings an innovative perspective and provides additional information by highlighting changes in some heterotrophic and mixotrophic communities that were not previously revealed in the analysis of the overall microeukaryote communities using generalist primers. Interestingly, the ciliates did not undergo a large turnover but rather our results indicate a spatial homogenization of the diversity with a reorganisation of the community structure (i.e. switch in dominance) (Figs. 2 and 3). This spatial homogenization of the ciliate communities was marked by the replacement of three clusters by one homogeneous community dominated by Haptoria across all lakes. Biotic homogenization of communities is a well documented phenomenon that has been observed in terrestrial and aquatic ecosystems largely in uenced by a reduction of environmental heterogeneity and the availability of diverse ecological niches [34][35][36] . In aquatic ecosystems, change in climate and productivity, and anthropogenic alteration of the watershed are the most prevalent causes of biotic homogenization 37 . Several of our studied lakes are exposed to similar stressors which includes nutrient enrichment, agricultural and urban development of watershed and climate change 38-41 . As such, our results support that these factors acted as deterministic lters selecting for a more homogeneous group of species dominated by generalist ciliates that display more exible life strategies.
A stronger restructuration of the ciliate communities was observed in the low elevation lakes demonstrating that environmental changes in lowland lakes impacted several trophic levels, including non-photosynthetic protist communities. The geographical variation in the amplitude of changes in diversity and community turnover of microorganisms associated with an elevation gradient have been previously demonstrated in terrestrial 42 and aquatic ecosystems 43 . These patterns are explained by the more pronounced human footprint in lowlands 42,44 , including nutrient-enrichment in freshwater ecosystems 6 . Supporting this trend, the present day trophic status of our studied lakes was signi cantly higher for lowland lakes than for high elevation lakes (Supplemental Fig. S3). As such, human-induced nutrient increase might have in uenced the observed changes in the ciliate community diversity of lowland lakes.
Changes in the dominant functional groups of ciliates in the 48 studied lakes highlighted modi cations of the physico-chemical conditions and biotic interactions of the pelagic and benthic zones.
The recent increase in the ciliates that can display mixotrophic life strategies whereby ciliates harbor algal endosymbionts or sequester plastid from their prey 45 , suggests that mixotrophic ciliates are directly or indirectly responding to new environmental conditions 45 . Interestingly, mixotrophic organisms tends to bene t from a more exible nutrition strategy and have been found to prone during transition phases between autotrophy-dominated and heterotrophy-dominated systems 46,47 . As such, they can easily adapt to the more frequent exposure to extreme events or highly-variable environmental conditions that have occurred in lakes over the last decades 48 . Moreover, Limnostrombidium (synonym of Strombidium, from marine waters) is a common freshwater genus regrouping specialists mixotrophic species that preferentially feed and use the chloroplast of picophytoplankton 31,45,49 . Their signi cant increase in recent time thus indicates that recent changes in autotrophic picoplankton dynamic and structure, that have been previously recorded in some of our studied lakes 16 , might have provided them with a competitive advantage. Importantly, however, more empirical studies are still needed in order to solidify inference made between the relative increase in mixotrophic life strategies and associated ecological conditions. Mixotrophic ciliates play a major role in foodweb structure of freshwater lakes 24,50 . The mixotrophic ciliate can contribute to the enhancement of the primary production 51 , while as prey, they represent a more direct transfer of the solar energy to the zooplankton 24 . Through this process, mixotrophic ciliates can enhance the e ciency of carbon transfer and energy ow in the food web 24 . Their recent increase thus suggests that lakes might have undergone important trophodynamic changes as mixotrophy is becoming an increasingly important pathway in aquatic food webs. These results also highlight the importance of integrating the mixotrophic component in food web modelling and carbon ow studies 52 .
Change in the benthic ciliate communities indicates that the benthic environment has also been impacted by recent environmental changes. The signi cant increase in the facultative or obligate anaerobic ciliate associated with the benthos, such as Metopus, suggests that the ciliate communities have been directly in uenced by the widespread deoxygenation of temperate lakes 5 (Fig. 5 and Supplemental Fig. S2). The signi cant decline in the benthic and hypolimnitic ciliates associated with well-oxygenated conditions further support that the habitability of the sediment-water interface has been declining for this particular group of ciliates. The depletion of oxygen concentrations in the profundal zone of freshwater lakes is a well recorded global phenomenon that can have a pervasive impact on the overall ecosystem functioning 53 . These changes have been associated with stronger and longer thermal strati cation, as well as a loss of water clarity, in part due to the increases in pelagic production 5 . Supporting this hypothesis, some of our studied lakes, where the strict anaerobe bacteriophage Metopus have been recently thriving, have experienced unprecedented episodes of eutrophication or cyanobacterial bloom over the last decades and subsequent periods of deep water hypoxia [54][55][56] .
Signi cant changes in several other functional groups of ciliate provide additional evidence of recent changes in aquatic food web structures and habitats. For instance, the signi cant increase in pelagic ciliates that preferentially lives in the strati ed epilimniun further support that lakes have been exposed to longer and stronger periods of strati cation in recent time. Furthermore, periphytic species, such as the sessile or sedentary forms Peritrichia and Suctoria, is consistent with more frequent pelagic blooms or macrophytes under higher nutrient and warmer conditions 57,58 . However, the lack of knowledge about the ciliate ecology, their biotic interactions and regulatory factors limit our capacity to decipher the relative importance of each potential stressors and the underlying mechanisms remain elusive.
Altogether, the diagnosis of the changes in the ciliate communities across the 48 studied lakes supports the use of ciliates as indicators of environmental changes 59 and provides evidence that the ciliate communities are strongly responding to the environmental changes that have occurred over the last century which includes widespread deoxygenation of deep waters, changes in thermal strati cation and nutrient-enrichment. Playing a key role in the metabolic pathways of aquatic ecosystems 51,60 , they can provide valuable insight into the functional ecology of lakes, and their strong response recorded in the sedimentary archives suggests important changes in the main pathways for the transfer of energy within the microbial food webs 61 . Although, molecular studies of protist communities and working with ancient DNA is quite challenging with several aspects that need to be taken into consideration for future studies (as summarized in Methods section of the Supplemental Material), the present study as well as several previous studies 27,51,62 have shown the great success of using such approaches to assess environmental changes in aquatic ecosystems. Their integration to environmental assessment using high-throughput sequencing and metabarcoding technologies is thus promising, providing a more holistic overview of the response of aquatic ecosystems to environmental changes. This is even more relevant as the science is moving toward ecosystem-wide food web modelling 20,51,63 , and protists, as key players of the microbial food-web, serve an important function of recycling carbon and energy in lakes.

Methods
Study sites, sediment core collection and subsampling A total of 48 lakes were used in this study (Fig. 1, Supplemental Table S3). The studied lakes were selected to be located along a large elevation gradient and displayed various physico-chemical characteristics (Supplemental Table S3). For all 48 lakes, the sediment cores were retrieved between 2010 and 2016 (cf. Keck et al. 18 for more detailed information on each lake) from the deepest part of the basin using a UWITEC gravity corer. Prior to subsampling, the core was air protected by a double layer of plastic wrap and stored in the dark at 4°C.
A top-bottom technique was used to provide a simpli ed assessment of the amplitude of change in the ciliate community diversity and composition. Brie y, a top sample corresponding to recent time and a bottom sample corresponding to the past were isolated from the cores. Importantly, the top samples (i.e. the last decade) were selected to be a few centimeters down the core to avoid bias associated with active benthic ciliate taxa and early diagenesis processes 12 . The depth of the bottom samples was chosen for each core in order to correspond to the pre-"Great Acceleration" period 2,3 (i.e. nineteenth century) which was determined using a combination of several approaches (i.e. X-ray uorescence, radiocarbon and radionuclides 210 Pb and 137 Cs; cf. the supplementary material and methods in Keck et al. 18 for more detailed information on each lake). In order to ensure that each sample covered at least 10-15 years, the thickness of the sediment samples was individually adjusted. Sub-sampling for DNA analysis was conducted in a clean and controlled environment using strict laboratory protocols to avoid contamination by modern DNA 15,18 .

Molecular Analysis
Two DNA extractions were performed on 0.5 g of wet sediment for each sample using the NucleoSpin® soil kit, according to the manufacturer instructions (Macherey-Nagel, Düren, Germany). The same DNA extracts as in Keck et al. 18 were used. Refer to the Method section of the Supplementary Material for more details regarding the laboratory protocol and conditions applied for the DNA extraction. A nested-PCR targeting the V7 region of the 18S rRNA gene was used to do the inventory of the ciliate community. In the rst step, a set of primers was used to target a speci c DNA region for ciliates of 800bp CS322F (5'-GATGGTAGTGTATTGGAC-3') and 1147R (5'-GACGGTATCTRATCGTCTTT-3') 64,65 . The rst PCR was performed in a total volume of 25 µL containing 1 µL of DNA extract, 2.5 µL of 10X NH4 reaction buffer, 2 µL of 50mM of MgCl2, 0.5 µL of 100mM dNTP, 1.25 µL of each primer at 10 pmol/µL, 2 µL of 10mg/mL BSA and 0.1 µL of 5Ci BioTaq (Bioline). The ampli cation cycle included an initial denaturation at 95°C for 10 minutes followed by 20 cycles of 15 seconds at 94°C, 15 seconds at 57°C and 30 seconds at 72°C. The amplicons were then subjected to a nal 10 minutes extension at 72°C. The second PCR was then applied on the products of the rst PCR using the general eukaryotic primers 960 F (5′-GGCTTAATTTGACTCAACRCG-3′) and NSR1438 (5′-GGGCATCACAGACCTGTTAT-3′), previously used by Keck et al. 18 , amplifying in DNA fragment of about 250bp. Molecular tails were added to the forward primer (5′-CTTTCCCTACACGACGCTCTTCCGATCT-3′) and to the the reverse primer (5′-GGAGTTCAGACGTGTGCTCTTCCGATCT-3′). The second PCR was performed in a total volume of 25 µL containing 0.8 µL of DNA from the rst PCR, 2 µL of 10X NH 4 reaction buffer, 1.6 µL of 50mM of MgCl 2 , 0.4 µL of 100mM dNTP, 1 µL of each primer with molecular tails at 10 pmol/µL, 1.6 µL of 10mg/mL BSA and 0.06 µL of 5Ci BioTaq (Bioline). The ampli cation cycle included an initial denaturation at 95°C for 2 minutes followed by 20 cycles of 30 seconds at 94°C, 30 seconds at 57°C and 45 seconds at 72°C. The amplicons were then subjected to a nal 10 min extension at 72°C. The nested-PCR protocol was applied on each DNA extractions separately, the full volume of the nal products resulting from the two DNA extracts of the same sample were then pooled and sent to GeT-PlaGe (Plateforme Génomique 31326 CASTANET-TOLOSAN Cedex) for amplicon puri cation, library preparation and paired-end (2 × 250 bp) sequencing on a MiSeq Illumina instrument (San Diego, CA, USA).
The reads were demultiplexed and R1/R2 reads assembled into contigs by the sequencing platform who provided one fastq les per sample. The high-throughput sequencing data were then cleaned in Mothur 1.45.1 66 . Filtering steps were used to conserve DNA sequences of 350±50 bp in length, with no ambiguities (N=0), 10 or less homopolymer (max homopolymer=10) and no mismatch was allowed in the primer sequence. The data was dereplicated in order to work with Individual Sequence Unit (ISU). ISUs were then aligned using an aligned version of the Silva 18S database restrained to the V7 region and ISUs that were not fully aligned to the Silva 18S barcode were removed. The detection and removal of chimera was done using Vsearch as implemented in Mothur with default parameters. The taxonomic assignment of the ISU was done using a curated version of the Protist Ribosomal Reference database PR2 67 "pr2_version-4.12.0_18S_cil_cur" (available on Zenodo repository system: https://doi.org/10.5281/zenodo.5163167) and using the command classify.seqs() and the method wang with a con dence score threshold of 75% and 100 iterations. Following this rst taxonomic assignment, the ISU represented with only one read or that were identi ed as "unknown" or "Eukaryota_unclassi ed" were removed. The ISUs were then clustered into molecular Operational Taxonomic Unit (OTU) using the furthest neighbour algorithm with a similarity threshold of 97%. Finally, the command classify.otu() was used to taxonomically assign the OTUs based on the rst taxonomic assignment of the ISUs with a con dence threshold of 80%. The OTUs that did not belong to the Phylum Ciliophora were removed. The taxonomic a liations were checked and harmonized manually using the classi cation from Gao et al. 68 . In order to study changes in the functional groups, the ecological preferences (i.e. preferred limnetic habitat or foraging traits; cf. Table 1) were indexed for the OTUs for which the taxonomic a liation was ne enough (at least assigned to the family rank); otherwise, the category "Unknown" was given. The association of OTUs to their functional traits was done through an exhaustive literature review, the foraging traits categories created were inspired from a combination of several previously published categories based on the feeding ecology of ciliates 69, 59,70,22 .  Ciliates that can in found in the littoral zone or well-oxygenated bottoms.

Commensals or Parasites
Freshwater ciliates that are either endocommensals of bivalves or ectocommensals of sh, or parasites of sh.

Facultative or obligate anaerobe
Facultative or obligate anaerobe living in the benthic environment but also includes some taxa able to live in anaerobic deep waters.

Pelagic
Ciliates that are found in high abundance in the epilimnion during the strati cation period and can be found in low abundance throughout the water column during spring and fall mixing.

Pelagic (hypolimnion)
Ciliates that preferentially lives in the hypolimnion or metalimnion during the strati cation periods and are generally found in low abundance during mixing periods. Ciliates from this group do not withstand anoxic conditions.

Sessile
Ciliates for which their life cycle includes a stage attached to a substrate. Usually stalked ciliates.
Categories associated with foraging strategies Algivores Herbivorous ciliates

Bacterivores
Ciliates that exclusively feed on bacteria, these ciliates are usually associated with the benthic environment or found in the metalimnion of highly productive lakes

Carnivorous
Ciliates that feeds on other ciliates or even small metazoans

Commensals and parasites
Regroup parasites, bacterivores and histophage ciliates. Commensals and parasitic ciliates were kept separated as they are more likely to be directly in uenced by the presence/absence of their hosts rather than in uenced by changes in the biotic and abiotic factors of the surrounding environment.
Filter-feeders (nonselectives) Passive lter-feeders and detritivores. Ciliates from this category are associated with highly productive aquatic environment or environment with a high concentration of dissolved organic matter (i.e. often found in polluted waters or waste water treatment plants).

Fungivorous
This categories contains only one species found in our samples Pseudoplatyophrya nana

Mixotrophs
Phagotrophic ciliates that harbor algal endosymbionts or sequester plastid from their algal prey. Mixotrophic ciliates thus tends to be also algivores.

Predators of small protists
Omnivorous ciliates feeding on algae and other small ciliates

Statistical analysis
Analyses were done using the R software version 3.11 71 using the vegan package 72 , the rpart package 73 and DESeq2 package 30 .
The normality and homogeneity of variance of the environmental variables were tested using a Shapiro- to adjust the p value for multiple testing (Supplemental Fig. S3).
In order to harmonise data and allow comparison between samples, reads were transformed into relative abundance. Changes in the β-diversity of the ciliates at the community level between the past and recent samples were investigated using a Bray-Curtis dissimilarity matrix that was built based on the relative abundance data at the subclass level. Results were visualized on a NMDS (Non-metric Multidimensional Scaling analysis). To evaluate the change in dispersion and thus in diversity between the past and recent samples, the distances between the samples and the geometric median for each group ("recent" and "past") were calculated. The difference between the median of each group was then tested using the Wilcoxon rank sum test in order to evaluate the overall displacement of the recent and past samples. A hierarchical cluster analysis on Bray-Curtis distances using unweighted pair group method with arithmetic mean (UPGMA) was used to identify if a clear separation could be observed between and within recent and past samples. A SIMPER analysis (SIMilarity PERcentage) 78 was performed on the relative abundance data to calculate the contribution of each subclass to the overall Bray-Curtis dissimilarity between the recent and past samples. The most abundant species can have a high contribution even when they do not differ among groups, as they tend to display the highest variance 72 , as such, the proportion in SIMPER contribution, in average change between recent and past samples, as well as the total number of reads were compared to each other.
In order to evaluate the relative importance of known physical characteristic of the lakes (i.e. continuous variables: elevation, maximum depth, surface area of the lake, surface area of the watershed) to the amplitude of changes in the β-diversity, a univariate regression tree analysis was applied on the Bray-Curtis dissimilarity between the recent and past sample of each lake. For categorical data (i.e. Trophic Status), an analysis of variance (ANOVA) was used to compare the Bray-Curtis dissimilarity coe cients between categories.
The difference in abundance between the recent and past samples was also evaluated for each functional group and at the Genus level using the DESeq2 framework applied on the raw count data 30 .  Dendrogram illustrating the results from the hierarchical cluster analysis. The ciliate community composition the dominant Order or Class levels (i.e. Class_NA, if the Order level could not be assigned) of each sample is illustrated on the right side of the dendrogram and expressed as relative abundance (%).
Each sample is named by their corresponding lake code (cf. Table S3) followed by the code "_T" or "_B" indicating the recent (i.e. "recent", in blue) or past (i.e. "bottom", in red) samples respectively. Distribution of the Bray-Curtis dissimilarity coe cients between recent and past samples (black dot) relative to lake elevation (m). Fitted regression tree model (n =48 lakes) identi ed a split at 1400m and is represented by black lines (mean values). Gray shading represents the 95% con dence intervals around means.

Figure 5
Amplitude of change in the functional groups between the recent and past samples according to the feeding ecology (A) and limnetic habitat preferences (B). Magnitude of change is expressed in log2 fold change, as estimated by the DESeq2 analysis (n= 48 lakes). Dark green bars represent groups for which the change was found signi cant according to the two-sided Wald test corrected with the Benjamini and Hochberg method (p-value < 0.05). Horizontal lines show the standard error. Note: OTUs that were