**Long-Term Changes in the Zooplankton Community of Lake Maggiore in Response to Multiple Stressors: A Functional Principal Components Analysis**

#### **Andrea Arfè 1,\*, Piero Quatto 2, Antonella Zambon 3, Hugh J. MacIsaac <sup>4</sup> and Marina Manca <sup>5</sup>**


Received: 25 January 2019; Accepted: 25 April 2019; Published: 8 May 2019

**Abstract:** We describe the long-term (1981–2008) dynamics of several physico-chemical and biological variables and how their changes may have influenced zooplankton structure in Lake Maggiore (Italy). Data was available for the 1981–1992 and 1995–2008 periods. Standardized time-series for temperature and total phosphorus (TP), chlorophyll-α, phytoplankton density (cel m<sup>−</sup>3), and cell size (μm3), as well as zooplankton structure (Copepoda, Cladocera, and Rotifera density, ind m−3) were smoothed using penalized B-splines and analyzed using Functional Principal Components (FPCs) to assess their dominant modes of variation. The first four FPCs explained 55% of 1981–1992 and 65% of 1995–2008 overall variation. Results showed that temperature fluctuated during the study period, particularly during 1988–1992 with a general tendency to increase. TP showed a declining trend with some reversions in the pattern observed in the years 1992, 1999, and 2000. Phytoplankton estimators and chlorophyll-α concentration showed a variable trend along the study period. Zooplankton groups also had a variable trend along the study period with a general increase in density of large carnivorous (mainly *Bythotrephes longimanus*) and a decrease of large herbivorous (mainly *Daphnia*), and a similar increase in the ratio of raptorial to microphagous rotifers. Our results suggest that the lake experienced a strong trophic change associated with oligotrophication, followed by pronounced climate-induced changes during the latter period. TP concentration was strongly associated with changes in abundance of some zooplankton taxa.

**Keywords:** B-Splines smoothing; Functional Data Analysis; limnology; monitoring ecological dynamics; oligotrophication; zooplankton; phytoplankton

#### **1. Introduction**

Lakes occur at a pivotal position in global landscapes, receiving inputs of both organic and inorganic matter, and generally reflect events occurring in their watersheds. Some of these events stress aquatic life, including soil and water acidification, soil erosion, loss of base cations, release of heavy metals or organic compounds, and application of essential nutrients capable of stimulating primary productivity [1–11]. Superimposed on these changes, climate warming directly impacts lakes via alteration of species' metabolic processes and indirectly by modifying food web interactions [12,13]. The scale of these stressors is influenced by many factors, including loading rate, basin morphometry and depth, and water residence time.

Many stressors interact in a manner that can be difficult to predict [14]. In part, this difficulty stems from the different possible responses by species or entire taxonomic groups to stressors, which may interact additively, synergistically, or antagonistically [2,15]. Mechanistically, negative co-tolerance by a species to interacting stressors yields either additive or synergistic responses with amplified consequences, whereas positive co-tolerance produces antagonistic responses with damped consequences [2]. The onset and termination of stressors often vary temporally, thus it may be possible to use long-term dynamics of lake food webs to discern changes owing to individual stressors or to interactions between them [16]. In a meta-analysis, Jackson et al. [15] determined that warming and nitrification resulted in additive net effects in freshwater systems, while warming and any other stressor was antagonistic.

While excellent long-term analyses exist for individual stressors, such as cultural eutrophication [17], there exists a relative dearth of similar studies to disentangle individual effects in systems subjected to multiple stresses. Long-term studies are essential to fully understand responses by lake ecosystems to perturbations. Schindler [18] warned against extrapolation of short-term or small-scale studies to reveal processes important to cultural eutrophication, noting that only long-term, whole-lake experiments and case histories of recovered lakes can reliably inform abatement policy. Despite this and other calls for long-term studies of lakes [6], there remains a surprising paucity of such studies (but see, for example, References [9–22]).

The objective of this study was to analyze the long-term dynamics of several physico-chemical and biological variables of Lake Maggiore to identify major trends and how this variation impacts zooplankton abundance and community structure. Lake Maggiore is a large, deep lake in the subalpine area of northern Italy, which sustained a wide-ranging series of major perturbations during the 20th century. These changes included cultural eutrophication, organic chemical contamination, commercial fishing, climate warming, and species introductions. To characterize the variation in time of several biological and physico-chemical parameters of Lake Maggiore, we adopted a Functional Data Analysis (FDA) approach [23] based on penalized B-spline smoothing and Functional Principal Components Analysis (FPCA). This approach, seldom considered within the ecological literature, has several features that make it especially suited for the analysis of limnological data.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Lake Maggiore is an oligomictic, subalpine (194 m above sea level) lake, largely contained within Northern Italy but shared in the Northern end of the basin with Switzerland. The lake has a surface area 212.5 km2 and a maximum (mean) depth of 370 m (177 m), a product of tectonic-glacial activity. Lake Maggiore is one of the best-studied lakes in the world, with a long-term record of physico-chemical and phytoplankton and zooplankton data dating back to the mid- or even early 1900s for some variables [24]. The drainage basin (including the lake area) covers 6599 km2, yielding a drainage basin/lake area ratio of 31:1. This high ratio along with unusually steep hillslopes substantially affect the lake's hydrology and many environmental variables. Geological features of the watershed have been influenced by alpine orogenesis and glaciation, and thus are complex. The major topographic feature is a narrow steepsided valley through the Piedmont and Lombardy regions of Italy.

As with many lakes in Europe, P-limited, oligotrophic, Lake Maggiore experienced cultural eutrophication from the late 1960s throughout the 1970s; peak total phosphorus (TP) and chlorophyll-α concentrations were detected at the end of 1970s, with values as high as 40 μg m−<sup>3</sup> of total phosphorus in-lake concentration at winter mixing (TPmix) of 1979, attesting to mesotrophic conditions of the lake [25]. Gradual implementation of sewage treatment plants, and a reduction of the phosphorus content in detergents, led to a gradual decline in TP and eutrophication reversal. A time lag was reported in the response of plankton communities to TP decline: as expected, changes in zooplankton and phytoplankton became evident only after TPmix concentration declined below a threshold of 15 mg m−<sup>3</sup> in 1988 [26]. According to Sas [27], this is the level at which the total phosphorus becomes limiting to phytoplankton growth. The trophic evolution of the lake is most strictly represented by Cladocera rather than copepods [28,29].

#### *2.2. Data Collection*

In this study, we describe dynamics of biological and physico-chemical variables of Lake Maggiore's ecosystem for the periods 1981–1992 and (separately) 1995–2008. Specifically, we considered physical and chemical variables, i.e., water temperature and total phosphorus, as well as biological variables, i.e., phytoplankton and zooplankton population density, biomass, and cell size (measured by cellular volume).

