**Changes in Vertical Distribution of Zooplankton under Wind-Induced Turbulence: A 36-Year Record**

#### **Mamoru Tanaka**

Atmosphere and Ocean Research Institute, The University of Tokyo, Kashiwa, Chiba 277-8564, Japan; tanaka.mamoru0@gmail.com

Received: 17 October 2019; Accepted: 22 November 2019; Published: 25 November 2019

**Abstract:** A multidecadal record of a local zooplankton community, stored in an open-access database, was analyzed with wind data to examine the impact of wind-induced turbulence on vertical distribution of zooplankton. Two major findings were made. First, the abundance of zooplankton assemblage (composed of copepods, cladocerans, etc.) in the upper layer (<10 m deep) decreased with increasing turbulence intensity, suggesting turbulence avoidance by zooplankton. Second, when focusing on each species, it was found that ambush (sit-and-wait) feeders showed statistically significant changes in response to turbulence, whereas suspension (filter) feeders did not. This is the first clear evidence that ambush feeders change vertical distribution in response to turbulence.

**Keywords:** white sea; arctic ocean; net tow; turbulence avoidance; feeding mode; National Centers for Environmental Information; European Centre for Medium-Range Weather Forecasts

#### **1. Introduction**

While sub-centimeter-sized zooplankton play important roles in the marine ecosystem [1], processes that control their spatial distribution are elusive. Microscale turbulence, a ubiquitous characteristic of the ocean environment [2], significantly affects zooplankton swimming, feeding, and escape behavior [3]. Encounter rates with prey and mates are enhanced by environmental turbulence [4], but, at the same time, turbulence obscures signs of approaching predators, increasing the risks posed by staying in highly turbulent regions [5]. A numerical physical–ecological simulation, which considered the trade-off between reproduction and predation, suggested that avoidance of high levels of turbulence is most advantageous for reproduction [6]. Indeed, turbulence avoidance by zooplankton has been observed in small tank experiments, in which turbulence intensity was controlled by an oscillating grid [7]. However, field studies that demonstrate turbulence avoidance by zooplankton are limited [8–11]. Moreover, their conclusions are based on relatively short-term campaigns (<10 days) [8–11]. In this study, the impact of turbulence on zooplankton distribution was examined using a multidecadal record of a local zooplankton community and long-term sea-surface wind data.

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

#### *2.1. Biological Parameters*

Zooplankton data were obtained from an open-access database provided by the National Centers for Environmental Information. A biological dataset "*36-Year Time Series (1963–1998) of Zooplankton, Temperature, and Salinity in the White Sea*" [12] was used in this study. Zooplankton assemblage was sampled at the White Sea Biological Station (66◦19.5 N, 33◦39.4 E; the Arctic Ocean) from 1963 to 1998, i.e., 36 years. The water depth was 65 m at the station. A standard juday net was vertically towed every 10 days, collecting samples from surface (0 to 10 m deep), middle (10 to 25 m), and deep (25 to

65 m) layers. Mesh size and mouth area were 168 μm and 0.1 m2, respectively. Zooplankton were fixed with 10% formaldehyde solution and classified to the species level. The net tows were performed during daytime. A total of 814 net tows were made.

Representative vertical position, mean depth distribution (MDD; m) was calculated for each net tow:

$$\text{MDD} = \frac{1}{N} \sum n\_i \mathbf{z}\_{i\prime} \tag{1}$$

where *N* is the total abundance of zooplankton in the water column, *ni* is the abundance in the *i*th depth bin (i.e., *i* = 1, 2, and 3), and *zi* is the average depth of the *i*th depth bin (i.e., *z*<sup>1</sup> = 5 m, *z*<sup>2</sup> = 17.5 m, and *z*<sup>3</sup> = 45 m). MDD was calculated for the entire zooplankton assemblage, as well as for each species. In case a certain species was not observed in any depth layer, MDD was not calculated for this species. MDD and zooplankton abundance in the surface layer (individuals m<sup>−</sup>2) were analyzed with turbulence intensity.

