Study lakes
The study was conducted in two water bodies created after aquatic restorations of mining pits, Lakes Milada (2.5 km2, high structural complexity; 50° 39′ N, 13° 58′ E) and Most (3.1 km2, low structural complexity lake; 50° 32′ N, 13° 38′ E), in the Czech Republic (Fig. 1). Aquatic restoration took place from 2001 to 2010 in Milada and from 2008 to 2014 in Most. Both lakes are medium-sized (surface area = 252 and 311 hectares, respectively), relatively deep (Zmean = 16 and 22 m, Zmax = 25 and 75 m), oligotrophic (mean summer total phosphorus < 10 and < 5 μg/L) and the Secchi depth varies between 3–10 m in Milada and 7–10 m in Most lake. The deeper Most has a well-oxygenated water column down to 50 m depth, whereas in Milada the profundal zone suffers from poor oxygen conditions below 20 m depth in summer17.
Locations and bathymetric maps of investigated lakes and positions of telemetry systems. Dots represent positions of each telemetry receiver and stars positions of temperature loggers. The map created using the Open Source QGIS version 3.18 (https://qgis.org/en/site/).
Macrophyte sampling was carried out prior to the study in September 2014 and May 201533. Macrophytes were dense only in the HSC lake where they covered 60–91% of 0–12 m deep inshore areas. In the LSC lake, there was only a sparse macrophyte coverage of 0.1–1.6% at 0–3 m depth. Dominant macrophyte species in both lakes were Potamogeton pectinatus, Myriophyllum spicatum and Chara sp.
Both lakes had similar fish community compositions. Roach (Rutilus rutilus) and perch (Perca fluviatilis) were dominant in both lakes, while ruffe (Gymnocephalus cernua), tench (Tinca tinca), European catfish (Silurus glanis), Northern pike and pikeperch (Sander lucioperca) were less abundant34. In addition, pelagic planktivorous maraena whitefish (Coregnus maraena) were present in the LSC lake.
Stocking of predatory species including pike were performed in both lakes for biomanipulation purposes. Stocking of pike was performed after filling of the lakes as a biomanipulation measure to improve and maintain high water quality. In HSC lake, young-of the-year pike were stocked in 2005 and 1–2 year old pike in 2005, while in LSC lake over 1-year old pike were stocked in years 2011–2013 (Detail information in Supplementary material1 Tab. A1). Stocked pike came from various breeding ponds in the Central Bohemia region.
The pike populations were monitored long-term in both lakes and all captured pike were PIT tagged and released back in to the lake (more details are given in35). The population size estimate (based on recaptures) was around 220 individuals (size range 60–120 cm) in HSC lake in years 2014–2016. Low recaptures numbers did not allow making asimilar population size estimate in LSC lake. By comparing efforts and catches in both lakes, we estimated that population in LSC lake accounted for only 30–40% of that in HSC lake population (Vejřík L., unpublished data).
Telemetry system
Two separate MAP positioning systems (Lotek Wireless Inc., Canada) were deployed in the HSC and LSC lakes to track tagged fish (see below). The systems consisted of 91 receivers (Lotek Wireless Inc., WHS3250; 44 receivers in HSC lake, and 47 receivers in the LSC lake) deployed in arrays with distances between the 3 nearest receivers ranging from 203 to 288 m (mean 251 ± 18 m) in the HSC lake and 191 to 341 m (mean 264 ± 33 m) in LSC lake (Fig. 1). The exact position of deployed receivers was measured using a high-precision GNSS-unit (Trimble Geo7x with a cm-precision RTK-service). Depth of receivers varied between 4.5 and 5.5 m. According to range testing done prior to the study, in these lakes (September 2014), such setting of receiver arrays should provide full coverage of both lakes under appropriate environmental conditions. Monitoring of the system accuracy was achieved by 20 stationary reference tags (10 tags in each lake; Lotek Wireless Inc., Canada, model MM-M-16-50-TP, burst rate 25 secs), located in 4 locations in each lake (2 open water locations at depths of 1, 5 and 13 m, and 2 nearshore locations at depths of 1 and 3 m). Further, testing of accuracy was performed by reference tags dragged below a boat after the deployment and before the final retrieval of the telemetry system (further information on positioning error are given in the Supplementary material 2). The systems were installed from April 2015 to March 2016 and the data was manually downloaded every two months. The period targeted in this study lasted from 27 May 2015 to 10 October 2015 to cover the summer period with the highest fish activity and development of macrophytes.
Fish tagging
A total of 30 pike individuals (15 in each lake) were captured by electrofishing (23 individuals), long-lines (6 individuals) and angling (1 individual). Mean total body length/mass was 79 cm/4.13 kg for the HSC lake and 86 cm/4.15 kg for the LSC lake, respectively (more details in Table 1). After capture, pike were anaesthetized by 2-phenoxy-ethanol (SIGMA Chemical Co., USA, 0.7 ml L−1, mean time in anaesthetic bath was 3.75 min), measured, weighed and tagged. A 1–1.5 cm long incision was made on the ventral surface posterior to the pelvic girdle and the transmitter (Lotek Wireless Inc., MM-M-11-28-PM, 65 × 28 mm, mass in air of 13 g, including pressure and motion sensors, burst rate 25 s, tag weight ranged between 0.1 and 1.7% of fish body weigh; Table 1) was inserted through the incision and pushed forward into the body cavity. The incision was closed using two independent sutures. Mean surgery time was 2.8 min. In addition, scales for age determination and stable isotope analysis (see below) were taken during anaesthesia. All pike were released immediately after recovery from anaesthesia on the same location in each lake regardless of their capture location. Fish were captured and tagged between 5 and 10 May 2015.
Macrophyte sampling
To obtain an assessment of macrophyte assemblage and coverage, two SCUBA divers visually assessed macrophyte occurrence along 25 (HSC lake) and 26 (LSC lake) transects at the end of June and the beginning of September 2015. Sampling considered both aquatic plants as well as submerged dead terrestrial plants (only present in the LSC lake). Transects were situated from the shore to a depth of 12 m or deeper when macrophytes were present there. The coverage of each macrophyte species, the uncovered bottom area, the percentage composition of each species and maximum and minimum height of each macrophyte species were measured at 1 m depth intervals. The height of macrophytes was measured using a measuring tape. Dead flooded terrestrial plants were mostly European elder Sambucus nigra and thus categorized as a single group. Structural Complexity Index (hereafter referred to as SCI) in each lake was calculated to compare habitat complexity between lakes and its development during the study period. Calculation of SCI was based on information from the 25 and 26 (HSC and LSC lakes, respectively) scuba diver macrophyte assessment transects (see above). Species coverage and macrophyte height were both considered for calculation of the index. Species coverage and macrophyte height along each transect were georeferenced in a GIS-environment, and overlaid the bathymetric map for the lake. Height of macrophytes was discretized into 100 bins (from 0 to 200 cm by 2 cm) and coverage was discretized into 100 bins by 1%. The structural complexity index SCI was then defined as discretized height multiplied by discretized coverage. This gave a SCI range from 0 (no occurrence of macrophytes) to 10 000 (macrophyte coverage 100% and height ≥ 2 m).
Temperature and oxygen measurement
To obtain abiotic parameters which can drive pike spatial distribution36,37, we monitored water temperature and oxygen concentration in both lakes. Water temperature was monitored using 60 data loggers (Onset, USA, HOBO Pendant temp/light 64 K). Data loggers were placed at two sites in each lake in order to cover east/west (HSC lake) and south/north (LSC lake) gradients (Fig. 1). At each site, data loggers were attached to a rope in 1 m intervals spanning from the surface to 13 m (14 data loggers) with one extra data logger located at a depth of 20 m. The rope was tied to a floating buoy anchored at 22 m depth. This setup ensured both dense coverage in depths of rapid temperature change and, with a measurement interval of 5 min, high spatiotemporal resolution of the temperature profile. Oxygen concentration was measured in each lake (once a month in the HSC lake, and once during the observed period in the LSC lake) by calibrated YSI 556 MPS probe (YSI Incorporated, USA). Measurements were performed close to the western (HSC lake) or northern (LSC lake) data logger station (Fig. 1).
Fish community sampling
To obtain data of pike prey distribution, we performed spatially stratified fish community sampling by gillnet and hydroacoustic surveys. Gillnet surveys were conducted in September 2014 and 2015 at two localities in each lake in benthic habitats and one central locality in each lake for pelagic habitats, using 30 m long standard European multi-mesh gillnets38. At each locality in each lake, one series of three survey nets were set in the benthic and pelagic habitats at depths 0–3, 3–6, 6–9 and 9–18 m. In the deeper LSC lake, series of three survey nets were also set at depths 18–24, 24–30 and > 30 m. Benthic and pelagic gillnets were 1.5 m and 3 m high, respectively. Gillnets were set overnight, i.e. installed 2 h before sunset and lifted 2 h after sunrise39. Only catches of fish older than young-of-the-year were considered for this study. YOY fish were excluded as their density can be largely underestimated in gillnet surveys40. Moreover, Vejřík et al.35 found only prey fish larger than 10 cm in the stomachs of pike caught from both study lakes. Catches were expressed as catch per unit of effort measured as number of fish caught per 1000 m2 gillnet area per night (NPUE), and as kilogram fish 1000 m−2 night−1 (BPUE).
The acoustic surveys were performed both during the day and night, using a calibrated Simrad EK 60 echosounder operating at frequency of 120 kHz and following a pre-set zig-zag cruise track. The transducer was mounted 0.2 m below water surface, beaming vertically downwards. Recorded data were analysed using Sonar5-Pro software version 6.0.341, using a Sv-threshold of − 62 dB (thresholded at 40 logR) and a target strength (TS) threshold of − 56 dB (corresponding to fish with an approximate total length of 4 cm42). Shoals were detected manually while fish tracking was used for individual fish. The following settings were used in the Sonar5-Pro auto-tracking tool to select fish tracks: minimum track length (MTL), maximum ping gap (MPG) and vertical range gating. MTL was set to 3 echoes, MPG was set to 1 and vertical gating to 0.15 m for the whole water column. Only areas deeper than 5 m were included in further analysis. The relative fish density was calculated for shoals (shoals ha−1) and individual fish (ind. ha−1) using the number of shoals or fish related to the surveyed area (calculated from the wedge formed by the distance, the range and the opening angle). The analysed water column was divided into depth layers (1 m thick, from 5 m below the surface to the bottom).
Data processing
Positions of individual fish were calculated using a proprietary post-processing software UMAP v.1.4.3, based on multilateration of the time-difference-of-arrival (TDOA) of the acoustic signal received at different telemetry loggers (Lotek Wireless Inc., Newmarket, Ontario, Canada). Positions calculated using UMAP software can contain position duplicates and the use of TDOA for positioning implies large errors in a proportion of the position estimates43. Therefore, a position filter was applied in order to remove duplicate positions and positions with large errors. A detailed description of the filtering procedure is given in the SM1 (sec. Filtering of positions estimated by the U-MAP software). The position estimates (unfiltered and filtered) and depth profiles of each individual fish was visually inspected. If both horizontal and vertical locations became constant without latter movement, this individual was considered dead or expelling the tag (2 ind. in the LSC lake and 3 ind. in the HSC lake) and removed from further analyses.
The detection probability of transmitter detection and thus horizontal positioning changed over time and space due to underwater structures and thermal stratification. Numbers of position estimates for each individual ranged from 0 to 2884 per day (from 0 to 83% of possible number of position estimates). In order to reduce potential bias resulting from the varying positioning rate44, the individual data was regularized by calculating mean position for every 15-min interval (q-position hereafter). The q-position mean depth was calculated as the 15-min mean depth for each detected transmission with sensor depth reading. Since this did not require trilateration, the q-position time series could have a depth value without having a position estimate. This gave a maximum of 96 q-positions per day for each individual. If no position estimate was obtained within a 15-min interval, the q-position was interpolated between q-positions close in time if several conditions were fulfilled; (a) the time difference between the two q-positions used for interpolation (last before and first after the interval(s) to be interpolated) had to be shorter than 2 h, (b) the distance between these q-positions was less than 100 m, and (c) the depth difference between the q-positions was less than 2 m. Under such conditions, we assumed that the fish movement was limited between the q-positions used for interpolation (likely due to resting in an area with poor positioning coverage), and that linear interpolation in order to regularize the time series was therefore justified. For each position estimate, distance to bottom was calculated as the bottom depth at that position subtracted by the sensor depth. Distance to bottom for each q-position was calculated as the mean of each distance to bottom estimate within the 15-min interval. Individual and total yield of positions and calculated q-positions during the study period is given in Table 2.
Depth of fish was measured by an internal tag sensor and transmitted together with an tag identification. Therefore, to get depth of fish, only detection of a single receiver was required (not three as for location of a horizontal position). Depth location was, therefore, obtained in much higher detail than fish horizontal location (fewer gaps). To synchronize depth and position with horizontal position, we calculated mean depth for the same 15-min intervals as we used for calculation of mean horizontal position from receiver detections at that interval.
Extent of horizontal area use was calculated using a 95% kernel utilization distribution (hereafter referred to as dH-KUD). dH-KUD was calculated for each day and individual separately and only for days with more than 12 daytime and 12 night-time q-positions. Night-time and daytime were defined as one hour after/prior to civil twilight periods39. Exact time of sunset and sunrise in each day was calculated using R-package “maptools”45 and dH-KUD was calculated using R-package “adehabitatHR”46. Parameters required for calculation of dH-KUD were set as follows:a simplified lake shape polygon was used as a boundary, a raster of a lake with a cell dimension of 10 × 10 m was used as a grid and the smoothing parameter h was set to value of 50. The extent of vertical movement was evaluated separately. Vertical space use (hereafter referred to as dV-KS) was calculated as a one-dimensional kernel estimate (95%) of Gaussian density function with the smoothing bandwidth parameter set to 0.4 fitted on distribution of utilized depths.
Activity of fish was calculated as horizontal swimming speed (expressed in body lengths per second, BL * s−1) and vertical swimming speed (depth change per second, m * s−1) between two consecutive q-positions.
To test the importance of an open water habitat for pike in both lakes, the proportion of time spent in open water was calculated (TOW). Each q-position was assigned to be either in the benthic (distance < 5 m from the bottom) or in the open water habitat (≥ 5 m from the bottom).
Daily water temperature and day length were abiotic factors considered to potentially drive pike behaviour in lakes36,37,47. Daily water temperature was calculated as the mean temperature from measurement of all data loggers at depths 0–3 m for each date during the study, separately for each lake. This parameter reflects both rapid daily and gradual seasonal temperature changes. Day length was calculated as time between sunrise and sunset.
Stable isotopes and growth
Stable isotopes are widely used in studies of food-web structure and function, as well as individual specialization among consumer populations48,49. Here, we used stable carbon (δ13C) and nitrogen (δ15N) isotopes to estimate the relative reliance of individual pike on littoral carbon (food) resources (hereafter abbreviated as littoral reliance, LR). The LR estimates were calculated using the two-source isotopic mixing model described in Post50, where the δ13C value of a consumer (measured from the outermost annual ring of pike scales) is compared to those of littoral and pelagic isotopic end-members. We used δ13C values of littoral macrophytes and benthic algae, and pelagic particulate organic matter (POM) as the littoral and pelagic isotopic end-members, respectively, because some pike individuals and herbivorous prey fishes had considerably higher δ13C values than littoral benthic invertebrates. For more details of sample collection and preparation for stable isotope analyses, see the previous studies of predatory35 and generalist17 fishes in the two study lakes.
Age determination and growth calculation for each individual were conducted using scale reading. Three scales were read for each individual, results were then averaged and used for back-calculation of size-at-age of each individual using the Fraser–Lee Equation51. Only the body increment during the year prior to tagging was used to test for the correlation of individual growth with the rates of horizontal and vertical movement, since this was the growth increment most relevant to the activity and body size during the study year.
Statistical analyses
Behavioural traits and use of the pelagic zone
The response variables analyzed were horizontal area use (dH-KUD), vertical space use (dV-KS), daily mean depth (hypothesis i); horizontal activity, vertical activity (hypothesis ii); and time spent in open water (TOW, hypothesis iii); all of which could be associated with either a foraging or cover context. To test whether individuals from both lakes exhibited consistent behavioural differences across time we fitted a series of random-effects models and selected the best one based on Akaike Information Criterion (AIC, hypothesis iv). These models included a different set of covariates previously selected through Random Forest (hereafter RF) regression analysis52 using the R package randomForest 53. RF is suitable for the identification of relevant interactions between variables with low to negligible effects in isolation (further description in SM1, sec. Identifying interactions with Random Forest), and thus it would allow us to constraint both the number and potential interactions between them.
For the analysis of dH-KUD, dV-KS, mean depth horizontal activity and vertical activity we first fitted full random-intercepts and random-slopes mixed-effects models (LMMs)54,55 with Gaussian error distribution using the package nlme54. These models all include the same relevant predictors and interactions from the RF analysis, which yielded quite consistent results between the analyzed variables. Fish body length and daily water temperature were included, prior transformation by z-score standardization56, as continuous covariates in the model. The factor Lake was used as a proxy of structural habitat complexity (SHC). To account for the individual variation in repeated-measures (see SM1 for more details), we added the unique identifier (tag identification number) as a random intercept. We included the time × Lake interaction with time both as a random slope effect and a fixed effect to model mean between-lake temporal trends. Time as a fixed effect represents the overall effect of time among individuals in the two lakes while time as a random intercept measures the variance of the temporal effects in response across individuals (i.e., repeatability)57. This random-effects structure was also always selected in alternative analyses using log-likelihood ratio tests58,59.
TOW was analysed using zero–one beta inflated models60,61 fitted with the function gamlss() in R package gamlss62 by specifying a BEINF distribution. GAMLSS models are a particular type of GAM for Location, Scale and Shape which allow mixed distributions of continuous observations in the range 0–1 and discrete values at 0 and 1 often representing probabilities ruled by different processes63. Models also allowed evaluation of the probabilities at the extremes. We assumed that the probability that the time spent in open water is zero (p0) is associated with a behavior involving inexistent use of the pelagic habitat. On the other hand, the probability at one (p1) would imply that fish fully use the pelagic habitat. In the range between 0 and 1, TOW is modelled as a i.i.d. Gaussian variable indicating variable levels of pelagic habitat use. These probabilities are therefore governed by different processes reflecting variations in behavior and are integrated in a single model through four distribution parameters (μ, σ, ν, τ) (see SM1 for more details on GAMLSS model parametrization; sec. GAMLSS model of pelagic habitat use). Each component was modelled as a function of the between-lake temporal trends (time × Lake interaction) controlling for the effect of dH-KUD and dV-KS (as well as body length) on TOW to account for the fact that large differences in vertical and horizontal range influence pelagic habitat use. Each continuous covariate was initially fitted using a penalized P-spline smoothing function, denoted by pb() in R syntax, which automatically selects the smoothing parameter (indicative of the degree of smoothness or complexity of the fitted curve). In the final models not all distribution parameters were necessarily modelled using an additive term as it depended on decisions made during the model selection procedure (see SM1 for further details on GAMLSS model parametrization, sec. GAMLSS model of pelagic habitat use).
To account for past error residual correlations and prevent pseudo-replication55, all models were fitted using an autocorrelation structure. This allowed us to further assess the differences in temporal trends by correctly partitioning inter- and intra-individual sources of variance54,55 whilst accounting for slowly fluctuating trait values or behavioural lags about those timelines. We selected the best autocorrelation structure by comparing models according to AIC64,65 and by previously visualizing ACF and PACF residuals to determine the starting lag value (rho)66. Along with the individual identity, time was included as a continuous variable to control for the correlation of residuals between any given days.
In a second step, once a full model was fitted, we created sets of candidate models with a different number of predictors and ranked them according to their AICc values (i.e., AIC with finite-size correction67, lower is better) and weights (higher is better), and finally compared. Here, the importance of the interaction terms detected in the RF analysis was further tested using likelihood ratio tests (for more details see SM1). The ranking and selection of LMMs was conducted using the R package AICcmodavg68. For GAMLSS models we used the different functions in the gamlss package to conduct a stepwise selection of appropriate terms for each distribution parameter (μ, σ, ν, τ) (see SM1 for a detailed description of the GAMLSS selection procedure utilized, sec. GAMLSS model of pelagic habitat use).
As a last step, models were re-fitted using restricted maximum likelihood (REML) and variance components assessed with the package MuMIn69. We calculated between- and within-individual variance components and estimated the repeatability (R) index, both on data from the two lakes and on subsets of data for each lake, to evaluate qualitative differences of repeatable patterns in temporal trend lines. Repeatability was calculated as:
$$frac{{{text{V}}_{{{text{tag}}_{text{ID}}}} }}{{left( {{text{V}}_{{{text{tag}}_{text{ID}}}} + {text{V}}_{{text{e}}} } right)}},$$
where (Vtag_ID + Ve) is the total phenotypic variance (Vp), resulting from adding the between-individuals variance (Vtag_ID) to the within-individual (residual) variance (Ve). A high R may indicate high between- and/or low within- individual variability in the outcome values (see SM1 for further details on repeatability index). As a separate measure of repeatability, we compared the individual ranks in horizontal space use between June-July, June–August, and June–September using the Spearman rank correlation test70. The same procedure was done for swimming activity.
Stable isotopes
We ran linear models to test if the trophic niche of pike was related to individuals’ behaviour (open-water use as a proxy), structural habitat complexity (SHC; lake factor as a proxy) or body length (i.e., total length). Prior to modelling, open-water use was logit-transformed and continuous explanatory variables were scaled (see above).
Growth
We used linear regression to find differences in individual growth (hypothesis v) in the previous season (log-transform of body increment in year 2014; dependent variable) between lakes and whether they were determined by the space use (mean daily dH-KUD and dV-KS), body length and age of fish (explanatory variables).
Ethics approval and consent to participate
This study complied with and was approved by the Animal Welfare Committee of the Biology Centre CAS (45/2014) according to § 16a of the Act No. 246/1992 Coll., on the protection of animals against cruelty, as amended. The study was carried out in compliance with the ARRIVE guidelines.
Consent for publication
This manuscript presents work that has not been published and is not under consideration for publication elsewhere. All authors involved in the manuscript have agreed to be listed and contributed to the research reported.
Source: Ecology - nature.com