Data for both periods were collected as part of the long-term monitoring of Lake Maggiore, funded since the late 1970s by the International Commission for the Protection of Swiss-Italian Waters (Commissione Internazionale per la Protezione delle Acque Italo-Svizzere, CIPAIS) and published in annual reports (available at www.cipais.org). Phytoplankton and zooplankton measurements were also obtained from the "Plankton and the pelagic food web" research project funded by the National Research Council (CNR). All measurements refer to a single sampling station (45◦58'30" N; 8◦39'09" E) at the lake's maximum depth, which is representative of the pelagic environment [30,31].

In more detail, monthly data on temperature at 0–20 m depth, the water layer representative of the euphotic zone, were based on annual reports of the CNR-ISE meteorological station [32–37]. Total phosphorus concentration was based on a maximum of 13 measures at different depths and was included as indicative of changes in lake trophic status. Data on phytoplankton density and cell size as well as chlorophyll-α were based on integrated phytoplankton samples collected every other week with a 1.5 L van Dorn bottle within the 0–20 m layer. Samples for chlorophyll analysis were filtered through GF/C glass fibre filters (about 1μm pore size). The filters were then stored on silica gel at −20 ◦C. After about 2–3 weeks the filters were mechanically ground and transferred to acetone 90%. The absorbance of the pigments was measured in a spectrophotometer (Perkin-Elmer Lambda 6); the chlorophyll-α and phaeophytin concentrations were calculated in accordance with Lorenzen [38]. Phytoplankton samples were preserved with acetic Lugol solution and counting was made following the Utermöhl method [39] at 400x in an inverted microscope [40]. Species identification and nomenclature followed the more recent monographs of the series Sußwasserflora von Mitteleuropa, established by A. Pascher (Gustav Fisher Verlag, and Elsevier, Spectrum Akademischer Verlag), specific manuals of the series Das Phytoplankton des Sußwassers, established by G. Huber-Pestalozzi (E. Schweizerbart'sche Verlagsbuchhandlung), and specific papers [41]. Zooplankton population abundance data (*Bythotrephes longimanus,* Leydig 1860; *Eubosmima longispina,* Leydig 1860; cyclopoid copepods, *Daphnia longispina-galeata,* Wagler, 1937; *Diaphanosoma brachyurum,* Liévin 1848; and *Leptodora kindtii.* Focke 1844; and Rotifera) were obtained monthly with a 76-μm nylon net Clarke-Bumpus plankton sampler, towed at a constant speed of ca 3 km h<sup>−</sup>1, along sinusoidal hauls from the surface to 50 m depth, i.e., in the water layer where zooplankton live [30]. Samples included at least 1000 L of water and were fixed in 4% buffered formaldehyde before counting and identifying the genus or species under a microscope at 6.3x. Taxa identification was based on volumes of the Identification Guides to the Plankton and Benthos of Inland Waters (formerly "Guides to the Identification of the Microinvertebrates of the Continental Waters of the World", H. J. F. Dumont coordinating Editor. Backhuys publisher). Further details on data collection are reported elsewhere [29,31].

#### *2.3. Smoothing Data through Penalized B-Spline Expansions*

Although biological and physico-chemical parameters are typically measured at discrete instants, their temporal dynamics can be better represented as a smooth function varying in some continuous time interval [*a*, *b*]. These smooth representations are of main interest in FDA, where they are typically recovered by penalized splines approximations. Briefly, a *spline basis* is a set of known functions, which can be used to approximate any other function with arbitrary precision. One such set frequently used in FDA consists of smoothly-joined piecewise polynomial functions called *B-splines*. These allow researchers to represent any collection of *m* time series *si*(*tik*) (indexed by *i* = 1, ... , *m* and observed at distinct time points *ti*1, ... , *tiKi* ) with smooth functions *xi*(*t*) of the form

$$x\_i(t) = \sum\_{j=-1}^{B} c\_{ij} b\_j(t),\tag{1}$$

where *t* belongs to the time-interval [*a*, *b*] and *B* is the considered number of B-splines *bj*(*t*). The number *B* is uniquely determined by the polynomial degree and the chosen number of *knots,* i.e. the points where the different polynomial pieces are joined [23]. The coefficients *cij*, which provide the best smooth approximation of the observed time-series, are obtained by minimizing the *Penalized Sum-of-Squares Error* (*PSSE*)

$$PSSE\_i \ = \sum\_{k=1}^{K\_i} \left[ s\_i(t\_{ik}) - x\_i(t\_{ik}) \right]^2 + \lambda\_i \int\_a^b \left[ \frac{d^2 x\_i(t)}{dt^2} \right]^2 dt. \tag{2}$$

*PSSEi* is the sum of two terms: the first represents the sum of squared approximation errors, which decreases as the B-spline approximation better fits the data; the second is linked to the so-called "*strain energy*" (much like stretched elastic bands, the more the curve *xi*(*t*) is "wiggly", the higher its strain energy). Hence, minimizing *PSSEi* results in a B-spline approximation whose fit to the data is a compromise between low approximation error and "wiggliness" of the resulting curve (thus avoiding overfitting). This trade-off is controlled by the *penalization coe*ffi*cient* λ*<sup>i</sup>* ≥ 0, which is usually chosen by the generalized cross validation method [23].

#### *2.4. Functional Principal Components Analysis (FPCA)*

FPCA is a generalization of standard Principal Components Analysis (PCA) to the situation where data is represented by smooth continuous curves (possibly obtained by penalization techniques), as in our setting [23]. Specifically, *Functional Principal Components (FPCs)* are square-integrable orthonormal functions *ej*(*t*) that best represent, in a least-squares sense (see the Appendix A), the *m* smooth curves *xi*(*t*) by means of expansions such as

$$\mathbf{x}\_{i}(t) \triangleq \overline{\mathbf{x}}(t) + \sum\_{j=-1}^{n} f\_{ij} \mathbf{v}\_{j}(t),\tag{3}$$

where *x*(*t*) = <sup>1</sup> *m m <sup>i</sup>* <sup>=</sup> <sup>1</sup> *xi*(*t*) is the mean function, *n* is the number of FPCs *ej*(*t*), and *fij* is the (*FPC*) *score* for the *i*-th time series with respect to the *j*-th FPC [23]. As in standard PCA, FPCs are chosen so as to capture the greatest possible variation in the observed data. More precisely, each FPC *ej*(*t*) (with *j* = 1, ... , *n*) is chosen in order to maximize the corresponding explained variation, i.e., the variance σ2 *<sup>j</sup>* of the *<sup>m</sup>* values *<sup>f</sup>*1*j*, ..., *fmj*, with <sup>σ</sup><sup>2</sup> <sup>1</sup> <sup>≥</sup> <sup>σ</sup><sup>2</sup> <sup>2</sup> ≥···≥ <sup>σ</sup><sup>2</sup> *<sup>n</sup>*. Thus, for example, the first and the second FPC *e*1(*t*) and *e*2(*t*) represent the mode of variation from the overall mean *x*(*t*) associated with, respectively, the greatest and the second greatest possible variation between the curves *x*1(*t*), ... *xm*(*t*). On the other hand, the score *fij* measures how well the *j*-th FPC captures variation over time of the *i*-th time series (scores of higher absolute magnitudes are associated to the FPCs capturing the greater amount of variation). In more detail, a small *fij* score signifies that the *j*-th FPC does not capture the time variation of the *i*-th time series. Instead, a high positive (negative) *fij* score signifies that the time variation of the *i*-th time series is captured by the same (opposite) trend than the *j*-th FPC. For example, if the *j*-th FPC shows an increasing global trend, then a positive (negative) *fij* means that the *i*-th time series increases (decreases) over time.

When the FPCs have been identified, their interpretation can be aided by plotting the two curves obtained by adding and subtracting from the overall mean *x*(*t*) the *scaled FPCs*, i.e., the curves *x*(*t*) ± 1.5σ*jej*(*t*) for each FPC *ej*(*t*). Furthermore, the VARIMAX strategy (commonly used in standard PCA) can be generalized to the FPCA setting. According to this strategy, the FPCs are rotated so that each component is associated to only a small number of high scores. Rotating the FPCs according to this strategy can therefore aid their interpretation [23].

#### *2.5. Describing Lake Maggiore's Dynamic by Functional Data Analysis*

The above-described data analysis methods were implemented to assess long-term temporal dynamics of Lake Maggiore as follows. First, the time-series for the considered variables were smoothed by means of penalized B-splines. Uniformly spaced knots (one about every 30 days) were considered for both the 1981–1992 or 1995–2008 periods, in each case for a total of B = 200 B-splines. This choice ensured that smoothed curves had enough flexibility to represent month-level changes in the underlying variables. Second, FPC analysis was performed on the B-spline smoothed time series of the standardized variables. Specifically, each time series *si*(*tik*) was standardized as (*si*(*tik*) − *si*)/τ*<sup>i</sup>* where *si* and τ*<sup>i</sup>* are respectively the mean and standard deviation of the *Ki* observations *si*(*ti*1), ... ,*si tiKi* . Standardization ensured time-series were expressed on the same dimensionless scale, allowing direct comparisons between them. Smoothed standardized time-series (using B = 200 B-splines) were then used to extract VARIMAX-rotated FPCs. All components that explained at least 10% of total variability of the smoothed standardized time series, i.e., the sum of the σ<sup>2</sup> *<sup>j</sup>* [23], were extracted. All analyses were performed in the statistical software R (version 3.0) using the "fda" library [23].

#### **3. Results**

#### *3.1. Long-Term Limnological Change*

Thermal conditions changes in Lake Maggiore over the study period. Mean summer temperature across the 0–20 m depth averaged 21.2 ◦C between 1981 and 1992 and 21.9 ◦C between 1995 and 2008, as also observed in previous studies [42]. Lake Maggiore underwent a gradual warming of the epilimnion in the decade 1988–1998, according to the trend of increasing heat content of the deep Italian lakes, pointed out by Carrara et al. [43]. The period was characterized by strong fluctuations, in particular between 1988 and 1992. Since 1999, the increase of water temperature was less evident, probably because of the occurrence of particular hydro-meteorological mechanisms reducing the heat content of the water column [43,44]. Lake Maggiore also experienced a dramatic reduction in total phosphorus (TP) concentration over the study period. Mean TPmix concentration declined from 10.9 μg L−<sup>1</sup> between 1981 and 1992, to 8.4 μg L−<sup>1</sup> between 1995 and 2008. An exceptionally high value, of 17 μg L<sup>−</sup>1, was detected in 1991, after the lake's complete winter overturn. Over a general increase, a second peak value was recorded in 1999, after another complete overturn. The overall increase to stable values of 11 μg L−<sup>1</sup> after 2000 could be mainly explained by meteorological-climatic conditions, rather than by an eutrophication reversal [45]. Eutrophication abatement results were mixed, however, as mean chlorophyll-α concentration fell sharply (from 4.7 to 3.4 μg L<sup>−</sup>1) between these periods, though results for phytoplankton biomass were less clear (1830 vs. 1181 mm3 m<sup>−</sup>3, respectively).

#### *3.2. Description of the Smoothed Time-Series*

The observed and smoothed time-series of the considered variables are represented in Figure 2. A visual inspection of the graphs in Figure 2 suggests that the smoothed curves have a good fit to the overall trend of the corresponding variables. Major variations in phosphorus concentration, phytoplankton cell size and chlorophyll-α concentration during 1981–1988 matched an increase in variability and mean density of *Bythotrephes longimanus,* and a decrease in between-year variability and mean annual population densities of *Bosmina* (mainly, *Eubosmina longispina*) and *Daphnia* (*D. longispina-galeata* gr.). The increased within-year variability in *Bythotrephes* observed during the 1990s was maintained across the second time period, whereas *Daphnia* and *Bosmina* variability and population density increased between 2002–2008, along with *Diaphanosoma brachyurum*. During 2002–2008, variability of diaptomid copepods seems to have decreased slightly relative to the periods

1981–1992 and 1995–2000. Rotifer (and to lesser extent) *Leptodora kindtii* density were characterized by increased variability during the last five years of the 1995-2008 period. In addition, rotifer communities changed dramatically through time, with the mean ratio of raptorial to microphagous species abundance decreasing from 1.5 ± 2.0 (mean ± st.dev.) during 1981–1992 to 0.5 ± 0.6 (mean ± st.dev.) during 1995–2002 [42,46].

Phytoplankton density was highly variable during the 1996–1998 interval. Phosphorus concentration, phytoplankton cell size, chlorophyll-α concentration, *Bosmina* and *Daphnia* population density all exhibited a more variable trend during 1981–1992 than during 1995–2008 using both observed and smoothed time-series data (Figure 2). Conversely, phytoplankton density, *Bythothrephes*, cyclopoid copepods, and *Diaphanosoma* population density were more variable between 1995–2008.

Panel (**B**): Phytoplankton-related variables.

**Figure 1.** *Cont.*

**Figure 1.** Observed and smoothed time series for all physico-chemical and biological variables. Continuous line = smoothed curve; circles = original data points. Panel (**A**): Physico-chemical variables. Panel (**B**): Phytoplankton-related variables. Panel (**C**): Zooplankton-related variables. Data were not available for 1993-1995. Points are monthly means, lines are the smoothed value.

#### *3.3. Extracted Functional Principal Components*

For both periods, the first four rotated FPCs were retained, and explained 55% and 65% of the overall variation of smoothed standardized variables for 1981–1992 and 1995–2008 periods, respectively. The effect of each FPC (denoted FPC1-4 in decreasing order of explained variation for both periods) on the trend of considered variables is represented in Figure 2. This shows, for the two periods, the overall

mean (continuous line) of the smoothed standardized variables, along with curves representing the effect of adding (dashed line) or subtracting (dotted line) each scaled FPC on the mean. The FPC1 corresponds to an overall increasing trend in both the 1981–1992 and 1995–2008 periods, as in each case the dotted (increasing) and dashed (decreasing) lines cross the graph of the overall mean as time increases (Figure 2a). In addition, the FPC1 for the 1981–1992 period also corresponds to more extreme of mid-year values during the 1990s. Similarly, the FPC1 for 1995–2008 corresponds to more extreme mid-year values during the second half of the 2000s, as well as very extreme values at the end of 2000. Fluctuations around the overall mean for FPC1 mirror the high variability of phosphorus concentration and of *Bythotrephes* density during 1981–1993.

Dynamics of the other principal components (FPC2-4) are highlighted in Figure 2b,c. Briefly, FPC2 corresponds to more extreme values during the 1980s (especially in 1981, 1982, 1983 and 1986) and during the 2000s. Strong deviations observed during the 1980s can be attributed to fluctuations of two parameters with the highest scores on this component: The density of *Bosmina* and the concentration of chlorophyll-α. FPC3 gives more weight to 1982 and 1984, as well as 1996, 1997, 2003, and 2005. The high variability, in some of these years, of chlorophyll concentration and *Daphnia* density explains much of the observed fluctuations around the overall mean. Lastly, FPC4 gives more weight to 1985 and the 2000s, especially 2005. The observed perturbations from the mean during the 2000s can be explained by variation in *Leptodora* population density.

**Figure 2.** *Cont.*

**Figure 2.** Temporal fluctuations from the mean of all standardized functional variables (solid line) obtained by adding (dashed line) or subtracting (dotted line) the extracted FPCs, each scaled by 1.5 times the square root of its explained variance (panels: **a**, FPC1; **b**, FPC2; **c**, FPC3; and **d**, FPC4).

#### *3.4. Functional Principal Components Scores*

The scores for the FPC1 show how, during 1981–1992, *Bythotrephes* population density (score 48.2) increased dramatically during the 1990s, while phosphorus concentration decreased over time (score −37.8) (Table 1). Phosphorus concentration maintained a similar trend during 1995–2008 (score −45.2), while water temperature increased during both considered periods (score of 30.5 in 1981–1992 and of 30.2 in 1995–2008), fitting FPC1 better than *Bythotrephes* density (score 19.8) (Table 1).


**Table 1.** Functional Principal Components (FPCs) scores for three different groups (in italics) of considered variables. The most extreme positive and negative scores are highlighted (italics) for each FPC. Pop. dens. = population density.

Scores for FPC2 during 1981–1992 show how this principal component had the most effect on *Bosmina* density (score 37.5), which reached higher values at the start of 1980 (Figure 2c). By contrast, rotifer density (score −21.1) was stable at the start of the 1980s and more variable at the end of the decade (Table 1). During 1995-2008, FPC2 was affected in a similar way by both *Bosmina* density (score 46.1) and rotifer density (score −21.8). Chlorophyll-α concentration (score −29.1), which reaches its lowest values in 2001 (Figure 2b), also becomes important for FPC2 in this period.

The scores for FPC3 during 1981–1992 show that this FPC were most affected by cyclopoid copepod density (score 41.8), which reached high values in both 1982 and 1984 (Figure 2c), and phosphorus (score −25.9) which reached a minimum concentration in 1984 as oligotrophication proceeded (Figure 2a). Between 1995–2008, FPC3 was most affected by *Daphnia* density (score 33.8), which reached high values in 2003 and 2005 (Figure 2b), and phytoplankton density (score −39.6), which had very high values in 1996–1997 (Figure 2b).

Lastly, in 1981–1992, FPC4 was most affected by phosphorus concentration (score −46.9), which had a minimum in 1985 (Figure 2a), and *Bosmina* density (score 26.0), which reached one of its highest values in the same year (Figure 2c). During 1995–2008, FPC4 was most affected by diaptomid copepod density (score 24.7), whose within-year dynamic changed substantially during the 2000s (Figure 2c), and *Leptodora* density (score −33.5), which reached high values in the summer 2005 (Figure 2c).

#### **4. Discussion**

We used Functional Data Analysis to describe long-term abiotic and biotic dynamics of Lake Maggiore. FDA has features that made it ideal for long-term studies like ours. FDA provides a useful framework to analyzing short and long-term dynamics of a lake ecosystem, its data smoothing respects the continuous nature of the underlying biological and physico-chemical processes while accounting for potential irregularities during measurement periods, and it allows analysis of dominant modes of variation in trajectories of measured parameters in a way that respects the time ordering of observations.

FPCA proved useful in identifying Lake Maggiore's responses to both trophic and climatic variability. Among the extracted FPCs, FPC1 illustrated the increasing trend of *Bythotrephes* population density and decreasing trend of phosphorus concentration at the beginning of oligotrophication, between 1981–1992. In 1996 and 1997, large blooms of ultraplankton—including small cyanobacteria—occurred in the lake, coincident with very low mean densities of *Daphnia*. However, as lake productivity declined owing to phosphorus abatement [25–27,42] below 15 mg m<sup>−</sup>3, the level at which TP mix becomes limiting for phytoplankton growth [27], *Bythotrephes* increased in density, while *Daphnia* declined [47,48]. A previous study has demonstrated that *Bythotrephes* occurs primarily in European lakes with low summer chlorophyll concentration [49], consistent with its resurgence during oligotrophication in Lake Maggiore. The decline in *Daphnia* (see Table 1) could be due to either reduced food supply (bottom-up response) or to increased predation by *Bythotrephes* (top-down response), or both. Available evidence suggests that food reduction may have caused the decline between 1983 and 1987, as declining phytoplankton stocks may have caused food limitation. However, beyond this point, predation by *Bythotrephes* likely accelerated the decline in *Daphnia* density, as the latter's death rate increased when *Bythotrephes* populations surged beyond 20 ind·m−<sup>3</sup> beginning in 1987 [47]. Larger and seemingly more potent *Bythotrephes* sharply curtail crustacean abundance when present at levels beyond ~2 ind·m−<sup>3</sup> in Ontario lakes [50]. Interestingly, *Leptodora* and *Bythotrephes* have inverse distributions in Ontario lakes with the former often replaced by the latter after it invades new systems [51]. These species show no such pattern in Lake Maggiore, where both species are native and presumably co-adapted to one another's presence. In addition, body size differences between the species are much smaller than in Ontario lakes, suggesting less physical dominance by *Bythotrephes*.

In addition to dramatic changes in population density, *Bythotrephes* has also experienced large shifts in seasonality, with peak density shifting from August to May (see Figure 2c) as the lake warmed [31,36]. This, in turn, may have increased predation pressure on *Daphnia* [47,48], which previously typically achieved population maxima in July but now peaks months earlier [52]. The second phase of *Daphnia* increase, during full and stable oligotrophy, resulted both from *Daphnia's* ability to cope with increased invertebrate predation by means of changes in population phenology and by a release from fish predation consequent to a decrease in coregonid (i.e., whitefish) abundance, which also favors consumption of *Bythotrephes* over *Daphnia* [53,54]. Previous research has demonstrated strong preference of coregonid fishes for *Bythotrephes* prey [55].

Oligotrophication was associated with a decline in phytoplankton density between 1995–2008, as clearly shown by both FPC2 and FPC3. However, this reduction in plankton food resources does not seem to have affected most cladocerans as FPC2 highlighted an increase of *Bosmina* whereas FPC3 captured the same trend for *Daphnia*. We hypothesize that the long-term increase in water temperature stimulated the growth of *Daphnia*—as demonstrated by peak population density achieved during the heat wave of summer 2003—possibly owing to reduced development time [56,57].

FPC analyses are compatible with a response by planktonic communities to both temperature increase and total phosphorus concentration decline. The former variable enhances primary production with its consequent indirect, positive effect on herbivores and negative effect on *Bythotrephes*, whereas the latter appeared to favor *Bythotrephes* and its strong top-down effects.

Raptorial dominance declined after 1988 concomitant with changes in cladocerans: *Bythotrephes longimanus* increased by an order of magnitude, *Daphnia* gr. *hyalina-galeata* showed a sharp decline and small cladocerans such as *Eubosmina longispina* and *Diaphanosoma brachyurum* increased [26,29]. An increase in *B. longimanus* releases rotifers, especially large *Conochilus hippocrepis* (or *C. volvox ehrbg*, Schrank, 1803) colonies, from competition with cladocerans [58]. The decrease in abundance and mean size of phytoplankton cells with lake restoration might have enhanced microphagous species able to ingest multiple food items of ~15–20 μm, while raptorials prefer larger particles [46,59–61]. Decreased abundance of competitors and decreased food size, therefore, might enhance microphagous species, while raptorial species were able to coexist with their competitors [46,59,62]. Large (≥1 mm diameter) *Conochilus* colonies have feeding rates comparable to *Daphnia* and lead to a seasonal increase in water clarity like that caused by *Daphnia* grazing [58].

The response of rotifers to changes in the macrozooplankton assemblage is probably more complex than previously hypothesized. Manca [63] reported a strong increase of *Conochilus hippocrepis* in Lake Maggiore when *Bythotrephes* became increasingly abundant in the late 1980s. However, the FPCA highlighted that the increase in rotifer density preceded that of *Bythotrephes* and the two populations seemed to follow divergent development trajectories during the periods examined. Therefore, the behavior of rotifers during oligotrophication of Lake Maggiore is different from that observed in Ontario lakes, where they benefited numerically from competitive release associated with reduced density of *Daphnia* competitors when *Bythotrephes* was abundant [64]. Once established, *C. hippocrepis* populations would be relatively immune to growing *Bythotrephes* populations owing to its gelatinous colonial matrix.

Exceptional meteoclimatic events can be detected over the long term [65]. As an example, the almost complete vertical mixing of February 1999 and the record flood of 2000 were also reported to impact zooplankton seasonal dynamics. As reported by Manca, Cavicchioni and Morabito [66], nutrient replenishment of surface waters from deeper waters consequent to vertical mixing triggered a cascading effect, involving all planktonic communities. Larger input of new nutrients stimulated phytoplankton productivity and higher food availability for herbivorous zooplankton, resulting in an increased density of rotifers, copepods, and cladocerans, as well as of *Daphnia* clutch size. Finally, an increase of predatory cladocerans (*Leptodora*) was recorded in July, five months after the complete mixing episode. Temperature can control nutrient supply, by affecting the extent of the mixing layer at the time of spring overturn: This mechanism is particularly effective in deep Italian lakes [67,68], where a strong link between winter temperature, mixing depth and was established. Similarly, the record flood of 2000 is well identified in the time series analyzed. Impact on zooplankton included detection in the pelagic of littoral taxa such as *Bosmina longirostris* and altered the seasonality of *Bythotrephes longimanus* [52].

#### **5. Conclusions**

Long-term data allow for disentangling the contribution of different stressors to changes in response variables. In Lake Maggiore, many stressors altered the system during the 20th century, prominent among them cultural eutrophication in the 1960s and 1970s. Phosphorus abatement began in 1977, with strong observed TP reductions observed thereafter. This change had pronounced effects

on numerous biological response variables, notably chlorophyll-α concentration and phytoplankton community size composition, and *Daphnia* population density. We also observed increased variability and enhanced population density of the invertebrate predator *Bythotrephes longimanus* during the initial study period, and by small-bodied cladocerans *Bosmina* and *Diaphanosoma* during the second period. The increase in population density and inter-annual variability of small cladocerans might be viewed as a response to decreasing *Daphnia* population density. *Bythotrephes* population densities increased during the second study period, resulting in enhanced predation pressure on *Daphnia* prey. Impact on *Daphnia* during this period was reduced by a temperature increase - which served to reduce its development time and increase its birth rate - and by earlier phenology such that it exhibited less temporal overlap with *Bythotrephes.* As previously observed, over the long term, both cladocera and rotifers responded more clearly to eutrophication reversal and to warming, than copepods. The latter seem to be more influenced by exceptional meteoclimatic events, such as mixing depth. This result agrees well with analyses of ecological roles of Lake Maggiore freshwater zooplankton by carbon and nitrogen stable isotope analyses. Copepods, particularly cyclopoids, occupy a distinct functional group from the other zooplankton secondary consumers, which is an observation also made in previous studies [54,69,70].

**Author Contributions:** Authors contributed as follows: conceptualization, M.M., H.J.M., A.Z.; methodology, A.A., A.Z., P.Q., M.M.; software, A.A.; validation, A.A.; formal analysis, A.A., A.Z., P.Q.; resources, M.M.; investigation, A.Z., P.Q., H.J.M.; resources, M.M.; Writing-Original Draft Preparation, A.A., A.Z., P.Q., H.J.M.; Writing-Review & Editing, A.A., A.Z., P.Q., H.J.M., M.M.; visualization, A.A., supervision, A.Z., P.Q., M.M.; project administration, A.Z., A.A.; funding acquisition, M.M.

**Funding:** This research was funded to M.M. by the "Long-term limnological research in Lago Maggiore of the International Commission for the Protection of Swiss-Italian Waters" (CIPAIS), in the frame of a cooperation agreement between Swiss and Italian Governments. H.J.M. was supported by the CNR Short Term Mobility Program year 2013, by a Canada Research Chair. A.A. was supported by a CNR cooperation contract on "Management of biological time series with advanced statistical methods.

**Acknowledgments:** We dedicate this paper in memory of Giuseppe Morabito, who studied and loved Lake Maggiore. We are grateful to three anonymous reviewers for very helpful comments.

**Conflicts of Interest:** The authors have no conflicts of interest to declare.

#### **Appendix A. Technical Details**

In FDA, data is assumed to be represented by a collection of smooth curves *x*1(*t*), ... , *xm*(*t*) defined on some interval [*a*, *b*], possibly obtained by smoothing some time series as described in the main text. These curves are then thought as single points *x*1, ... , *xm* belonging to the Hilbert space *L*2[*a*, *b*] of square-integrable functions.

In order to clarify how FPCs provide the best representation of the variation in the data *x*1, ... , *xm* in a least-squares sense, we introduce the following definitions.

First, the inner product ·, · and norm · in the space *<sup>L</sup>*2[*a*, *<sup>b</sup>*] are, respectively, defined as

$$\langle f, g \rangle = \int\_{a}^{b} f(t)g(t)dt \text{ and } \|f\| = \sqrt{\langle f, f \rangle} \tag{A1}$$

for all functions *<sup>f</sup>*, *<sup>g</sup>* <sup>∈</sup> *<sup>L</sup>*2[*a*, *<sup>b</sup>*].

Second, a set of functions *f*1, ... , *fn* is defined to be *orthonormal* if *fi*, *fi* <sup>=</sup> 1 and  *fi*, *fj* = 0 whenever *i j* for all *i* and *j* between 1 and *n*.

Third, the functions *x*1, ... , *xm* in *<sup>L</sup>*2[*a*, *<sup>b</sup>*] are defined as *xi* <sup>=</sup> *xi* <sup>−</sup> *<sup>x</sup>* for each *<sup>i</sup>* <sup>=</sup> 1, ... , *<sup>m</sup>* (where *x* = *m*−<sup>1</sup> *<sup>m</sup> <sup>i</sup>*=<sup>1</sup> *xi*), so that *x*1, ... , *xm* represent the centered data.

*Water* **2019**, *11*, 962

In FPCA, we seek orthonormal elements *e*1, ... ,*en* of *L*2[*a*, *b*] and scores *fij* such that the approximations

$$\widetilde{\mathbf{x}\_i} \cong \sum\_{j=1}^n f\_{ij} \mathbf{e}\_j \tag{A2}$$

yield the lowest possible *Sum-of-Squares Error (SSE)*

$$SSE = \sum\_{i=1}^{m} \|\widetilde{\mathbf{x}\_i} - \sum\_{j=1}^{n} f\_{i\bar{j}} x\_j\|^2. \tag{A3}$$

Chosen any specific orthonormal set *e*1, ... , *en*, the *SSE* can be written as

$$\begin{split} \sum\_{i=1}^{m} \|\widetilde{\mathbf{x}\_{i}}\|^{2} - 2\sum\_{i=1}^{m} \sum\_{j=1}^{n} f\_{\overline{ij}} \langle \widetilde{\mathbf{x}\_{i}}, \mathbf{c}\_{j} \rangle + \sum\_{i=1}^{m} \sum\_{j=1}^{n} f\_{\overline{ij}}^{2} \\ &= \sum\_{i=1}^{m} \|\widetilde{\mathbf{x}\_{i}}\|^{2} - \sum\_{i=1}^{m} \sum\_{j=1}^{n} \langle \widetilde{\mathbf{x}\_{i}}, \mathbf{c}\_{j} \rangle^{2} + \sum\_{i=1}^{m} \sum\_{j=1}^{n} \left( \langle \widetilde{\mathbf{x}\_{i}}, \mathbf{c}\_{j} \rangle - f\_{\overline{ij}} \right)^{2} \end{split} \tag{A4}$$

which achieves its minimum when *fij* <sup>=</sup> *xi*,*ej*.

With this choice of coefficients and since *fij* = *xi* − *x*,*ej* = *xi*,*ej*−*x*,*ej*, the *SSE* becomes

$$\sum\_{i=1}^{m} \left\| \mathbf{x}\_{i} - \overline{\mathbf{x}} \right\|^{2} - \sum\_{i=1}^{m} \sum\_{j=1}^{n} \left( \langle \mathbf{x}\_{i}, e\_{j} \rangle - \langle \overline{\mathbf{x}}, e\_{j} \rangle \right)^{2},\tag{A5}$$

which attains its minimum value when each quantity σ<sup>2</sup> *<sup>j</sup>* <sup>=</sup> <sup>1</sup> *n n j*=1 *xi*,*ej*−*x*,*ej* 2 , corresponding to the variance of the scores *f*1*j*, ... , *fmj*, is maximized. Since FPCs are selected in order to maximize these variances, FPCs provide the minimum *SSE* among all possible orthonormal set *e*1, ... , *en*.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Cladoceran (Crustacea) Niches, Sex, and Sun Bathing—A Long-Term Record of Tundra Lake (Lapland) Functioning and Paleo-Optics**

#### **Liisa Nevalainen 1,\*, E. Henriikka Kivilä 1, Marttiina V. Rantala 1,2 and Tomi P. Luoto <sup>1</sup>**


Received: 2 September 2019; Accepted: 24 September 2019; Published: 27 September 2019

**Abstract:** Under fundamental ecosystem changes in high latitude lakes, a functional paleolimnological approach may increase holistic understanding of lake responses and resilience to climate warming. A ~2000-year sediment record from Lake Loažžejávri in the tundra of northern Finnish Lapland was examined for fossil Cladocera assemblages to examine long-term environmental controls on aquatic communities. In addition, cladoceran functional attributes, including functional diversity (FD), UV absorbance (ABSUV) of *Alona* carapaces, and sexual reproduction (ephippia) in *Bosmina* and Chydoridae were analyzed. Cladoceran communities responded to a major change in benthic habitat quality, reflected as elevated (increasingly benthic) sediment organic matter δ13C signal since the 17th century. FD fluctuations showed association with climate oscillation, FD being generally higher during warm climate periods. These ecological changes were likely attributable to diversification of littoral-benthic consumer habitat space. ABSUV, irrespective of increases during the Little Ice Age (LIA) due to higher UV transparency of lake water, was lower under increasing autochthony (benthic production) suggesting establishment of physical UV refugia by the benthic vegetative substrata. *Bosmina* ephippia exhibited a decreasing trend associated with increasing benthic production, indicating favorable environmental regime, and, together with chydorid ephippia, transient increases during the climate cooling of the LIA driven by shorter open-water season.

**Keywords:** autochthony; cladocera; functional ecology; organic carbon; paleolimnology; tundra lakes; UV radiation

#### **1. Introduction**

Impacts of recent climate warming are emphasized in small and shallow high latitude lakes being driven by higher air and water temperature and longer open-water season [1]. In a long-term perspective, these major physical drivers shift arctic and subarctic aquatic ecosystems toward an unprecedented ecological status [2,3]. Identification and characterization of the new high latitude lake trajectories is significant since they connect lake food webs, by aquatic-terrestrial coupling, to global-scale biogeochemical processes [4]. Northern ecotonal tree line lakes and their sedimentary environmental archives act as sentinels for estimating current ecosystem functioning and organization with respect to natural variability over the course of the Holocene [5].

The balance between autochthonous (in-lake produced) and allochthonous (catchment-originated) organic matter in aquatic systems is of high importance to the global carbon cycle, since lakes store and transfer carbon and act across the atmospheric-terrestrial-aquatic boundaries [4]. A common character of high latitude lakes is high water transparency, i.e., low amount of terrestrially derived

dissolved organic matter (usually estimated as dissolved organic carbon (DOC)) due to lack or scarcity of catchment vegetation [6]. While typically systems with DOC concentrations over 5 μg L–1 are estimated to be heterotrophic [7], the high latitude lake food webs with low (approximately 1–3 μg L–1) DOC are principally dependent on autochthonous production. High water transparency induces low attenuation of sunlight and UV radiation and organisms must therefore cope with high UV exposure [8]. Furthermore, food webs and secondary production may be supported primarily by benthic primary production that is fueled by the ample light [9]. However, under climate warming, i.e., due to the advancing tree line, expanding tundra vegetation or thawing permafrost, previously autochthonous and transparent lakes may experience increases in input of allochthonous organic matter from the catchment causing fundamental changes in light climate and energy pathways [5,10,11]. Alternatively, climate warming may enhance autochthony through lengthening of the open-water season and increasing habitat availability for benthic primary producers [12,13].

Knowledge of ecological and biogeochemical functions of aquatic organisms, e.g., feeding guilds, habitat preferences, and photoprotective mechanisms, is important when examining lake food webs and the organic matter cycle. Paleolimnological research focusing on long-term distribution of microscopic aquatic fossils (e.g., crustacean cladocerans or diatom algae) has slowly shifted toward biodiversity sciences [14–16] but an actual functional approach, where functional rather than taxonomic classification of fauna and flora is assessed in relation to natural and anthropogenic ecosystem variability, still remains rather rare. The functional approach may allow a holistic understanding of changes and drivers occurring in lakes and their surroundings since ecosystem functions are not dependent on taxonomic identity [17,18]. In addition to the use of functional diversity (FD) as an index for assessing functional trait distribution in fossil biological assemblages, functionality may be estimated, for example, by patterns in feeding, habitats, reproduction, and morphology [19–22].

In the current study, our focus was on Cladocera (Crustacea) communities and their functioning over the past millennia in a subarctic tundra lake. More specifically, we analyzed fossil cladoceran communities and functional attributes including functional diversity (FD index), reproduction patterns (Chydoridae and *Bosmina* fossil ephippia), and photoprotection (melanization based on carapace UV absorbance) in a 2000-yr sediment core and compared these data to previously available paleoclimate and biogeochemical proxies. The objectives were to track long-term changes in cladoceran community functioning and lake water bio-optics, and to identify interconnections of functional ecology and past climate-driven limnological changes.

#### **2. Materials and Methods**

Lake Loažžejávri is located in the subarctic (mean July temperature +12.3 ◦C, mean annual temperature −2 ◦C) shrub tundra of Finnish Lapland (69◦53 N, 26◦55 E) at an altitude of 255 m a.s.l. (Figure 1). The lake has an area of ~3 ha with a catchment size of ~200 ha and it is shallow (water depth 1.2 m), oligotrophic (total phosphorus (TP) 5.9 μg L<sup>−</sup>1, chlorophyll-α (chl-*a*) 1.6 μg L<sup>−</sup>1) and transparent (dissolved organic carbon (DOC) 3.4 mg L<sup>−</sup>1, UV attenuation coefficient at 305 nm (*K*dUV) 11.1 m<sup>−</sup>1). A 38-cm sediment sequence was cored from the basin in late summer 2014 with a Limnos gravity corer and subsampled at 1-cm intervals and dated with the 14C dating method. The core covers approximately the past two millennia. The lake has previously been investigated as a part of regional limnological survey [23], and its sediments for late Holocene paleoclimate [24] and carbon utilization of aquatic macroinvertebrate communities [25], from where detailed information on limnological characteristics, catchment and sedimentology, chronology, and biostratigraphy can be found. The age-depth model for the core was based on two plant macrofossils at sediment depths 20 (1390–1440 C.E.) and 25 cm (970–1025 C.E.) analyzed for AMS (accelerator mass spectrometry) 14C dates. Additionally, 137Cs and 210Pb were analyzed but their concentrations remained low and therefore insufficient to add further temporal resolution to the age-depth model [24].

**Figure 1.** Location of Lake Loažžejávri in northern Finnish Lapland (**red dot in left column**) and its catchment characteristics (**right column**).

Fossil Cladocera were analyzed from the sediment subsamples (1-cm resolution) to examine past community changes. The samples were prepared by heating and stirring ~2 cm3 of fresh sediment in 10% potassium hydroxide (KOH on a hot plate following the standard protocol [26]. The samples were then sieved through a 51-μm mesh and the residues were centrifuged for 10 min at 4000 rpm. Permanent microscope samples were mounted in glycerol gelatin stained with safranine on a hot plate. The micropscope slides were examined for fossil cladoceran remains (carapaces, head shields, postabdomens, claws, ephippia) with a light microscope at 100–400× magnifications and a minimum of 70 individuals were identified [27] from each sample (53 individuals from sample at 3 cm due to low amount of remains).

In addition to standard community analysis described above, functional attributes of sexual reproduction, functional diversity, and melanization were analyzed. Total bosminid (Bosminidae, *Bosmina*) and chydorid (Chydoridae) ephippia (chitinous envelops for diapausal resting eggs, indicative of sexual reproduction) were enumerated during the standard counting and were used for estimates of sexual reproduction (ratio of sexual to asexual reproduction) by enumerating *Bosmina* and chydorid ephippia (indicative of sexual reproduction) and carapaces (indicative of asexual reproduction) [28]. Functional diversity of the cladoceran community was expressed as Rao's FD index, i.e., Rao's quadratic entropy [29]. For this index, each fossil cladoceran taxa was assigned with selected qualitative functional characters (traits, Table 1) 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 habitat (pelagic, benthic, attached to vegetation) following [21]. To examine melanin in cladoceran remains as an index for UV exposure and photoprotection [30], fossil cladoceran carapaces from a large chydorid species *Alona a*ffi*nis* (Leydig) were extracted from sediment subsamples with fine forceps under a binocular microscope. The carapaces were then measured under UV wavelengths 340 and 305 nm for their UV absorbance (ABSUV) to indicate the degree of melanization. The absorbance measurements were performed with a UV/VIS spectrophotometer (UV-1800, Shimadzu Corporation, Kyoto, Japan). Seven carapaces were measured for a mean UV absorbance value per sample omitting highest and lowest absorbance values. UV absorbance measurements are expressed as anomalies from the mean absorbance value of the data series.

The data were further analyzed with numerical methods to elucidate relationships between cladoceran community and functional attributes and paleoenvironmental variation. Redundancy analysis (RDA) was utilized to examine impacts of sediment biogeochemistry (indicative of habitat quality), i.e., elemental and isotopic composition (C%, N%, C/N ratio, δ13C, and δ15N) and amount of organic matter (OM%, data from [25]), on cladoceran community succession though time. In addition, external atmospheric forcers, i.e., sun spot numbers (SSN) and summer (June–August) air temperature (T-JJA) as 10-year averages [following 31, 32, respectively] were used in the RDA as environmental factors. SSN reconstruction is based on dendrochronologically dated radiocarbon concentrations and physics-based models [31] and T-JJA reconstruction is based on maximum latewood density tree-ring chronologies [32]. Species data were square-root transformed prior to data analyses. Paleoenvironmental parameters with high inflation factors were omitted to include a set of variables with inflation factors <10 (C/N, δ13C, δ15N, OM%, SSN, T-JJA). Forward selection of environmental parameters was performed with 999 permutations and significant parameters were assigned as *p* < 0.05. Selected bivariate environmental correlations for functional attributes were examined with Pearson's correlation coefficient and linear regression. RDA was performed with Canoco 5 software [33] and linear regressions with PAST3 software [34]. In addition to OM%, organic matter δ13C, and SSN, mean carbon source contributions in the benthic food web (based on fossil Chironomidae δ13C), modeled with a three-source (benthic, planktonic, and terrestrial) Bayesian mixing model previously applied in [22], were utilized as paleoenvironmental reference data. Here, the modeled planktonic and benthic carbon contributions (to Chironomidae diet) were used in a planktonic to benthic (P/B) ratio. Further, autochthonous (planktonic + benthic) and allochthonous (terrestrial) carbon sources were used in autochthonous to allochthonous ratio. These ratios were calculated based on the Chironomidae δ13C based source data [25].

**Table 1.** Functional traits of cladoceran taxa encountered from Lake Loažžejávri sediment core. Characterization is based on body size (S = small, M = intermediate, L = large), body shape (G = globular, O = oval, E = elongated), feeding (F = filterer, S-D = scraper-detritivore, P = predator, including parasitism), and habitat (P = pelagial, S = sediment, V = vegetation).


#### **3. Results**

A total of 27 cladoceran taxa were identified from the sediment subsamples. The most frequent taxa (occurring in every subsample) included *Alonella nana* (Hill's N2 36.5, mean relative abundance 42.5%), *Bosmina longispina* (36.4, 32.5%), *Alona a*ffi*nis* (34.5, 7.9%), *Chydorus sphaericus*-type (24.3, 6.1%), and *Paralona pigra* (25.1, 1.5%). In the cladoceran stratigraphy (Figure 2), *B*. *longispina* and *A*. *a*ffi*nis* dominated with ~20%–60% abundances throughout the sediment sequence, *B*. *longispina* peaking between 1400 and 1500 C.E. and A. nana between 1600 and 1800 C.E. Many littoral-benthic taxa, including *C*. *sphaericus*-type, *Acroperus harpae*, *Alonella excisa*, and *Eurycercus* spp. increased slightly during 600–1200 C.E. and others, including *Alona quadrangularis* and *A*. *guttata*, increased or emerged at the top sequence after 1600 C.E.

**Figure 2.** Relative abundance of most common (Hill's N2 > 3.5) cladoceran taxa in the Loažžejávri sediment core indicated with a gray silhouette, where black horizontal lines represent sediment sub-samples. Five-fold exaggeration curves of less abundant taxa are shown with gray lines. Temporal extensions of the cold Little Ice Age (LIA) and warm Medieval Climate Anomaly (MCA) are approximated.

RDA and forward selection of environmental variables resulted in eigenvalues of 0.14 and 0.07 for axes 1 and 2, respectively, explaining 35.2% of variation in species data (Figure 3). It identified δ13C (41.7%, F = 6.2, *p* < 0.001), OM% (20.1%, F = 3.2, *p* < 0.001), and SSN (11.5%, F = 1.9, *p* = 0.035) as the most significant environmental parameters explaining variation in cladoceran species data, together accounting for 73.3% of the explained variation. Most recent samples (12–0 cm, ~1650 C.E.—present) had positive axis 1 scores in relation to increasing δ13C, OM% and SSN. *Alonella nana*, *Alona werestschagini*, *Alona quadrangularis*, and *Alona guttata* were associated with increasing δ13C. *Rhynchotalona falcata*, *Drepanothrix dentata*, *Alona intermedia*, and *Daphnia* spp. were associated with increasing OM% and SSN having more positive scores for RDA axis 2. Older than ~1650 C.E. samples had negative axis 1 scores and the abundant taxa, such as *Bosmina longispina* and *Chydorus sphaericus*-type, and vegetation-associated taxa *Alonella excisa*, *Acroperus harpae*, and *Eurycercus* spp. were related to the negative end of RDA axis 1 (i.e., low δ13C).

Of the cladoceran-based functional indices (Figure 4), chydorid and *Bosmina* sexual reproduction both increased around 700–800 C.E. and peaked later on during 1400, 1800, and 1900 (chydorids) and 1300 (*Bosmina*) C.E., after which *Bosmina* sexual reproduction declined. Minimum values occurred around 1700 and 1900 C.E. in chydorids and *Bosmina*, respectively. Mean abundance of chydorid and *Bosmina* ephippia encountered in the samples was high, 33 (min. 6, max. 65) and 18 (min. 1, max. 50) ephippia, respectively. Cladoceran FD was highly variable, but exhibited increases during 600–1100 C.E. and after 1700 C.E. until the present. Mean ABSUV for the sediment core was 1.2 AU

(absorbance units). ABSUV exhibited several peaks, indicative of more intensive melanization (i.e., higher UV exposure) during the early core at ~300, 600, and 1000 C.E. and in the younger layers ~1500 and 1700–1800 C.E.

**Figure 3.** Redundancy analysis ordination diagram for the cladoceran taxa and samples in relation to significant environmental variables δ13C of sediment organic matter (benthic production), organic matter content (as loss-on-ignition), and solar activity (as reconstructed sun spot numbers). Pre ~1650 C.E. samples (13–37 cm) are indicated by white dots and the recent samples (0–12 cm) by black dots. Data for sediment geochemistry are from [25] and for solar activity from [30]. Species abbreviations include four letters of the genus and three from the species name.

**Figure 4.** Cladoceran based functional attributes; functional diversity, UV absorbance of *Alona* carapaces, ratio of ephippia (sexual reproduction) to carapaces (asexual reproduction) of Chydoridae and *Bosmina*, and cladoceran planktonic to littoral ratio. In addition, ratios of mean source contributions (planktonic to benthic and autochthonous to allochthonous) in the benthic food web (based on Chironomidae δ13C and Bayesian mixing modeling), sediment organic matter δ13C, organic matter content (as loss-on-ignition), and solar intensity (as reconstructed sun spot numbers) are shown. Data for sediment geochemistry and source contributions are after [25] and for solar intensity after [31]. Temporal extensions of the cold Little Ice Age (LIA) and warm Medieval Climate Anomaly (MCA) are approximated.

Significant Pearson's correlations were found between T-JJA and FD (r = 0.48, *p* = 0.002; Figure 5a), sediment C/N ratio and carapace UV absorbance (r = 0.46, *p* = 0.004; Figure 5b), and sediment carbon content and *Bosmina* ephippia/carapaces ratio (r = −0.60, *p* = < 0.001; Figure 5c). In addition, ABSUV had significant (*p* = <0.05) correlations with C% (r = −0.37) and N% (r = −0.33) and *Bosmina* ephippia with N% (r <sup>=</sup> <sup>−</sup>0.56), C/N (r <sup>=</sup> 0.32), <sup>δ</sup>13C (r <sup>=</sup> <sup>−</sup>0.59), <sup>δ</sup>15N (r <sup>=</sup> 0.64), and OM% (r <sup>=</sup> <sup>−</sup>0.59).

**Figure 5.** Bivariate relationships (*r* = Pearson's correlation, *r*<sup>2</sup> = coefficient of determination in linear regressions, *p* = significance) of (**a**) cladoceran functional diversity and reconstructed summer (June–August) mean temperature, (**b**) carapace UV absorbance anomaly and sediment carbon to nitrogen ratio (C/N), and (**c**) ratio of ephippia (sexual reproduction) to carapaces (asexual reproduction) in *Bosmina* and sediment carbon content. Black dots represent sediment core subsamples, black dashed lines indicate linear regression models and gray lines 95% confidence intervals. Sediment geochemical data originates from [25] and summer mean temperature reconstruction from [32].

#### **4. Discussion**

#### *4.1. Lake Functioning*

While recent research has highlighted the role of northern waters in the cycling of terrestrial organic matter, there also exists contrasting evidence suggesting insignificant contribution of terrestrial inputs to aquatic systems [35–37]. For example, a trend of increasing autochthony rather than allochthony has been observed in small lakes from northern Finnish Lapland [13,25]. Benthic autochthonous production prevailed in Lake Loažžejávri over the past two millennia supporting aquatic consumers such as macrobenthic Chironomidae (chironomid) larvae [24,25]. Of the autochthonous organic material consumed, planktonic carbon component, inferred and modeled from δ13C of chironomid head capsules, was slightly elevated (increase in P/B, Figure 4) during the warm Medieval Climate Anomaly (MCA) at around 900–1300 C.E. This increase is likely indicative of higher algal production and its more efficient utilization by consumers [25]. Benthic production and benthic carbon component in chironomid diet (inferred from δ13C of chironomid head capsules) clearly dominated in the lake after ~1600 C.E. and may have been a reflection of reduced ice-cover period after the LIA climax creating wider benthic niche space for autotrophic organisms generally increasing the relative importance of the benthic habitat over the pelagic one [38].

The production in the lake being mostly benthic, the main driver for long-term cladoceran community shifts in Lake Loažžejávri sediment profile was related to the benthic habitat quality. According to the RDA, δ13C signature indicating the origin of organic matter, i.e., autochthonous (planktonic or benthic) or allochthonous, was the most significant environmental parameter explaining cladoceran assemblages (Figure 3). The isotopic carbon signature separated the sediment core subsamples along the RDA axis 1 to pre- and post-1650 C.E. sample clusters (Figure 3). Organic matter δ13C varied relatively little during the first 1500 years of the sequence but started to increase at ~1600 C.E. indicative of more pronounced benthic primary production [25]. In concert, the cladoceran community diversified with the inclusion and increase of several benthic taxa, e.g., *Alona werestschagini*, *A*. *quadrangularis*, *A*. *guttata*, *R*. *falcata*, and *Drepanothrix dentata* until the 21st century (Figure 2). Chydorid

benthos has been shown to utilize unselectively periphyton and detritus for food in similar lake environments in the region [13] and the assemblage change was thus likely driven by a diversification of microhabitats among the benthic vegetative substrata [39]. Alternatively, the diversification of the benthic habitats may have allowed species with more specialized feeding strategies to increase in the lake following the benthic succession. In common, consistent increases of specialized cladoceran taxa in relation to climate warming have been reported from arctic Canada [40].

Despite its shallowness, the lake supported a planktonic food web, as evidenced by the high abundance of euplanktonic *Bosmina longispina*, presence of *Ceriodaphnia*, and scattered occurrence of *Daphnia* throughout the core (Figure 2) as well as presence of some free living Chironomidae (*Ablabesmya* and *Procladius*, [24]). In the current results, *Bosmina* was associated with more negative δ13C values of organic matter (Figure 3), indicative of an association with phytoplankton production. *Bosmina* is known to feed selectively on phytoplankton in small and shallow lakes in northern Finland [13,41]. In contrast to *Bosmina*, some planktonic cladocerans, e.g., *Daphnia* and *Ceriodaphnia* may assimilate benthic food particles grazing directly over the sediment microbial active layer [9]. In the RDA, *Daphnia* was associated with elevated δ13C of organic matter suggesting a connection with the benthic production. Related to the shift in autochthonous production (i.e., increase in benthic production), *Bosmina* decreased slightly but still remained as the dominant planktonic taxon, while *Daphnia* increased at the top of the sediment core (Figure 2). Accordingly, it is possible that the establishment and increase of *Daphnia* ~1850 C.E. onward was related to its ability to utilize benthic resources [42]. In all, planktonic cladocerans, *B. longispina* as the governing taxon, showed rather constant presence in the sediment sequence aside from a transient increase at ~1400 C.E. (Figure 2) This abrupt event was likely associated with a short-term (~100 years) hydroclimatic event of the early Little Ice Age (LIA) causing higher lake level and relative increase of the planktonic habitat.

Despite the shift toward benthic production, *Bosmina* showed high resilience [41] as it remained abundant up to the surface (Figure 2). Though, a change in its reproduction occurred with a decreasing trend in sexual reproduction after a peak around 1300 C.E. suggesting less environmental stress (lower sexual reproduction) during the recent centuries. Sexual reproduction in Cladocera is evidenced in the fossil record through their ephippia and has previously been investigated in relation to climate oscillation, i.e., open-water season length, which dictates the relative importance of asexual (parthenogenetic) vs. sexual (gamogenetic) reproduction [28,43]. Accordingly, it has been suggested that sexual reproduction is less significant (lower sexual reproduction) during warm climate conditions as asexual reproduction prevails under long open-water season. As cold climate of the LIA started to prevail in northern Finland already during the 15th century, the long-term trend in *Bosmina* reproduction was not apparently strictly dictated by open-water season length but also other environmental stimuli counteracted in the reproduction patterns. *Bosmina* sex ratios had a strong negative relationship with sediment carbon content (and organic matter content, Figure 5c) possibly indicative of the past low productive conditions (i.e., low sediment carbon content) being more unfavorable to *Bosmina* than the new, i.e., post-1650 C.E., regime with increasing benthic autotochthonous production. Chydorid sexual reproduction increased during the 15th century, and later peaks ~1800 and 1900, that were likely attributable to reduced length of the open-water season driven by cold climate events of the LIA [44,45].

FD combines the observations of biodiversity and ecological functions to detect trends and patterns in ecosystem functions [17,46]. In the current study, cladoceran FD was generally lower during the early core and between 1200 and 1800 C.E. and higher between ~900 and 1200 C.E. and in the top sequence, from ~1850 onward (Figure 4). The pristine Lake Loažžejávri is mostly governed by natural forcers such as climate due to its remote location in the tundra [25] and, in agreement, there existed a positive relationship between FD and reconstructed T-JJA (Figure 5a) as FD was higher during the MCA and during the post-LIA climate warming (Figure 4). Based on previous paleoecological records [21,47], it has been suggested that long-term development of cladoceran FD is connected with lake productivity, although dependent on the size and eutrophication history of the system in question. We did not find

any statistical relationships between FD and biogeochemical proxies related to productivity but the established relationship between FD and T-JJA likely reflects the ultimate climatic (i.e., temperature) control on habitat development and productivity in this cold tundra environment through dictating the length of ice-cover period [48].

#### *4.2. Paleo-Optics*

Natural variability in light and UV attenuation in forest and tundra lakes may be strongly related to coupling processes with temperature and precipitation via surrounding vegetation and permafrost [5,30]. However, in aquatic systems located on barren catchments terrestrial inputs of UV screening compounds are low and autochthonous organic matter, e.g., algal biomass acts as an important component in UV attenuation [49]. The tundra catchment of Lake Loažžejávri is currently covered with only ~5% of wetlands [23] but the lake shore is directly connected with paludified shorelines and water channels up in the nested drainage basin (Figure 1). Even though Loažžejávri features low lake water DOC (3.4 mg L<sup>−</sup>1) and chl-a (1.6 μg L<sup>−</sup>1), it has been estimated that only a few percent of UV penetrates below 0.5 m water depth [50]. However, the shallow littoral and benthic habitats of the lake are exposed to UV and aquatic biota are thus likely to respond to UV. The ABSUV record of Loažžejávri, indicating melanization and past UV exposure of the benthic *Alona*, was highly variable in the past suggesting large natural variations in underwater UV regimes (Figure 4).

In arctic and subarctic lakes, long-term changes in UV exposure have been mostly governed directly, or indirectly through photodegradation of OM, by solar radiation intensity [51–53]. Although solar forcing interacted partly in cladoceran assemblage succession, SSN reconstruction explaining ~12% of the species variation (Figure 3), it did not have substantial association with the current ABSUV record (Figure 4). In an adjacent lake, which differed from Loažžejávri by its location in the mountain birch woodland at a lower altitude, millennial variation in ABSUV was in clear connection with solar intensity [13,53]. In Loažžejávri, highest carapace ABSUV was recorded in the early core (~300 and 500 C.E.) and later on during 1400–1500 and 1700–1800 C.E. The high variability of ABSUV in the early core does not seem to clearly correlate with changes in SSN reconstruction, sediment geochemistry or cladoceran communities (Figures 2 and 4). However, the latter ABSUV peaks were likely related to a more transparent water column during the LIA, related to constrained organic matter flux from the catchment and reduced in-lake production driven by longer ice-cover period. A similar pattern of increased aquatic UV exposure during the LIA has occurred in lakes of northern and eastern Finland and in the Alps [30,53]. The pre MCA period at the middle of the first millennium of the Common Era has been characterized to be hydroclimatically dry in Fennoscandia. This period overlaps with the Dark Ages Cold Period characterized by noticeable climatic fluctuations [54,55], that may have been driving the early fluctuations in ABSUV (higher peaks ~300 and 500 C.E.) through climate-driven catchment coupling.

There existed a conspicuous long-term relationship between sediment C/N ratio and ABSUV indicating lower UV exposure under increasing autochthony (Figure 5b). This result contradicts a previous investigation from a regional lake set across the subarctic tree line suggesting an intrinsic control of terrestrial DOM (wetland origin) on underwater UV exposure and carapace UV absorbance [50]. As said, Lake Loažžejávri is located on a barren tundra catchment with little terrestrial organic matter contribution throughout its examined history, although connected with some wetland impact [18]. Since ABSUV was lower under prevalently autochthonous conditions, it is possible that UV screening properties of phytoplankton, or specifically those of phytobenthos in this benthic dominated system, have impacted carapace melanization of the benthic *Alona*. In subarctic and arctic lakes, benthic phototrophic communities, including cyanobacteria, algae, and hetero- and chemoautotrophic micro-organisms grow to form thick and stratified biofilms or mats [56], which support benthic and pelagic secondary production and bacterial planktonic production [57,58]. As such, these massive biofilms may act as physical UV refugia for benthic invertebrates allowing them to crawl deeper into the shelter of the mat, since top layers of the benthic mats are abundant in algal pigments

screening UV [59,60]. At the same time as acting as a UV screen, the benthic mats likely provided fertile grazing surface for the microbenthos. It has been suggested previously that UV responses of shallow water benthic communities are principally driven by their access to physical UV refugia [61]. Accordingly, we propose that the generally lower carapace UV absorbance values in the top of the core, with the exception of the two high peaks apparently related to the LIA, were caused by increased benthic autochthonous production providing a physical UV shelter for benthic invertebrates.

#### **5. Conclusions**

Climatic fluctuations of the cold Little Ice Age and after were seen in the studied tundra lake as increases in benthic autochthonous production due to reduction in the length of ice-cover period. This centennial ecosystem scale succession toward dominance of benthic production has altered community structure, induced higher functional diversity, promoted relative importance of asexual reproduction, and reduced UV exposure in cladocerans. Main ecological mechanisms were related to diversification of benthic niche space, likely as a development of benthic microbial mats.

Fossil cladoceran communities in the studied tundra lake, consisting of planktonic and benthic species, seem to have been relatively resilient to climatic fluctuations until their habitat structure was disturbed. Periodically highly abundant sexual reproduction as a way for dormancy or better fitness (genotypic variation) has likely contributed to cladoceran community resilience. The diversification of the benthic niches induced functionally richer cladoceran communities including the keystone planktonic grazer *Daphnia* and specialized benthic species suggesting that functional diversity is coupled with lake productivity. In addition, the establishment of benthic mats likely provided physical UV refugia for benthic cladocerans emphasizing the fundamental importance of habitat quality for these microinvertebrates.

**Author Contributions:** Conceptualization, L.N. and E.H.K.; methodology, L.N., E.H.K., M.V.R. and T.P.L.; software, L.N.; validation, L.N., E.H.K., M.V.R. and T.P.L.; formal analysis, L.N.; investigation, L.N.; resources, L.N.; data curation, L.N.; writing—original draft preparation, L.N.; writing—review and editing, E.H.K., M.V.R. and T.P.L.; visualization, L.N. and E.H.K.; supervision, L.N.; project administration, L.N., E.H.K. and M.V.R.; funding acquisition, L.N.

**Funding:** This research was funded by Academy of Finland, grant number 287547 (VIOLET-project) and 308954 + 314107 (SCUM-project).

**Acknowledgments:** Annukka Galkin and staff of the Kevo Subarctic Research Station are thanked for assistance with the field work. We thank two anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