#### *2.2. Physical Parameters*

Turbulence intensity at the biological station was estimated with wind data provided by the European Centre for Medium-Range Weather Forecasts. The historical reanalysis dataset "*ERA-40*" [13] was used in this study, covering the entire period of the biological data. A sequence of wind speed at 10 m above the sea surface at the biological station was downloaded with a temporal resolution of 6 h. Wind events may need to be sustained for several hours to affect underwater distribution of zooplankton [8–11]. Hence, representative wind speeds *U* (m s<sup>−</sup>1) during the net tows were obtained by averaging wind speeds at 09:00 and 15:00 in local time (*U*<sup>09</sup> and *U*15, respectively). The wind speed, *U*, was rejected when the difference between *U*<sup>09</sup> and *U*<sup>15</sup> exceeded half of *U*. Consequently, 105 of 814 data points (about 13%) were rejected.

Turbulent kinetic energy (TKE) dissipation rate ε (W kg<sup>−</sup>1), which is a typical parameter to quantify turbulence intensity, was estimated by the "law of wall" method. This study employs an empirical equation valid for the layers shallower than 10 m deep, which was provided by MacKenzie and Leggett (1993) [14]:

$$
\log\_{10}(\varepsilon\_{\rm ML}) = a \cdot \log\_{10}(\rm lJ) + b \cdot \log\_{10}(z) + c,\tag{2}
$$

where εML is the volume-based TKE dissipation rate (W m−3) and *z* is the depth (m). The empirical coefficients were *a* = 2.688, *b* = −1.322, and *c* = −4.812, respectively [14]. As suggested by Equation (2), turbulence intensity due to wind greatly varies along the vertical coordinate. To determine the representative turbulence intensity for the surface layer (0 to 10 m), the εML values were averaged over the surface layer, i.e., <sup>1</sup> 10 10 *<sup>z</sup>*=<sup>1</sup> <sup>ε</sup>ML(*z*). Then, the averaged <sup>ε</sup>ML (W m−3) was converted to <sup>ε</sup> (W kg−1) based on the typical seawater density of 1028 kg m−<sup>3</sup> [2]. Since Equation (2) is not applicable for ice-covered periods (November to May), ε values were rejected for those periods (296 of 814 data points; 36%). Finally, 419 pairs of biological and physical parameters were analyzed.

#### *2.3. Statistical Analyses*

Long-term environmental changes, such as slow climate changes, could induce long-term trends both in the physical and biological parameters, resulting in spurious correlations between them [15]. Hence, such long-term trends were examined by linear regression models. Seasonal changes, which have the same problem, were also examined by analysis of variance (ANOVA) tests and autocorrelation analyses. The ANOVA tests were performed for the data divided into monthly intervals (i.e., June, July, August, September, and October), whereas the autocorrelation analyses were applied to monthly averaged data. As autocorrelation is applicable for equally spaced data points, the winter data (i.e., November to May), which were not used for any other analyses, were included in the autocorrelation analyses.

Then, potential effects of ε on surface abundance and MDD of the entire assemblage were examined by linear regression models and ANOVA tests. Logarithmic values of ε were used for the linear regression models. The ANOVA tests were performed for the data divided into 8 intervals of different turbulence levels, which were equally spaced in logarithmic scale. The relationships between ε vs. the surface abundance and MDD were also examined based on the datum averaged over each cycle (i.e., June to October) to eliminate potential bias due to seasonal changes. The effects of ε on each zooplankton species were examined by two-sample t-tests, where the data were divided into 2 levels of turbulence: one for low (ε < 10−<sup>7</sup> W kg<sup>−</sup>1) and the other for high (ε > 10−<sup>7</sup> W kg<sup>−</sup>1).

#### **3. Results**

