Effect of temperature on zooplankton vertical migration velocity

Zooplankton diel vertical migration (DVM) is an ecologically important process, affecting nutrient transport and trophic interactions. Available measurements of zooplankton displacement velocity during the DVM in the field are rare; therefore, it is not known which factors are key in driving this velocity. We measured the velocity of the migrating layer at sunset (upward bulk velocity) and sunrise (downwards velocity) in summer 2015 and 2016 in a lake using the backscatter strength (VBS) from an acoustic Doppler current profiler. We collected time series of temperature, relative change in light intensity chlorophyll-a concentration and zooplankton concentration. Our data show that upward velocities increased during the summer and were not enhanced by food, light intensity or by VBS, which is a proxy for zooplankton concentration and size. Upward velocities were strongly correlated with the water temperature in the migrating layer, suggesting that temperature could be a key factor controlling swimming activity. Downward velocities were constant, likely because Daphnia passively sink at sunrise, as suggested by our model of Daphnia sinking rate. Zooplankton migrations mediate trophic interactions and web food structure in pelagic ecosystems. An understanding of the potential environmental determinants of this behaviour is therefore essential to our knowledge of ecosystem functioning.


Introduction
Zooplankton diel vertical migration (DVM) is an ecologically important process, driven by internal and external drivers, and affecting nutrient transport and trophic interactions in lakes and oceans (Hays, 2003;Ringelberg, 2010). In a typical DVM, organisms migrate from deep and cold waters towards the warmer surface at sunset, and they descend at sunrise (Ringelberg, 2010;Williamson et al., 2011). Predation, light, food availability, and temperature are recognised as the main DVM drivers. However, avoidance of visually hunting juvenile fish and food requirements are generally considered the most significant factors (Ringelberg, 1999;Hays, 2003; Van Gool & Ringelberg, 2003;Williamson et al., 2011). According to this reasoning, zooplankton hide during the daytime in the deep and dark layers of the lake in response to chemical cues released by fish (kairomones) (Dodson, 1988;Neill, 1990;Loose & Dawidowicz, 1994;Lass & Spaak, 2003;Boeing et al., 2003;Beklioglu et al., 2008;Cohen & Forward, 2009). At night, when mortality risk from visual predators is low, zooplankton migrate to surface waters to graze on abundant phytoplankton. Food availability can modify migration behaviour as well ( Van Gool & Ringelberg, 2003); if a deep chlorophyll maximum is present, zooplankton may not migrate (Rinke & Petzoldt, 2008). Finally, DVM can still take place in very transparent and fish-less lakes where zooplankton hide to prevent damage from surface UV radiation (Williamson et al., 2011;Tiberti & Iacobuzio, 2013;Leach et al., 2014).
Several laboratory (Daan & Ringelberg, 1969; Van Gool & Ringelberg, 1998a, b;Ringelberg, 1999; Van Gool & Ringelberg, 2003) and field observations (Ringelberg & Flik, 1994) suggest that the speed at which zooplankton migrate is influenced by the drivers affecting the DVM. Variations in this velocity have been used as a proxy to infer possible reactions and behavioural changes in zooplankton populations. This speed has two components: (1) the swimming velocity (SV), which is the organism's instantaneous velocity during a reactive phase (Ringelberg, 2010); and (2) the displacement velocity (DV), which is the vertical displacement of the organism divided by the time taken to perform the movement (Gool & Ringelberg, 1995). The DV is smaller than SV because it combines various animal reactions, including latent periods.
The dependence of the zooplankton SV on the presence of predators in the water column is species specific. Beck & Turingan (2007) showed that the swimming speed of brine shrimp (Artemia franciscana Kellog, 1906) and two copepod species (Nitokra lacustris Schmankevitsch, 1875 and Acartia tonsa Dana, 1849) increased when zooplankton are exposed to water with larval fish. However, this was not the case for the rotifer Brachionus rotundiformis Tschugunoff, 1921. Dodson et al. (1995 hypothesised instead that Daphnia may adopt a conservativeswimming behaviour, because faster swimming velocity would increase predator-prey encounter rates and therefore predation risk. In their experiment, they did not observe increases in SV of Daphnia pulex Leydig, 1860 when exposed to Chaoborus-enriched water. However, Weber & Van Noordwijk (2002) later showed that the swimming response of Daphnia galeata Sars, 1864 to info-chemicals can be clonespecific and that Perca kairomone can positively affect their speed. Finally, Van Gool & Ringelberg (1998a, b, 2003 observed that the Daphnia swimming response to changes in light intensity (S) can be enhanced by chemical signals from Perca fluviatilis Linnaeus, 1758 and by food concentration as well. Field observations, by Ringelberg et al. (1991) and Ringelberg and Flik (1994), showed instead that the DV at dusk of Daphnia galeata x hyalina was well correlated with S measured in the water column at the time of the migration. The response of the animals was strongly light-driven but not influenced by food concentration or water transparency. Moreover, acceleration and deceleration in S can also play a role in determining the Daphnia swimming response (Van Gool, 1997).
Laboratory observations have shown that organism swimming speed can also be greatly enhanced upon exposure to higher water temperatures. Temperature has a twofold effect. Warmer water enhances the biological and metabolic activity of organisms, increasing energy available for locomotion and therefore the power available for thrust generation (Beveridge et al., 2010;Moison et al., 2012;Humphries, 2013;Jung et al., 2014). At the same time, a higher water temperature reduces fluid kinematic viscosity (m) and drag on zooplankton beating appendages, so that they can swim faster (Machemer, 1972;Larsen et al., 2008;Larsen & Riisgård, 2009;Moison et al., 2012). Since this temperature-dependant behaviour has never been demonstrated in the field and during the DVM, it is not known whether this may be relevant during the vertical ascent.
With regard to the downwards DVM at sunrise, the zooplankton velocity may vary greatly depending on the swimming behaviour adopted by the organism. Descent can occur by active swimming, when organisms orient downwards, or by passive sinking (Dodson et al., 1997b;Ringelberg, 2010) with lower velocities controlled by organism's buoyancy and gravity. To date, little is known about which parameters really affect the choice of a swimming behaviour and velocity. The only available study is by Ringelberg & Flik (1994) who reported that Daphnia swim downwards only when the light stimulus S is high but this threshold is currently not known. However, no evidence was provided that the organisms really sank.
Available measurements of zooplankton swimming velocity during the DVM in the field are rare. The first estimates of zooplankton DV were made by Ringelberg & Flik (1994) using successive net hauls to measure the zooplankton position in the water column. This method only offers a coarse resolution of the organisms position and is time-consuming. More recent studies by Lorke et al. (2004) and Huber et al. (2011) instead measured the DV using velocity data and variations in the backscatter strength signal from an Acoustic Current Doppler Profiler (ADCP). However, their objective was not to analyse variations in time series of the zooplankton DV.
The objectives of this study were (1) to explore and explain for the first time seasonal variability in the DV of diel migration in the field; and (2) to understand which environmental parameters really drive the rate of zooplankton migration under real field conditions, when combined DVM drivers act at the same time.
Existing studies compared only the effect of one parameter at the time on the DV and assessed the zooplankton behaviour only with light-induced swimming responses. In this study, we continuously measured the velocity of the migrating layer (bulk velocity or mean DV) at sunset and sunrise along with the following DVM drivers measured in the field: water temperature, turbulence, chlorophyll-a concentration, light conditions, and zooplankton concentration and size during the DVM. We quantified v up as the bulk velocity at dusk, when zooplankton actively swim to reach the surface, and employed a correlative approach to infer the likely dependence of v up upon potential DVM drivers. At sunrise, when organisms usually sink towards the aphotic lake layer, the mean DV is referred as v down . This velocity was modelled and correlated with zooplankton density and size measured in the laboratory to verify whether organisms sank.

Study site
Measurements were taken in 2015 and 2016 in Vobster Quay, a shallow man-made lake located in Radstock, UK. The lake has an average depth of 15 m and maximum depth of 40 m (see Fig. 1). It has a simple bathymetry, with very steep shores and a flat bottom. The lake is oligotrophic, with an average chlorophyll-a concentration of about 1 lg l -1 . During summer stratification, the surface temperature reaches a maximum of 21°C, and the bottom temperature is approximately 9°C. The metalimnion usually extends from 5 to 17 m. The water transparency is very high, with an average Secchi depth of 10 m over the summer. The lake was stocked in August 2004 with a population of P. fluviatilis and Rutilus rutilus Linnaeus, 1758.

Acoustic measurements
Acoustic measurements were employed to track the position of the zooplankton layer during the DVM (e.g. Lorke et al., 2004) and to estimate their DV at sunset and sunrise using the acoustic backscatter strength. An acoustic Doppler current profiler (ADCP, Signature 500 kHz by Nortek) was bottom-deployed at location ''A'' in Fig. 1. The device was set up to record the acoustic backscatter strength (BS) from the vertical beam every 180 s for 90 s with a frequency of 0.5 Hz. The BS 1 m below the surface and above the bottom was removed due to surface reflection and the ADCP blanking distance. Data were available from 7 July to 27 July 2015, 24 June to 7 July 2016, and 21 July to 19 August 2016.
The BS was then converted to the relative volume backscatter strength (VBS) to account for any transmission loss of the intensity signal, using the sonar equation (Deines, 1999): where P dbw is the transmitted power sent in the water, L dbm the log10 of the transmit pulse length P ¼ 0:5 m (Nortek, pers. comm.) and R the slant acoustic range. a is the acoustic absorption coefficient estimated following Francois (1982) and using the temperature profiles from the thermistors chain (''T'' in Fig. 1).

Bulk velocity estimation
The bulk velocity of the migrating layer is defined as the slope of the zooplankton layer during the DVM. Fig. 2d shows an example of the migration on 2 July 2016, where the VBS in black depicts the zooplankton in the water column. When the zooplankton start swimming upwards after sunset (21:45 local time), a line can be fit to the migrating layer, whose slope is constant throughout the ascending phase (red line in Fig. 2d). This slope is the upward bulk velocity (v up ) and provides the mean DV of the organisms. The bulk velocity during the reverse DVM is referred as v down .
To objectively and automatically estimate v up and v down from the VBS data, an image-detection algorithm was developed in MATLAB to detect the acoustic shape of the zooplankton layer and to estimate its slope. Two input parameters need to be provided by the user: (1) the box (target area) containing the layer and delimited by Z WINDOW and T WINDOW limits (see Fig. 2a); when the DVM begins and ends and its limits can be immediately identified from the acoustic image; (2) the VBS range V range ¼ ½V min ; V max of the zooplankton in the migrating layer. This range is characteristic of the zooplankton for that day and is affected by the organism abundance and size. Because the VBS always reaches the maximum in the migrating layer, the algorithm can be controlled just by adjusting the lower limit V min and setting V max to the observed maximum of 85 dB.
Once these two groups of parameters are set, to detect the layer, the algorithm picks the points in the target area where V min VBS 85 dB (see red dots in Fig. 2b). The acoustic shape then can be identified by selecting the outer points of the layer (Fig. 2c). A line can be finally fitted to the red points to estimate v up . Fig. 2d shows v up;81 ¼ 2.43 mm s -1 for V min ¼ 81 dB.
Some of the VBS data, representing isolated patches of zooplankton or environmental noise within the target area, had to be manually removed. An  Fig. 2b, where the three blue patches show dense aggregations of zooplankton with high VBS not belonging to the migrating layer. These points would be selected by the algorithm and affect the velocity computation in the layer.
Because the choice of V min can be arbitrary, the algorithm was run by changing V min from 75 dB, which was the lowest observed value, to V min ¼ V max ¼ 85 dB, generating 11 different layer detection and fitting results. Two examples are shown in Fig. 3 for V min ¼ 75 dB and V min ¼ 85 dB. Each image can be then inspected manually to remove spurious results and bad fits, when at least one of the following conditions are met: (1) the algorithm selects points within the target area but outside the layer; this occurs when V min is too low, such as the case in Fig. 3a; (2) the fitted line cuts across the migrating layer rather than following its centroid; this condition is met when V min is too low or high. When V min ¼ 85 dB (Fig. 3b), too few points, which lack alignment with the centre of the target layer, were picked to correctly estimate the slope.
After excluding the fits that meet the rejection criteria previously defined, the velocity of the accepted fits can be bootstrapped to estimate the mean bulk velocity v up and its 95% confidence interval. For example, by choosing V min from 78 to 83 dB for the 2 July migration, v up ¼ 2:57 AE 0:1 mm s -1 (Fig. 3c). This procedure was applied for both 2015 and 2016, producing a time series of sunset migration velocities v up and sunrise velocities v down (Fig. 4). The ''Electronic Supplementary Material'' document reports in detail the results of the layer detection algorithm for each analysed day (Online Resources 1-7).

Temperature measurements
The vertical temperature stratification was assessed with a thermistor chain deployed in location ''T'' (Fig. 1). The chain consisted of 10 RBR thermistors (Solo T, RBR Ltd.) with a temperature accuracy of AE 0:002°C and 5 Hobo loggers (TidbiT v2, Onset Computer) with an accuracy of AE 0:2°C. Data were sampled every 5 min. Data were available from June to October 2015 and from May until mid-August 2016 (Fig. 5). The mean temperature that the zooplankton experienced during the ascent phase of the DVM can be obtained by overlapping the thermistor chain data with the position and timing of the zooplankton acoustic shape for each migration date (Fig. 6).

Turbulence estimation
To assess the effect of turbulence on zooplankton swimming velocity, profiles of temperature fluctuations were acquired at location ''S'' ( Fig. 1) with a Self-Contained Autonomous MicroProfiler (SCAMP), a temperature microstructure profiler, following the sampling procedure in Simoncelli et al. (2018). After splitting each profile into 256-point (25 cm) segments, turbulence was assessed in terms of turbulent kinetic energy (TKE) dissipation rates e. Dissipation rates were estimated for each segment using the theoretical spectrum proposed by Batchelor (1959): where q ¼ 3:4 is a universal constant (Ruddick et al., 2000), . By fitting the observed spectrum S obs to S B ðkÞ, it was possible to estimate k B and then e. Bad fits, that provided invalid e, were rejected and removed using the statistical criteria proposed by Ruddick et al. (2000).
On each date, an average of six turbulence profiles were acquired to ensure a statistically reliable estimation of e. Profiles were then time-averaged to obtain The daily-averaged profiles of e are reported as ''Online Resources 8'' in the ''Electronic Supplementary Material''. Profiles were then depth-averaged throughout the whole water column to obtain the time series of the dissipation rates of TKE (Fig. 7). Data from the profiles were not averaged based on the zooplankton position because turbulence data were not always available on the same days as the acoustic measurements (see black dots vs. green areas in Fig. 7).

Chlorophyll-a concentration
Chlorophyll-a concentration profiles were acquired at sampling station ''S'' (Fig. 1). Data were collected from approximately 0.5 m from the lake surface with a depth resolution of 2 mm between 5 h and 5 min before dusk. Chlorophyll-a was sampled five times in July-August 2015 and ten times between May and August 2016 (see black triangles in Fig. 8). At least five profiles were sampled on each date.
Profiles were acquired using a fluorometer mounted on the SCAMP and manufactured by PME. The sensor excites the water with a blue L.E.D. emitter at 455 nm in a glass tube to avoid interference with ambient lighting. It then measures the voltage from the light emitted from the chlorophyll with a photo-diode. On each date, voltage profiles were binned every 30 cm and time-averaged to obtain one profile. Voltage data were finally converted to chlorophyll-a concentration using a linear calibration curve provided in Online Appendix A.
Potential food availability during the DVM was estimated by vertically averaging the chlorophylla concentration from the available chlorophyll-a profiles (Fig. 9). Since chlorophyll-a concentration was not available for all the dates with the acoustic measurements (green areas in Fig. 9), a linear interpolation was used to estimate missing chlorophyll concentration data between the available observations (red line in Fig. 9). This is the final time series used in subsequent correlations.

Underwater light conditions
To characterise the underwater light conditions, surface solar radiation I 0 was measured with a frequency of 5 minutes by a meteorological station located in Kilmersdon (Radstock, UK), 2 km from the site. Data were available for the studied period with the exception of 4 July 2016. Photosynthetically active radiation (PAR) profiles I PAR were also collected at the same time as chlorophyll-a concentration profiles. The PAR was measured with a Li-Cor LT-192SA underwater radiation sensor mounted facing upwards on the SCAMP. On each sampling date, I PAR profiles were filtered with a moving average filter to remove noise and then binned and time-averaged to obtain a mean PAR profile I PAR (Fig. 10). Profiles collected on 16 July 2016 were discarded because no valid data were recorded. A 20-cm Secchi disc was employed to measure the Secchi depth z SD in June-July 2015 and April-August 2016 (Fig. 11). The light extinction coefficient K was then estimated from PAR (K PAR ) and Secchi disc data (K SD ). The I PAR profiles were fitted the Lambert-Beer equation to estimate K PAR (Armengol et al., 2003): where I PAR;0 is the surface PAR (Fig. 12, black circles). The ambient attenuation coefficient K SD was determined from the Secchi disc depth z SD (Fig. 12, empty triangles) using the following empirical relationship (Armengol et al., 2003): Finally, the underwater light conditions at the depth where the zooplankton started the DVM can be described in terms of relative change in light intensity S (Ringelberg & Flik, 1994). This quantity can be estimated by using the definition provided by Ringelberg (1999) and assuming that the light intensity I in the water column decayed exponentially with the depth z: where t is the time, I 0 the surface solar radiation and K the light extinction coefficient. If K is time-independent, S is depth-independent and reduced to: By employing a linear correlation to the time series of the natural log of the solar radiation (ln I 0 ðtÞ) during the DVM, a slope can be estimated that provides the mean value of S (Fig. 13).

Zooplankton concentration
Zooplankton samples were also collected to characterise the zooplankton population density and vertical length approximately 50 cm) was used for this purpose. The water column was sampled at location ''S'' (Fig. 1) before the DVM, by towing the net through consecutive layers of thickness h ¼ 3 m. Three samples (n ¼ 3) were collected and pooled for each stratum. Samples were then stored in a 70% ethanol solution to fix and preserve the organisms and prevent decomposition (Black & Dodson, 2003). Four tows were acquired in June-July 2015 and 10 tows from   Four zooplankton groups were distinguished: Daphnia spp., adult and copepodite-stage copepods, small Cladocera and copepod nauplii. The concentration of each group was then estimated by dividing the abundance by the filtered water volume V ¼ p=4 Á d 2 Á h=n (Fig 14). The complete profiles of zooplankton concentration are reported in the ''Electronic Supplementary Material'' as ''Online Resources 9 and 10''.
Net tows do not, however, allow estimation of zooplankton concentration within the migrating layer and during the DVM with adequate time and vertical resolution. The VBS can be used to provide suitably resolved data on zooplankton concentration and size (Sutor et al., 2005;Record & de Young, 2006;Hembre & Megard, 2003;Huber et al., 2011;Barth et al., 2014;Inoue et al., 2016) and was available during the DVM. For each migration dates, the VBS was depth-averaged within the zooplankton layer to provide the mean zooplankton concentration during the vertical ascent (Fig. 15).

Zooplankton position
The VBS data can be further used to extract more information about the zooplankton behaviour in relation to the DVM. The acoustic shape (Fig. 2c) from the image-processing algorithm also provides the zooplankton position in the water column prior the DVM and when the migration ends. By taking the average of the upper and lower limits of the shape when the DVM begins (Z s, lower and Z s, upper ), it is possible to infer the average depth range within which the organisms were found before moving (Fig. 16, black dots). By repeating the same operation when the DVM stops, the two depths Z e, lower and Z e, upper once averaged provide the mean location where the organisms completed the migration and started spreading out near the surface (Fig. 16, empty dots).

Modelling of v down
The zooplankton sinking velocity U sink was modelled with Stokes-like equations and by using organism size and density measured from laboratory experiments. The estimated value of U sink can then be compared to the downwards velocity v down to independently verify if zooplankton sank during the descending phase of the DVM. We limited our attention to Daphnia because its hydrodynamics can be easily described by mathematical models, in contrast to other migrating genera. No data exist on the rate of sinking of live copepods. Moreover, they sink frequently with the body orientated in a horizontal or oblique position (Paffenhöfer & Mazzocchi, 2002) and that angle is not known. Daphnia sink instead mostly vertically (Dodson et al., 1995).

Model
If zooplankton passively descend during the DVM, U sink depends only on gravity, fluid characteristics and individual size and shape (Ringelberg & Flik, 1994;Kirillin et al., 2012). Moreover, assuming that zooplankton individuals do not interfere with each other while sinking, U sink can be approximated by using the average sinking rate of a single organism. The sinking rate U sink;s of a spherical Daphnia, as a function of the water temperature T, can be found by balancing the Stokes drag force exerted by the fluid (F ¼ 6plRU s ) with the body buoyancy (B ¼ Vg q d À q w ð Þ =q w ): where V is the animal volume, R the body radius, g the gravitational acceleration, q d the Daphnia density, and q w ðTÞ and lðTÞ the water density and dynamic viscosity, respectively. However, Daphnia are not spherical and Eq. 7 was thus improved by assuming that Daphnia are an ellipsoid with width w ¼ 2a e and length l ¼ 2b e (Fig. 17b). The new body drag is: where b e ¼ b e =a e . The new mean velocity is: Accounting for the additional drag from the antenna, the two appendages can be modelled as two cylindrical needles with dimensions a c and b c (Fig. 17b). The new equation is: where and b c ¼ b c =a c .

Daphnia size
In order to parametrize the above model, morphometric data on twenty specimens were measured under a Leica M205C dissecting microscope. The measured  (Fig. 17c) are: b e is half the body length l and a e half the body width. The arm length is b c and its width a c .

Daphnia density
The density of Daphnia carcasses q d was assessed using the density gradient method proposed in Kirillin et al. (2012). Nine different fluid densities were created in different cylinders by mixing freshwater with different volumes of sodium chloride solute (1 g ml -1 ). The cylinder densities were verified with a portable densimeter (DMA 35 by Anton Paar) after the water reached the equilibrium with room temperature.
Twenty specimens were then individually tested. Each organism was gently injected and moved from the cylinder with the lowest density (q w ¼ 1000 kg m -3 ) to the one with the highest density (q w ¼ 1050 kg m -3 ) until the cylinder where the animal reached the neutral buoyancy was found. When the organism stops sinking, it can be assumed that the carcass has a density between the fluid density of the cylinder where it floats and the previous cylinder where it sank. The results are reported in Fig. 18.

Bulk velocity
In 2015 the velocity at sunset ranged from a minimum of 2.23 mm s -1 on 12 July to a maximum of 3.55 mm s -1 on 22 July (Fig. 4, upper panel). In 2016 the time series of v up showed similar values with a lowest of 2 mm s -1 at the beginning of June and a maximum of 4.5 mm s -1 in mid-August (Fig. 4, lower panel). Both time series showed a positive increasing trend over time. In contrast, the bulk velocity at sunrise showed a weak trend over time with a mean of 1.78 ± 0.08 mm s -1 in 2015 and 1.71 ± 0.1 mm s -1 in the following year.

Temperature measurements
The daily average temperature (Fig. 5) showed that the lake was already stratified in June 2015 and the epilimnion started to deepen at the end of the same month. In 2016 the water temperature varied little with depth in early May and stratification began to form in mid-May. During the summer months, the water temperature reached 20°C.
The mean temperature of the migrating layer during the ascent phase (Fig. 6) showed a seasonally increasing trend in both years. There was also short-term temperature variability due to changing weather conditions. Finally, the 95% confidence interval (shaded area) is small with an average of 0.34°C in 2015 and 0.37°C in 2016. This indicates that temperature variations experienced by zooplankton during the DVM were little.

Dissipation rates of TKE
The mean dissipation rate in 2015 (Fig. 7a) was 10 À8 W kg -1 , while turbulence was lower in 2016 with an average of 10 À9 W kg -1 (Fig. 7b). The data also show that there is no systematic trend over time in the mean dissipation level in the water column. Chlorophyll-a concentration The concentration of chlorophyll-a was very low for both years (Fig. 8). In 2015 the concentration in the water column peaked at 2.5 lg l -1 near the lake bottom at the end of June. In 2016 the concentration showed a maximum at the beginning of the time series (early May) and it reached a deep maximum of 5 lg l -1 in July below 10 m. The time series of depthaveraged chlorophyll-a concentration (Fig. 9) showed little evidence of a temporal trend throughout the study period.

Underwater light conditions
PAR profiles (Fig. 10) showed a linear trend on a semi-logarithmic scale, with the exception of the data above 1 m, indicating an exponential decay with depth of the ambient light in the water column. Some profiles provided I PAR \50 W m -2 because they were collected near dusk. The Secchi disc depths (Fig. 11) showed a mean depth of 10.4 ± 2 m in 2015 with a decreasing trend after the beginning of July despite a constant and low chlorophyll concentration in the water column (Fig. 8). In 2016 z SD was on average 10.1 ± 1 m and never fell below 8 m. The time series of K PAR and K SD showed a good agreement (Fig. 12). The coefficient K showed a positive trend in 2015, with a mean value of 0.19 ± 0.03 m -1 , and stationary conditions in 2016 with a mean of 0.18 ± 0.01 m -1 . It never exceeded 0.22 m -1 , with the exception of K SD at the end of July.
Since K PAR and K SD were time-independent, the relative change in light intensity S has been estimated from Eq. 6. The result (Fig. 13) shows that S was always negative and in most of the cases was within the range S ¼ ½À0:1; À0:04 min -1 , suggested for inducing phototactic behaviour in Daphnia (Ringelberg, 2010).

Zooplankton concentration
The time series of the zooplankton vertical concentration in 2015 (Fig. 14a) showed that Daphnia was dominant, accounting for up to 80% of the mean concentration in the water column in June. The concentration remained almost constant until the end of the same month and then it started decreasing reaching a minimum of 7 org. l -1 in July. The concentration of adult, copepodite and naupliar stages of copepods was very small in June ( $ 5 org. l -1 ), but in mid-July they became as abundant as Daphnia with an average density of up 10 org l -1 . The maximum concentration of zooplankton (Fig. 14c) followed the same trend with a peak of 35 org. l -1 of Daphnia at the end of June and a maximum of 20 org. l -1 of copepods (including nauplii).
Data in 2016 displayed a more complete picture of the zooplankton density dynamics. From late April until June, copepods dominated with an average concentration of 27 org. l -1 (Fig. 14b). The Daphnia population increased throughout May peaking at 25 org. l -1 in July. The density of small cladocera remained very low and negligible with respect to the other species. For the remainder of the summer the total zooplankton abundance decreased, until the end of August where a maximum of 5 org. l -1 of Daphnia and 11 org. l -1 of copepods were observed. The maximum abundance in 2016 (Fig. 14d) was very similar in magnitude to that observed in 2015 with the exception of copepod nauplii whose abundance reached 42 org. l -1 in late May. The concentration declined quickly with a minimum of 12 org. l -1 at the end of the time series.
The trend in the VBS of the migrating layer at the time of the migration in 2016 (Fig. 15) is similar to that observed in the time series of the zooplankton concentration from net tows (Fig. 14c, d): VBS peaks at 81.5dB at the beginning of July and at 80 dB at the end of the same month. In 2015 the observations period of the VBS is too short to be compared with the zooplankton concentration. Local variations in VBS can also be observed which were driven by daily variations of zooplankton concentration and thickness of the migrating layer.

Correlation between v up and potential drivers
To investigate potential drivers of the observed variability in v up , the upward bulk velocity was linearly correlated with the temperature (Fig. 6), chlorophyll-a concentration (Fig. 9), the light conditions during the DVM (Fig. 13) and the VBS as a proxy for the in situ concentration and size (Fig. 15). The dissipation rates e (Fig. 7) were not used in the correlation, as turbulence values could not be interpolated because of their patchy and unsteady nature. Prior to further analysis, we removed one sampling date due to the absence of data on S in 2016.
The predictors to be included in the multiple linear regression model were tested to assess any possible multicollinearity. The correlation matrix of the variables (Table 1) shows that the VBS was negatively correlated with the temperature T (R VBS;T ¼ À 0:645) and weakly related to the chlorophyll-a concentration (R VBS;Chl ¼ À 0:375). Chlorophyll concentration was also correlated with the light availability with R S;Chl ¼ 0:649. However, the variance inflation factor (VIF) index for each of the predictors (Table 2) shows that the collinearity among the three variables was weak and unimportant (VIF VBS ¼ 2:532, VBS Chl ¼ 2:123 and VIF S ¼ 2:123\5). Variations in the chlorophyll-a concentration were in reality very small, and any small interdependence was likely a numeric artefact.
A model fitting dataset was created from the data by extracting 46 random observations (about 80% of the original sample). The remaining 12 observations constituted the validation dataset that was later used to validate the prediction of the fitted model. A multiple linear regression was applied to the 46 observations with a significance level a ¼ 0:01. The model results indicated that there was no significant statistical correlation between v up and the relative change in light intensity S (F statistic ¼ À0:927, p value ¼ 0:358), the food availability (chlorophyll-a) concentration (F ¼ À 0:2420, P ¼ 0:810), and VBS (F ¼ À 2:0328, P ¼ 0:0482). The upward velocity was, instead, significantly correlated with changes in temperature T (F ¼ 4:933 and P ¼ 1:26 Â 10 À5 ). The linear regression model was further simplified by excluding the other predictors, and the result is reported in Fig. 19a with P ¼ 2 Â 10 À12 and a coefficient of determination R 2 ¼ 0:66. The final regression model also predicted the observations from the validation dataset reasonably well (Fig. 19b) with P ¼ 1 Â 10 À4 and R 2 ¼ 0:79.

Modelling of v down
The parameters employed in Eqs. 7-10 are reported in Fig. 18 and Table 3 for the Daphnia density q d and sizes, respectively. The mean density of a single Daphnia was 1019.2 kg m -3 (red line) with CI 95 ¼ ½1016:8; 1022:5 kg m -3 (black error bar). The organism's mean body length l was 0.886 mm, while its width w was 0.378 mm, with b e ¼ 2:34, clearly indicating a deviation from the spheric shape. The antenna length b c was 0.425 mm and its width a c ¼ 0:039 mm. The dependence of the Daphnia sinking rate on the water temperature T changes according to the equation used to model the organism body shape (see Fig. 20). The model for U sink;s consistently provides larger sinking rates with respect to the other two models with large variations as the water temperature T increases. The difference between U sink;e and U sink is about 0.5 mm s -1 , however, U sink;e slightly diverges from U sink for warmer temperatures.

Discussion
We used for the first time acoustic measurements to quantify the time series of zooplankton displacement velocity during the DVM in an oligotrophic lake. Based upon environmental measurements, and independent physical calculations and modelling, we   aimed to infer the environmental drivers that might influence the rate of zooplankton migration in the field.

Upward bulk velocity
In our correlation model, we hypothesised that v up may depend on temperature, zooplankton density, food availability and light intensity. Because of lack of sufficient data, fish abundance and turbulence were not included in the correlation analysis. The effect of temperature on zooplankton swimming velocity has been demonstrated by laboratory studies of single swimming organisms (Ringelberg, 2010). When temperature increases, zooplankton metabolic activity rates change (Heinle, 1969;Gorski & Dodson, 1996;Beveridge et al., 2010), and they can also swim faster because of the reduced drag exerted on their swimming appendages (Larsen & Riisgård, 2009). Temperature may also enhance phototactic response in Daphnia ( Van Gool & Ringelberg, 2003); however, this was not observed in our results. A stepwise linear regression, with pairwise interactions of all DVM drivers, showed that the interactive effect of S and T on DV was not significant (F ¼ 1:0672, P ¼ 0:307). In principle, if zooplankton are present at very high population densities while migrating, organisms would have to reduce the velocity they can swim at to avoid encounters with other animals (i.e., the surrounding space for moving would be limited). The instantaneous speed and the resulting bulk velocity would be smaller. In the opposite scenario, with a low abundance within the migrating layer, organisms would be free to swim unobstructed, reaching the maximum thrust they can propel at, resulting in a higher bulk velocity. Our dataset showed, however, that the VBS was weakly correlated with v up probably because the zooplankton concentration was not high enough to yield such swimming behaviour. This effect might have been appreciable if measurements had been available also during the peak in zooplankton abundance in June (Fig. 14c, d). The lack of other studies in the literature precludes us from corroborating this hypothesis.
Food availability can both positively and negatively affect Daphnia migration velocity depending on environmental conditions and genetic differences (Dodson et al., 1997a). Our analysis, however, did not provide evidence of this effect because of the stationary chlorophyll-a concentration during the period over which the migration velocities changed.
Surprisingly, the bulk velocity was not affected by the average change in light intensity. Changes in light intensity have been shown to influence zooplankton displacement velocity at the scale of a single DVM (Ringelberg et al., 1991;Ringelberg & Flik, 1994). However, this relationship was not observed at a multiple-DVM scale, in our study. It is likely that this is due to among-study differences in field conditions, and in the impact of additional drivers and mediators of migratory behaviour.
The lack of data on fish abundance and kairomone concentrations did not allow direct inclusion of predation pressure in our analyses. However, some speculation can be made about its possible effect. While we observed an increasing near-linear temporal trend in v up , predation often has a seasonal peak coinciding with the time at which larval fish hatch and begin actively feeding upon zooplankton; this leads to a stronger zooplankton response to light within that specific period ( Van Gool & Ringelberg, 2003). The  Table 3. The grey box represents the observed range of the downwards velocity (v down ) in the field zooplankton layer depth prior the DVM (black dots in Fig. 16) suggests that the organisms did not change their daytime vertical position during the observation period and that they resided on average at the same depth in both years. If a strong predation pressure had occurred over a specific time period, a repositioning of the population towards deeper layers in the daytime would have been expected during this period (Dodson, 1990;Leach et al., 2014). However, the near-constant migration amplitude suggests that visual predation pressure was present and consistent throughout the whole study period, and thus a sudden increase in that pressure appears unlikely. Turbulence was not included in our model either. Available measurements could not be interpolated to fill the gaps of missing observations. However, e was constant and very low (Fig. 7), approaching the background dissipation level. The lack of a clear temporal trend in turbulence suggest that this could not be an important driver of the observed trend in zooplankton migration velocity.
The correlation model showed that temperature was a primary factor in controlling migration speed in the field with respect to the other DVM drivers. Because zooplankton swim in a low Reynolds number environment, temperature-dependent viscosity is a crucial mediator of their motility. Assuming that the viscosity mainly controls the organism swimming response (Larsen & Riisgård, 2009), our dataset showed that a viscosity variation of 13% (from m ¼ 1:8 Â 10 À6 m 2 s -1 at the start of the dataset to m ¼ 1:1 Â 10 À6 m 2 s -1 in August) led to an increase in velocity of about 61%. The same variation was observed by Larsen et al. (2008) in the laboratory, but over a wider viscosity range (about 36%) for the instantaneous velocity of the copepod Acartia tonsa Dana, 1849. Moison et al. (2012) found instead an increase of 29% in the speed of Temora longicornis Müller O.F., 1785 for a viscosity variation of 15%. Because we were dealing with a mean velocity during the migration and no further studies were available, it was not possible to make any direct comparisons with other experiments.

Downwards bulk velocity
The velocity U sink (blue line in Fig. 20) from our new model matched the observed range of v down in the field (grey box) within the 95% confidence interval (blue shaded area). Our analysis also indicates that spherical models are not suitable for reproducing Daphnia sinking rates; U sink;s (orange line) consistently provided velocities above the falling rates reported in the literature for narcotised animals (Hutchinson, 1967;Gorski & Dodson, 1996;Ringelberg & Flik, 1994;Dore et al., 2009), because the genus is not sphereshaped. Moreover, U sink;e (green line) still overestimated the observed velocity (grey box), because accounting for additional drag from the Daphnia antenna is important when modelling sinking rates (Gorski & Dodson, 1996;Dore et al., 2009). The good agreement between U sink and the field data strongly indicates that organisms sank during the DVM because organism's buoyancy and gravity are the only governing parameters of the reverse migration (Ringelberg & Flik, 1994). Since v down shows no trend over time (Fig. 4, negative values), no correlation exists between this velocity and any DVM drivers in the field. The only available study for Daphnia by Ringelberg & Flik (1994) suggests that the genus actively swims only when the light stimulus S is high. According to the meteorological station measurements, zooplankton started migrating several minutes before sunrise when the solar radiation was still zero. This may suggest that S was not a causal factor for active downwards DVM (Ringelberg & Flik, 1994). However, the distance of the station from the study site might have prevented us from measuring small light changes that could have initiated the DVM.
Our model also shows that the velocity changed very little with the temperature, as seen in Kirillin et al. (2012). This explains why v down was constant during the observations and did not correlate with T. Results from the literature cannot be directly compared with our model and observations because organisms have different lengths or density (i.e. eggs in brood chamber), and experiments were performed under different conditions of T.
Since Daphnia always sank during the observational period, our results indicate that these zooplankton may favour passive sinking over active swimming during the reverse DVM to preserve energy and avoid generating hydrodynamic disturbances detectable by predators.