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 Czech Republic (Fig. 1). Aquatic restoration lasted 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 4–9 m. 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 summer 32.
Macrophyte sampling was carried out prior to the study in September 2014 and May 2015 33. Macrophytes were dense only in the HSC lake where they covered 60-91% of 0-12 m deep inshore areas. In the LCS 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 abundant 34. 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 details are in the Supplementary material, herein referred as SM, sec. Stocking of pike). Stocked pike came from various breeding ponds in Central Bohemia region.
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, 47 receivers in LSC lake) deployed in arrays with distances between 3 nearest receivers ranging from 203 to 288 m (mean 251± 18 m) in 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-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 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 1, 5 and 13 m, and 2 nearshore locations at depths 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. The systems were installed from April 2015 to March 2016 and data were 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 highest fish activity and development of macrophytes.
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 HCS and 86 cm/4.15 kg for 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, 65x28 mm, mass in air of 13 g, including pressure and motion sensors, burst rate 25 s) 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. Further, 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 to 10 May 2015.
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 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 (herein referred as SCI) in each lake was calculated to compare habitat complexity between lakes and its development during study period. Calculation of SCI was based on information from the 25 and 26 (HSC and LSC lake, respectively) scuba diver macrophyte assessment transects (see above). Species coverage and macrophyte height were both considered for calculation of the index. A detailed calculation procedure can be found in the SM (sec. Calculation of Structural Complexity Index).
Temperature and oxygen measurement
To obtain abiotic parameters which can drive pike spatial distribution 35,36, we monitored water temperature and oxygen concentration in both lakes. Water temperature was monitored using 60 data loggers (Onset, USA, HOBO Pendant temp/light 64K). Data loggers were placed at two sites in each lake in order to cover east/west (HSC) and south/north (LSC) gradients (Fig.1). At each site, data loggers were attached to a rope in 1 m intervals spanning from the surface to 13 meters (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 minutes, high spatiotemporal resolution of the temperature profile. Oxygen concentration was measured in each lake (once a month in HSC, and once during observed period in LSC) 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 gillnets 37. 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 sunrise 38. Only catches of fish older than young-of-the-year were considered for this study. 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 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.3 39, using a Sv-threshold of -62 dB (thresholded at 40 logR) and a target strength (TS) threshold of -56 dB (corresponding to fish of an approximate total length of 4 cm, 40). 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 whole water column. Only areas deeper than 10 m were included in further analysis. The relative fish density was calculated for shoals (shoals km-1 cruise track) and individual fish (ind. km-1 cruise track).
Positions of individual fish were calculated using a proprietary post-processing software UMAP v.1.4.3, based on multilateration of 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 estimates 41. Therefore, a position filter was applied in order to remove duplicate positions and positions with large errors. A detailed description of filtering procedure is given in the SM (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 LSC lake and 3 ind. in HSC lake) and removed from further analyses.
To obtain representative and balanced daily number of horizontal positions and depth locations for each fish, mean 15-minutes position (herein referred as q-position) and depth were calculated. Detailed description of the calculation is given in the SM (sec. Calculation of q-position and Calculation of depth location). Extent of horizontal area use was calculated using a 95 % kernel utilization distribution (herein referred 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 periods 38. Exact time of sunset and sunrise in each day was calculated using R-package “maptools” 42 and dH-KUD was calculated using R-package “adehabitatHR” 43. Parameters required for calculation of dH-KUD were set as follows: simplified lake shape polygon was used as a boundary, raster of a lake with dimension of cells 10x10 m was used as a grid and smoothing parameter h was set to value of 50. Extent of vertical movement was evaluated separately. Vertical space use (herein referred as dV-KS) was calculated as one-dimensional kernel estimate (95%) of Gaussian density function with 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 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 < 5m 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 lakes 35,36,44. Daily water temperature was calculated as mean temperature from measurement of all data loggers at depths 0-3 meters 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 populations 45,46. 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 Post 47, where the δ13C value of a consumer (here 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 predatory 48 and generalist 32 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 equation 49. Only the body increment during the year prior to tagging was used to test for 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.
Behavioural traits and use of pelagic zone
We tested whether individuals from both lakes exhibited consistent behavioural differences across time regarding horizontal area use (dH-KUD), vertical space use (dV-KS), daily mean depth, horizontal activity, vertical activity and time spent in open water (TOW), which could all be associated with either foraging or cover contexts. Before model building, we selected potentially important predictors influencing variables being modelled. Except for the analysis of TOW, which included dH-KUD, dV-KS and body length as predictors, the other variables were analysed using fish body length (i.e., total length) and environmental factors such as daily water temperature included as continuous covariates. Lake was considered used as a proxy of structural habitat complexity (SHC) and included as a categorical factor. Continuous predictors were normalized using z-score standardization (set to a mean of 0 and standard deviation of 2, 50) to ease comparison between continuous predictors and with the untransformed factor Lake.
First, we investigated interactions between potential predictors by running a Random Forest (hereafter RF) regression analysis 51, which is suitable for the identification of relevant interactions between variables with low to negligible effects separately (further description in SM, sec. Identifying interactions with Random Forest).The RF analysis was conducted using the R package randomForest 52. After selection of important interaction terms, we fitted a generalized additive mixed model (GAMM) in order to test whether predictors had a linear or non-linear relationship with response and whether smooths terms (e.g., cubic splines) were necessary to be implemented for some covariates 53. This model included relevant interaction terms from previous analysis as fixed factors and a simple random-intercept at the individual level as unique identifier (tag_ID) to account for individual repeated-measures. An analysis of variance (ANOVA) on the GAMM output allowed discerning of whether the relationship between response and predictors was linear 54. To fit GAMM models, we used the bam function (R package mgcv, 55), which is specifically designed for large datasets 56.
If linearity was assumed, we ran separate linear mixed models (LMM) 57,58 including the time × Lake Interaction, with time (as a daily continuous variable) modelled both as a random slope effect and as a fixed factor along with a random intercept of tag ID (see SM for more details). This random effects structure was always selected in preliminary analyses using likelihood ratio tests 59,60 and models mean between-lake temporal trends rather than individual temporal trend lines which are more related to an individual by time interaction (i.e., reaction norms61). Time as a fixed effect represents the overall effect of time among individuals in the two lakes, which is a function of the grouping structure (i.e., the slope may vary over time but not necessarily between individuals), while time as a random intercept measures the variance of the temporal effects in response across individuals (i.e., repeatability) 62. For the analysis of continuous variables, we fitted linear random-intercept and random-slope mixed-effects models with Gaussian error distribution using the R package nlme 57. TOW was analysed using zero-one beta inflated models 63,64 fitted with the function gamlss in R package gamlss 65 specifying distribution BEINF. 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 zero and one often representing probabilities ruled by different processes 66. We assumed that the inexistent use of the pelagic habitat (i.e., stationarity) was associated with the probability at zero (p0) while the complete use of it with the probability at one (p1) (Equations 3, 4 in SM, sec. GAMLSS model of pelagic habitat use). On this basis, we modelled the four distribution parameters (µ, σ, ν, τ) as a function of the between-lake temporal trends (time × Lake interaction), behavioral and phenotypic predictors (dH-KUD, dV-KS and body length) and, where appropriate, smoothing functions. However, it should be noted that all distribution parameters were not necessarily modelled using additive terms and it depended on decisions made during the model selection procedure (see SM for further details on GAMLSS model parametrization, sec. GAMLSS model of pelagic habitat use). Since GAMLSS models already test the assumption of linearity between dependent and explanatory variables, we skipped the preliminary fitting and analysis of a GAMM model.
In addition to the random effects, models were fitted using an autocorrelation structure to account for past error residual correlations and prevent pseudoreplication 58. This would allow further unbiased evaluation of the differences in temporal trends by correctly partitioning inter- and within-individual sources of variance 58,67 while accounting for slowly fluctuating trait values (i.e., behavioral lags) about those timelines. We compared models with different autocorrelation structures using Information theory 54,68 and through preliminary inspection of their ACF and PACF residuals to determine the starting values (rho parameter) 69. Along with the individual identity (tag_id), time was included as a continuous variable to control for the correlation of residuals between any given days.
The base model was further used to determine the overall support to our hypotheses with respect to variation explained by the fixed effects. We added predictors one by one, thus creating sets of candidate models nested within the previously fitted model. We used maximum likelihood to obtain unbiased estimates of significance of the fixed effects. Next, we used information theory to rank models from best to worst average representation according to their AICc values (i.e., AIC with finite-size correction, 70) (lower is better) and AICc weights (higher is better). Lastly, we compared those models and weighed the importance of the interaction terms detected in the RF analysis using likelihood ratio tests (for more details see SM1). The selection and classification of models was conducted using the R package AICcmodavg 71. For TOW, we used the different functions in the gamlss package to perform stepwise selection of appropriate terms for each distribution parameter (μ, σ, ν, τ) (see SM for a detailed description of the GAMLSS selection procedure utilized, sec. GAMLSS model of pelagic habitat use).
Models were re-fitted with restricted maximum likelihood for the random-effects and autoregressive terms. We used the R package MuMIn 72 for assessing variance components and to estimate marginal and conditional R2 values. From the fitted models, we estimated inter- and within-individual variance components and calculated individual repeatability (R) both on data from both lakes and on subset data for each lake to evaluate qualitative differences of repeatability in temporal trend lines across lakes. Repeatability was calculated as the proportion of variance explained at the individual level relative to the total variance following the expression A high R value either indicates high inter-individual variation or low within-individual variability in the values of response, i.e., consistency in individual repeated-measures across time (see SM for further details on repeatability estimates). As a separate measure of repeatability, we compared the individual ranks in horizontal space use between June-July, June-August, and June-September using Spearman rank correlation test 73. The same procedure was done for swimming activity.
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).
We used linear regression to find differences in individual growth 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).