Wind speed reached 10 m s−<sup>1</sup> during the analysis period. TKE dissipation rate, ε, ranged from 10−<sup>8</sup> to 10−<sup>6</sup> W kg−1. The average was <sup>ε</sup> = 2 <sup>×</sup> 10−<sup>7</sup> <sup>±</sup> 2 <sup>×</sup> 10−<sup>7</sup> W kg−<sup>1</sup> (mean <sup>±</sup> standard deviation), typical in the surface layer [2]. Zooplankton assemblage was dominated by copepod species (Table 1), typical in the White Sea [16]. Abundance of entire zooplankton assemblage was highly concentrated in the surface layer (37,000 <sup>±</sup> 46,000 individuals m−2) relative to the middle and deep layers (17,000 <sup>±</sup> 23,000 and 5000 <sup>±</sup> 10,000 individuals m–2, respectively). Average MDD of the entire assemblage was 12.9 ± 5.2 m.

**Table 1.** List of zooplankton species. Feeding mode definitions are based on the literature [17–19]. The effects of turbulence intensity on surface abundance (0 to 10 m deep) and mean depth distribution (MDD) were examined by two-sample *t*-tests, where the data were divided into two levels of turbulence: one for low (ε < 10−<sup>7</sup> W kg<sup>−</sup>1) and the other for high (ε > 10−<sup>7</sup> W kg<sup>−</sup>1). Boldface italics indicate statistical significance (*p* < 0.02).


Linear regression models showed that long-term trends were not statistically significant in ε (*p* = 0.639), surface abundance of entire assemblage (*p* = 0.324), nor MDD of the entire assemblage (*<sup>p</sup>* <sup>=</sup> 0.822) (Table 2). Yearly change in <sup>ε</sup> was <sup>−</sup><sup>4</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> W kg−<sup>1</sup> yr<sup>−</sup>1, which corresponds to an overall decrease of 10−<sup>8</sup> W kg−<sup>1</sup> for 36 years (Figure 1a). This is one order smaller than the standard deviation of ε. Similarly, overall changes were calculated as +7380 individuals m−<sup>2</sup> for the surface abundance and −0.2 m for the MDD (Figure 1b,c), much smaller than the standard deviations of those parameters. In contrast, the ANOVA tests showed that seasonal trends were significant for the surface abundance (*p* < 0.001) and the MDD (*p* < 0.001) (Table 3). The surface abundance has a peak in August, while the MDD increased from early summer to fall (Figure 2b,c). Such seasonal cycles are also suggested by the autocorrelation analyses (Figure 3b,c). Seasonal cycles in ε were found in the autocorrelation analysis (Figure 3a) but not in the ANOVA test (*p* = 0.073; Table 3; Figure 2a).

**Table 2.** Results of linear regression models to examine long-term trends in turbulent kinetic energy dissipation rate (denoted as "ε"), surface abundance of entire assemblage (0 to 10 m deep) ("Abundance"), and mean depth distribution of the entire assemblage ("MDD"). Data from Figure 1.


**Figure 1.** Time series for turbulent kinetic energy dissipation rate (ε) (**a**), surface abundance of the entire assemblage (0 to 10 m deep) (**b**), and mean depth distribution (MDD) of the entire assemblage (**c**). Filled circles denote raw data. Red lines denote linear regression models. Results of the regression models are summarized in Table 2.

**Table 3.** ANOVA table to examine the effects of month (denoted as "Month") on turbulent kinetic energy dissipation rate ("ε"), surface abundance of entire assemblage (0 to 10 m deep) ("Abundance"), and mean depth distribution ("MDD") of the entire assemblage. Data from Figure 2.


**Figure 2.** Seasonal trends in turbulent kinetic energy dissipation rate (ε) (**a**), surface abundance of the entire assemblage (0 to 10 m deep) (**b**), and mean depth distribution (MDD) of the entire assemblage (**c**). Filled circles denote raw data. Horizontal bars denote averages over month categories (i.e., June, July, August, September, and October). Error bars denote standard error. Results of the ANOVA tests are summarized in Table 3.

