Functional-based approaches are increasingly being used to study aquatic ecosystems as an alternative to traditional taxonomy-based approaches. Functional diversity is a biodiversity measure based on the ecological role of the species present in a community. Species-specific functional traits, or “what they do” , allows species to be defined by their interactions within an ecosystem  in terms of their ecological roles and how they interact with the environment and with other species .
Many recent ecological studies [2,4,5] suggest the importance of species ecological roles, and not just the number of taxonomic species, in the relationship between biodiversity and ecosystem functioning. This is a central concept if we are to understand and predict the resilience of a community to perturbations. If two species are deemed to be functionally alike and to occupy a similar trophic niche [6,7], the loss of one of those species is not likely to have an impact on the resource pool, as the other will increase its activity accordingly. The loss of functional diversity in the ecosystem is mitigated, as the species lost does not possess unique functional traits. Thus, the sum of organism functional traits within an ecosystem can be said to represent an indirect measure of its functional diversity [1,8]. Although the importance of functional diversity is widely recognized, there is no consensus on how to quantify functional diversity within a community, as relationships between the various indices have not yet been established .
Freshwater zooplankton play a key role in aquatic ecosystems in the transferal of biomass and energy from phytoplankton to top predators, e.g., [10,11]. While studies on phytoplankton have increasingly highlighted the importance of functional traits and functional classification in ecological studies, only a few have attempted this approach with zooplankton. Barnett et al.  applied the measure of functional diversity  to crustacean zooplankton communities. Quantitative functional traits considered were the C:N ratio, mean body size, and preferred food size range, while qualitative traits described the preferred habitat, food selectivity, and trophic position of each species.
Analysis of δ13C and δ15N stable isotopes is widely used to quantify food sources, trophic positions, and the interactions of organisms, as the δ13C of a consumer can infer the assimilated source of dietary carbon [13,14] and δ15N the trophic role [14,15]. Considering a functional-based perspective, δ13C and δ15N could identify ecological relations among taxa, such as competition and predation, ecological niche, habitat preference, and taxa redundancy, and help to define functional groups in an ecosystem. In this study, we propose the use of δ13C and δ15N stable isotope analysis to quantify some of the “qualitative functional traits”  for pelagic crustacean zooplankton taxa in Lake Maggiore, a large, deep subalpine lake in Italy. As species composition, diversity, and biomass of zooplankton can change significantly seasonally, especially in temperate lakes , the seasonal variation of δ13C and δ15N was used to define seasonal patterns in habitat preference, food sources, and taxa trophic position. This allowed an interpretation of taxa redundancy and a hypothesis of bottom-up and top-down mechanisms potentially driving the observed changes.
One use of stable isotopes is to determine the proportional contributions of several sources in a mixture. An example of a source proportion calculation includes the determination of various food sources in an animal’s diet . Linear mixing models are used to estimate proportions for two sources using isotopic signatures for a single element (e.g., δ13C), or for three sources using isotopic signatures for two elements (e.g., δ13C and δ15N; ). In this study, we have used a linear mixing model [18,19] in order to discriminate the relative contribution of different preys in the diet of zooplankton consumers such as calanoid and cyclopoid copepods and the predatory cladoceran Leptodora kindtii. This further contributed to the understanding of food preference and taxa trophic position in the zooplankton community of Lake Maggiore.
2. Materials and Methods
Lake Maggiore (45°57′ N 8°32′ E 3°47′ W) is the second deepest (dmax 370 m) and largest (area 212.5 km2, volume 37.5 km3) subalpine lake in Italy. Being phosphorus-limited (TPmax ca. 10 μg L−1), the lake is oligotrophic and has recovered from eutrophication of the late 1970s [20,21].
Except for September, vertical zooplankton hauls from the surface to a depth of 50 m were collected. This follows the standard routine sampling for deep subalpine lakes, in which samples are collected within the upper 50 m depth, as previous research on the vertical distribution of zooplankton showed that this is the water layer in which zooplankton live . Monthly samples were collected from April to November 2009, when total zooplankton biomass was ≥3 mg·m−3, using a wide-mouth 450 μm mesh zooplankton net of diameter 0.58 m, filtering 13 m3 lake water from three pelagic stations (G: 45°58′30″ N 8°39′09″ E, B: 45°54′28″ N 8°31′44″ E, L: 45°49′70″ N 8°34′70″ E) .
Zooplankton samples for the quantification of biomass were collected with a Clarke-Bumpus plankton sampler of a 126 μm mesh size and fixed in ethanol 96% to estimate the taxa-specific population density (ind·m−3) and standing stock biomass (dry weight, mg·m−3) . Organisms for isotopic analyses were kept overnight in filtered (1.2 μm GF/C filters) lake water for gut clearance, before sorting into taxa and quantities suitable for isotopic analyses. The taxa analyzed were Daphnia longispina galeata gr., Eubosmina longispina, Diaphanosoma brachyurum, Bythotrephes longimanus, Leptodora kindtii, adults of the calanoid copepods Eudiaptomus padanus and Eudiaptomus gracilis, and of the cyclopoid copepods Mesocyclops leuckarti and Cyclops abyssorum.
Zooplankton δ13C and δ15N isotopic signatures were referred to that of the pelagic baseline, which was expressed by the primary consumer, Daphnia longispina galeata gr. This choice came from previous stable isotope studies in Lake Maggiore , showing that Daphnia’s δ13C signature in the different seasons was closely correlated with the signature of seston (r = 0.86; p < 0.01; N = 13), confirming that Daphnia was an appropriate proxy for the pelagic baseline against which the carbon isotopic signals of other zooplankton can be compared. Δ13C was used to detect seasonal changes in taxa specific feeding behavior and assess the origin of carbon sources fueling the pelagic food web.
p = (δ13Cpredator − δ13Cprey2)/(δ13Cprey1 − δ13Cprey2); q = 1 − p
p = ((δ15Nprey3 − δ15Nprey2)(δ13Cpredator − δ13Cprey2) − (δ13Cprey3-δ13Cprey2)(δ15Npredator − δ15Nprey2))/
((δ15Nprey3 − δ15Nprey2)(δ13Cprey1 − δ13Cprey2) − (δ13Cprey3 − δ13Cprey2)(δ15Nprey1 − δ15Nprey2));
Nutrient enrichment of freshwaters is a worldwide challenge . In combination with intensified climate warming, these anthropogenic changes threaten aquatic biodiversity and ecosystem services [2,3]. Ecological communities and their responses to environmental stressors, such as eutrophication in aquatic systems, can be investigated with a functional approach, where interest is put on species’ ecological roles. For example, feeding traits, habitat preferences, reproduction, or morphological attributes (e.g., body size) can reflect certain ecological functions . The concept of functionality in ecological communities allows a comprehensive understanding of how environmental changes alter ecosystems through biological functions rather than just taxonomic composition . This kind of mechanistic approach on aquatic systems may reflect important ecosystem level processes, for example, changes in productivity and trophic structure of lakes. In the functional approach, functional diversity (FD) is a biodiversity measure, which takes into account the variety of biological functions of species and may allow for a more holistic understanding of environmental changes and ecosystem responses [4,6].
With major losses in global biodiversity , paleolimnological data sets can aid our understanding of the relationships among functional diversity and ecosystem productivity, climate change, and trophic dynamics [8,9,10]. Since key members of aquatic communities (or their traces) are preserved as fossils in lake sediments, paleolimnology can be used to evaluate ecosystem functions and ecological resilience in lakes. Further, lake sediment archives are advantageous for functional classification and ecosystem level responses to environmental changes, because time lags often prevent detection during short-term observations [11,12].
Here, we continue the application of paleolimnological research on Lake Maggiore in northern Italy (Figure 1) to investigate the lake as a “natural laboratory” with its well-documented history of eutrophication and re-oligotrophication during the last century. This subalpine lake is Italy’s second largest and deepest lake and part of the long-term limnological monitoring in Europe . Naturally oligotrophic and phosphorus-limited Lake Maggiore eutrophied in the 1960s as a result of nutrient loading from the catchment and wastewater discharge [14,15]. Phosphorus concentrations started to increase and the peak of nutrient enrichment occurred during the late 1970s when total phosphorus at winter mixing (TPmix) was 31 µg L−1 and the lake became mesotrophic (Figure 2). After that, due to enhancements in wastewater treatment, recovery of the lake proceeded and TP decreased back to oligotrophic levels (~10 µg L−1) during the early 1990s. Previous research on Lake Maggiore plankton has indicated that the aquatic communities are highly responsive to nutrient status, trophic dynamics, and climate [15,16,17,18]. For example, Lake Maggiore zooplankton has exhibited major changes in taxonomic composition, body size, and population density under eutrophication and re-oligotrophication [15,17,18].
The aim of the current research was to examine responses to eutrophication and the subsequent limnological recovery of Lake Maggiore by the cladoceran communities. We analyzed fossil cladoceran communities for their taxonomic composition, functional characterization, and functional attributes of FD and Daphnia ephippia length (DEL) in a sediment core covering the years of pre-eutrophication, eutrophication and recovery (1944–2010). We aimed to identify the main environmental forcings on the long-term succession of cladoceran communities, functional assemblages, and functional attributes, and discuss the roles of bottom-up versus top-down controls in ecosystem functioning during the eutrophication and re-oligotrophication of Lake Maggiore.
2. Materials and Methods
A modified Wilco box corer (liner internal surface area 28 cm2, shaft depth 50.5 cm) was used to collect a 34-cm core in the Pallanza Basin of Lake Maggiore (45°54.76 N, 8°32.96 E, z = 98 m) on 11 February 2010 (Figure 1). The intact core was stored dark at 4 °C until processing. The core was vertically sliced at 1-cm intervals and the outer edge of each segment was discarded. A 10 g subsample was removed from each interval for 137Cs dating at University of Applied Sciences (Holland). Increases in the 137Cs activity were interpreted as increased fallout of 137Cs from nuclear activity (Figure 3a). An age–depth model was created by interpolating the surface sediment sample (0 cm, year 2010), increased 137Cs activities of Chernobyl fallout (1986) at 14 cm, and the nuclear weapons testing horizon (1963) at 25 cm (Figure 3b). As the model is based on constant sedimentation rate between these time horizons it should be interpreted cautiously.
The remaining portion of each sediment interval was prepared for fossil Cladocera analysis following the standard methods . First, the samples were heated in 10% KOH to deflocculate sediment and then sieved through a 76-µm mesh, which is adequate to retain cladoceran zooplankton fossil remains and the remains of smallest taxa of Chydoridae (e.g., Alonella nana and Coronatella rectangula, Table 1). The residue was stained pink with safranine. The fossil Cladocera remains (carapaces, headshields, postabdomens and ephippia) were identified and enumerated with a light microscope (magnifications 100–400X) and the most abundant body part was chosen to represent the number of individuals of each species. A minimum of 100 fossil individuals were counted from each subsample. Relative abundances of individual taxa were used to determine the community composition in the core samples. In addition, to examine functional characterization of the community (i.e., functional assemblages) taxa were assigned to the groups (Table 1): predators, large filter feeders, small and intermediate filter feeders (hereafter small filterers), globular epibenthos, and oval epibenthos based on a previous functional grouping of Cladocera . Size (i.e., body length) from base to apex (spine excluded) of encountered fossil Daphnia ephippia was measured with a Zeiss microscope at 100X magnification equipped with a camera and analyzed with Image pro express 5 software to estimate the mean Daphnia ephippia length (DEL). Ephippia were measured from 32 sediment samples and number of size measurements per sediment subsection varied from 1 to 73 (mean 13).
Principal component analysis (PCA) was used to summarize temporal succession of cladoceran taxonomic communities (compositional gradient <1.5 SD, standard deviation units) and functional assemblages (<1.0 SD). The response data were square root transformed for PCA. In addition, redundancy analysis (RDA) was used to analyze relationships between functional assemblages and cladoceran taxonomic communities and limnological variables during the period of continuous environmental monitoring (since 1978). Chlorophyll-a concentration (chl-a), total phosphorus at winter mixing (TPmix), water temperature of the euphotic zone 0–20 m (T0–20m), Bythotrephes longimanus abundance, and total pelagic fish catch were included as environmental variables (Figure 2). Environmental variables were forward selected and the significance (p ≤ 0.05) of each variable was tested with Monte Carlo permutations (999). Cladoceran functional diversity (FD) was evaluated with Rao’s FD index , i.e., Rao’s quadratic entropy. For the index, each cladoceran taxa was assigned with qualitative functional character including body size (small < 500 µm, intermediate 500–1000 µm, large > 1000 µm), body shape (elongated, oval, globular), feeding type (filterer, scraper-detritivore, predator) and microhabitat (pelagic, benthic, attached to vegetation, Table 1) . This characterization was based on ecological data available for cladoceran taxa [21,22] and the characters were inserted as functional character present (1) and absent (0). Multivariate analyses (PCA and RDA) and analysis of FD were performed with Canoco 5 software . Tukey’s pairwise comparisons were utilized to indicate differences in FD and DEL during the reference (pre-1960), eutrophication (1960–1985), and recovery (post-1985) periods. These analyses were performed with PAST software . Segmented regression analysis was utilized to detect statistically significant breakpoints (minimum confidence level of 95%) in FD and DEL. The best breakpoint was selected based on maximizing the statistical coefficient of explanation and performing tests of significance with SegReg program .
We detected 30 cladoceran taxa in the Lake Maggiore sediment core (Table 1). The most abundant taxa were Bosminidae, including Eubosmina longispina-type (n = 34, mean percent abundance 63.5%), Eubosmina coregoni-type (n = 26, 5.7%), and Bosmina longirostris (28, 2.6%) and Chydoridae (chydorids), including Chydorus cf. sphaericus (n = 31, 11.3%) and Alona affinis (n = 34, 6.7%). In the early core (until 1960s; reference period) the communities were dominated by A. affinis and E. longispina-type (~60%) together with lower abundances of B. longirostris, Leptodora kindtii, Sida crystallina and several less common chydorids (e.g., Alona quadrangularis, Paralona pigra, Eyrycercus spp. and Acroperus harpae; Figure 4). Chydorus cf. sphaericus started to increase in the early 1960s and reached maximum abundance (~40%) and dominance during early 1970s. Bythotrephes longimanus, which were abundant in zooplankton samples starting in the 1990s (Figure 2), were not detected in the fossil community (Table 1). In functional assemblages, succession proceeded from oval epibenthos during pre-1960s to globular epibenthos (~40–50%) between 1960 and 1975, and small filterers (80–90%) after that until the top core (Figure 4). Predators and large filterers were scarce (<10%) and occurred mostly prior to 1960s.
PCA for cladoceran taxonomic communities resulted in eigenvalues 0.514 for PCA axis 1 and 0.217 for PCA axis 2. Cumulative percentage of variance explained by the PCA were 51.4% and 73.2% for PCA axes 1 and 2, respectively. Samples in the reference period (pre-1960) had positive PCA axis 1 scores with increasing PCA axis 2 scores (Figure 5a). Samples from the eutrophication period (1960–1985) had reducing scores along PCA axes 1 and 2. Recovery period samples (post-1985) clumped together at the negative end of PCA axis 1. PCA for functional assemblages had eigenvalues of 0.763 for PCA axis 1 (cumulative % of variance 76.3) and 0.169 for axis 2 (93.2%). Sample scores drifted from positive axis 2 values to negative prior to 1960 and from positive to negative axis 1 values during eutrophication (Figure 5b). Most recent samples of the recovery period had negative axis 1 values and close to zero axis 2 values.
RDA for functional assemblages resulted in eigenvalues 0.222 for RDA axis 1 and 0.1585 for axis 2 and all the environmental variables explained 40.6% of the variance in the assemblages. RDA identified TPmix (38.3%) as the single significant (p < 0.05) environmental factors explaining variance in functional assemblages (Table 2). RDA for taxonomic communities resulted in eigenvalues 0.2158 and 0.0778 for axes 1 and 2, respectively and forward selection did not results in any statistically significant results (Table 2). Bythotrephes longimanus abundances, however, explained most variance in the communities with 33.4% (p = 0.0600).
FD varied between 1.7 and 3.2 in Rao’s FD index (Figure 6). The highest FD occurred in the early core until 1960 after which FD decreased until the early 1980s. FD increased slightly between 1990 and 2000 but was reduced post-2000. DEL measurements varied between 470 and 690 µm (Figure 6). DEL was largest (~650 µm) in the early to mid-core until 1970s and then started to decrease and remained consistent (~550 µm). Tukey’s test indicated significant differences in FD and DEL between reference and eutrophication periods, and reference and recovery periods (Table 3). No significant differences in FD or DEL were found between eutrophication and recovery periods. SegReg identified single significant breakpoints for FD at 15.64 cm (early 1980s) and DEL at 20.14 (early 1970s, Figure 6).
Bottom-up (food, habitats) and top-down (predators) controls drive cladoceran communities in Lake Maggiore, based on trophic dynamics inferred from sediments dating from 1943 to 2002 . Here we present more refined functional analysis and include the most recent decade of fossil records (1944–2010). Cladoceran communities succeeded from Alona affinis to Chydorus cf. sphaericus, and later to Eubosmina longispina-type dominance with a functional shift from large filterers (e.g., Sida, Daphnia) and oval epibenthos (e.g., A. affinis) to globular epibenthos (C. cf. sphaericus) and most recently to small filterers (Eubosmina, Figure 4). The biotic changes occurred gradually and even prior to the eutrophication as small filterers and epibenthos exhibited shifts prior to 1960s but predators and large filterers responded promptly at ~1960 (Figure 4). This suggests that there was a difference in timing of the response among the functional groups and species, likely related to changed resources just prior to major nutrient loading. The PCAs indicated a clear temporal shift in the taxonomic communities and functional assemblages through the reference–eutrophication–recovery phases with a new stable state during the recent years of re-oligotrophication (Figure 5). The opposite drift of PCA axis 2 scores in taxonomic (from negative to positive scores, Figure 5a) and functional assemblages (from positive to negative, Figure 5b) during the reference–early eutrophication may suggest that diverse control mechanisms drove changes in taxonomic vs. functional communities. In agreement, RDA forward selection statistics suggest that functional assemblage shifts were associated with bottom-up controls (Table 2), since they were mainly predicted by TPmix during the time period of the core covered by limnological monitoring (since late 1970s, Figure 2). Taxonomic communities were best explained, although only marginally significant, with top-down control, i.e., Bythotrephes abundance (Table 2). Bythotrephes is known to regulate zooplankton, especially Daphnia as its main prey item, in Lake Maggiore and impact their distribution, abundance and phenology [18,26,27].
In general, caution should be taken when applying the ecological patterns reported here as whole-lake patterns, as the current sediment core was sampled from the Pallanza Basin (Figure 1). Nevertheless, the trends reported here are consistent with the previous research from another core . Fossil remains (e.g., mandibles, caudal spines, resting eggs) of Bythotrephes longimanus, although abundant in zooplankton samples in Lake Maggiore (Figure 2, also in the Pallaza Basin), were not detected from the current core (Table 1) or in the previous sediment-based research  and the reason for that remains unexplained. Typically, Bythotrephes remains (caudal spines and sometimes resting eggs) are well recovered from sediments and their accumulation has been used to estimate presence and abundance [28,29,30,31]. Bythotrephes remains have been found from sediments in large and deep subalpine lakes, such as Lake Garda in northern Italy , in similar type of geo-limnological settings than Lake Maggiore. The use of a larger (76-µm) than the commonly used ~50-µm sieve  for processing sediment samples for fossil cladoceran analysis in the current study is not a reasonable explanation for the lack of Bythotrephes remains, since the previous studies have used even a 250-µm aperture sieve [28,29,30]. Accordingly, it may be that the sediment sampling site in the Pallanza Basin was not a representative location for Bythotrephes fossils to preserve. It is further possible that some still unidentified limnological factor in Lake Maggiore prevents Bythotrephes fossils for preserving in the sediments.