**Figure 3.** Results from the autocorrelation analyses on turbulent kinetic energy dissipation rate (ε) (**a**), surface abundance of the entire assemblage (0 to 10 m deep) (**b**), and mean depth distribution (MDD) of the entire assemblage (**c**). Red lines denote 99% confidence intervals.

The surface abundance of the entire assemblage decreased with increasing ε (Figure 4a). The linear regression model shows a correlation coefficient of r = −0.12 (*p* = 0.012) between the surface abundance and ε (Table 4). Additionally, the MDD increased (deepened) with increasing ε (Figure 4c). The regression model shows r = 0.14 (*p* = 0.003; Table 4) for the MDD. Such trends against ε can be clearly seen in bar graphs (Figure 4b,d). The ANOVA tests showed statistically significant differences among the different levels of ε (*p* = 0.004 for the surface abundance and *p* = 0.010 for the MDD; Table 4). The analyses for the data averaged over seasonal cycles also suggested negative and positive slopes for the surface abundance (r = −0.33) and the MDD (r = 0.29), respectively (Figure 5), consistent with those for the raw data (Figure 4a,c). However, those trends were not statistically significant (*p* = 0.052

for the surface abundance and *p* = 0.093 for the MDD; Figure 5). This is probably due to the small range of the average ε (Figure 5).

**Figure 4.** Surface abundance (0 to 10 m deep) (**a**,**b**) and mean depth distribution (MDD) (**c**,**d**) of the entire assemblage vs. turbulent kinetic energy dissipation rate (ε). The panels on the left show raw data. Each circle corresponds to a vertical net tow sample (n = 419). Solid lines denote linear regression models. The panels on the right show averages for different levels of turbulence. Error bars denote standard error. Results of the linear regression models and ANOVA tests are summarized in Table 4.

**Table 4.** Results of linear regression models and ANOVA tests to examine the effects of turbulent kinetic energy dissipation rate (denoted as "ε") on surface abundance of the entire assemblage (0 to 10 m deep) ("Abundance") and mean depth distribution ("MDD") of entire assemblage. *r* denotes the correlation coefficient. Data from Figure 4.


**Figure 5.** Same figures as in Figure 4 (left panels) but for the data averaged over seasonal cycle (i.e., June to October). Horizontal and vertical lines denote standard error. The first year, 1963, was excluded since net tows started from September. The number of samples is n = 35.

When focusing on each species, the feeding mode was found to be associated with sensitivity to turbulence. Ambush feeders, such as calanoid copepod *Oithona similis* and chaetognath *Saggita elegans*, showed a statistically significant decrease in surface abundance and increase in the MDD in response to increased ε (*p* < 0.02; Table 1). In contrast, suspension feeders, such as calanoid copepod

*Temora longicornis*, cladoceran *Evadne nordmanni*, and appendicularian *Fritillaria borealis*, showed no significant changes in surface abundance or MDD (Table 1). Calanoid copepod *Acartia longiremis*, which can switch between ambush and suspension feeding, exhibited no significant changes (Table 1). No changes were found for harpacticoid copepod *Microstella norvegica*, a particle feeder.

#### **4. Discussion**

Analysis of a multidecadal record of a local zooplankton community revealed avoidance of wind-induced turbulence by zooplankton, whereas significant response to turbulence was found only in ambush feeders (Table 1). This is consistent with laboratory experiments that demonstrated that ambush feeding is hindered by high levels of turbulence, while suspension feeding is less dependent on turbulence intensity [20,21]. Such turbulent effects on ambush feeding are also demonstrated by a theoretical model, which is designed to predict gut contents of ambush feeders in turbulent water [22]. In the literature, the model results were compared with those from field campaigns and showed that gut contents of ambush feeders decreased with increasing turbulence intensity [22]. Results from this study are consistent with those from the experimental and theoretical works [20–22].

In contrast to ambush feeders, pure suspension feeders exhibited no significant changes in response to turbulence. Given that suspension feeders are able to adapt to relatively high levels of turbulence [6], they probably place less priority on changing their position. Additionally, *Acartia longiremis*, which have multiple feeding modes, may switch their feeding mode to suspension feeding when in turbulent waters, rather than seek out low levels of turbulence (a similar discussion is seen in [21]). Although the particle feeder *Microstella norvegica* showed no significant changes in response to turbulence, this species, in another study, exhibited significant migration to deeper depths in response to wind-induced turbulence [11]. The reason for the difference between this study and the literature [11] remains unclear.

Physical processes, such as wind [23] and surface cooling [24], frequently disturb surface waters, producing vertical gradients in turbulence intensity in the water column [25]. Hence, downward migration is generally optimal behavior to seek lower levels of turbulence. However, turbulence is also generated by other processes, such as bottom stress associated with barotropic tides and swells [26], and shear stress associated with internal gravity waves [27]. Hence, actual turbulence intensities in the ocean would be different from those estimated by the simple model in Equation (2). This will result in errors in estimation of TKE dissipation rate (ε). Additionally, the reanalysis dataset would include potential errors in wind speed, which consequently induce additional errors in ε The European Centre for Medium-Range Weather Forecasts does not quantify the magnitude of the errors, but it could be substantial.

The r value of surface abundance vs. ε was r = −0.12 (Table 4), which corresponds to the coefficient of determination, r2, of 0.0144. This means only 1% of the fluctuation is explained by ε. Zooplankton generally exhibit highly intermittent distributions in the ocean, resulting in high levels of inter-sample variability in zooplankton density [28]. This means that single/instantaneous samples could have substantial (generally large) differences from the true density [28]. Such an inter-sample variability would be a source of potential errors (or biological noise) in surface abundance (Figures 1b, 2b and 4a).

Here, I compare turbulent velocity to typical swimming speeds of zooplankton. While turbulent flow speeds are highly variable in space and time, the representative velocity scale near the boundary is the friction velocity (*u*∗) (m s−1). The friction velocity (*u*∗) is a function of <sup>ε</sup> and the depth (*z*) (m) (e.g., Equation (4.9) in [2]):

$$
\mu\_\* = (\varepsilon \ z \ \kappa)^{1/3},
\tag{3}
$$

where κ is the von Karman constant (=0.41). Assuming ε = 10−<sup>7</sup> W kg−<sup>1</sup> (at which the significant changes in surface abundance and MDD were found) and *<sup>z</sup>* = 10 m, we obtain *<sup>u</sup>*<sup>∗</sup> = 0.7 cm s−1, which is on the same order as the turbulent velocity scales in the seasonal thermocline [29]. This is comparable with typical swimming speeds of sub-centimeter-sized zooplankton, but much lower than

their instantaneous escape jumps >10 cm s−<sup>1</sup> [29,30]. Hence, I suspect that they can perform oriented swimming even in highly turbulent waters to seek optimal levels of turbulence.

Despite the potential errors inherent in the open-access database, statistically significant trends were found between turbulence and ambush feeders. Turbulence estimation based on wind speed is valid for the surface layer and provides information about physical processes at scales of individual plankton, which is generally absent in biological samplings. I hope this study encourages other researchers to examine the reproducibility of the observed trends.

**Funding:** This research received no external funding.

**Acknowledgments:** The author gratefully acknowledges discussions with Amatzia Genin, Rubens M. Lopes, Gregory N. Ivey, Yoshinari Endo, J. Rudi Strickler, and Hidekatsu Yamazaki.

**Conflicts of Interest:** The author declares no conflicts of interest.

#### **References**


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