**Dynamics of** *Calanus* **Copepodite Structure during Little Auks' Breeding Seasons in Two Di**ff**erent Svalbard Locations**

#### **Kaja Balazy \*, Emilia Trudnowska and Katarzyna Błachowiak-Samołyk**

Institute of Oceanology, Polish Academy of Sciences, Powsta ´nców Warszawy 55, 81-712 Sopot, Poland **\*** Correspondence: kaja@iopan.pl; Tel.: +48-58-7311696

Received: 15 May 2019; Accepted: 3 July 2019; Published: 9 July 2019

**Abstract:** Populations dynamics of key zooplankton species in the European Arctic, *Calanus finmarchicus* and *Calanus glacialis* (hereafter defined as *Calanus*) may be sensitive to climate changes, which in turn is of great importance for higher trophic levels. The aim of this study was to investigate the complete copepodite structure and dynamics of *Calanus* populations in terms of body size, phenology and their relative role in the zooplankton community over time in different hydrographic conditions (two fjords on the West Spitsbergen Shelf, cold Hornsund vs. warm Kongsfjorden), from the perspective of their planktivorous predator, the little auk. High-resolution zooplankton measurements (taken by nets and a laser optical plankton counter) were adapted to the timing of bird's breeding in the 2015 and 2016 summer seasons, and to their maximal diving depth (≤50 m). In Hornsund, the share of the *Calanus* in zooplankton community was greater and the copepodite structure was progressively older over time, matching the little auks timing. The importance of *Calanus* was much lower in Kongsfjorden, as represented mainly by younger copepodites, presumably due to the Atlantic water advections, thus making this area a less favourable feeding ground. Our results highlight the need for further studies on the match/mismatch between *Calanus* and little auks, because the observed trend of altered age structure towards a domination of young copepodites and the body size reduction of *Calanus* associated with higher seawater temperatures may result in insufficient food availability for these seabirds in the future.

**Keywords:** population dynamics; size; match-mismatch; Spitsbergen; laser optical plankton counter

#### **1. Introduction**

*Calanus finmarchicus* and *Calanus glacialis* (hereafter *Calanus*) co-exist and dominate the mesozooplankton biomass in the shelf waters of the Arctic Ocean [1–3]. Originally, both species have different centers of distribution and thus are adapted to different environmental conditions and consequently adopt different life-history strategies [4–7]. As dynamic changes in the environmental conditions and in timing of primary production in the Arctic [8] cause high variability in the development rate and length of life cycles of the *Calanus* species [9], comprehensive phenological studies of these copepods are essential, especially in the context of their availability to planktivores in the warming Arctic. The important role of *Calanus* spp. in marine ecosystem functioning is based on the transfer of omega-3 fatty acids (long-chain PUFAs produced by marine algae), which are crucial for the growth and reproduction of all marine organisms [10]. *Calanus* spp. are full of lipids (up to 50–70% in dry mass) [11–13], and this makes them extremely nutritious for planktivores. Until recently, the amount of lipids was suggested to be species-specific, with a higher lipid content attributed to *C. glacialis* compared to the smaller *C. finmarchicus* [4,5,13]. This assumption emphasized the different role of these species in the Arctic ecosystem. However, in most of the studies the identification between these two morphologically very similar species was based on their size [13–15], while according to recent

molecular studies the size ranges of *C. glacialis* and *C. finmarchicus* overlap and are generally much broader than previously assumed [16–20]. Therefore, the problem with a correct separation of these species could lead to substantial confusion dealing with their life history, population dynamics and ecological roles. Because the lipid content was found to depend on body size rather than species, the analysis of the *Calanus* spp. size seems to be more appropriate to understanding the process of energy transfer from *Calanus* to higher trophic levels [20].

The individual size of *Calanus* spp. is tightly coupled with the ambient seawater temperature, with smaller individuals observed in warmer seawater conditions [16,21,22]. Temperature increase affects not only intra-species size variability but leads also to alteration of the zooplankton community size structure, which is considered to be more important than shifts in biomass [23,24]. The scenarios of pelagic food web modifications due to the increased seawater temperature indicate that smaller boreal species, due to faster reproduction, will have increasingly important roles in the higher latitudes [15,20,25]. A northward expansion of organisms, better adapted to warm water conditions, has already been observed [26–28]. Fluctuations in planktonic production may lead to a disturbance in interactions between predators and prey (match/mismatch) [29,30] and as a consequence may severely disrupt the functioning of the whole ecosystem [31–34]. However, new hypotheses have also emerged, suggesting that energy transfer to higher trophic levels may be more efficient than previously assumed, because of the accelerated zooplankton development [20]. Unfortunately, such promising scenarios may not apply to strictly specialized visual predators, actively selecting larger *Calanus* spp. individuals, such as the little auk [35,36], the keystone planktivorous seabird in the Arctic [37]. The little auk requires their prey to occur in high proportion in relation to other zooplankters in the seawater [14,37], because high abundances of small zooplankters may hinder the detection of their preferred large, energy-rich prey. The reproduction period of little auks is tightly matched in time with particular phase of *Calanus* development, because they actively select mainly the lipid-full fifth copepodite stage (CV) [38–43]. A recent study of *Calanus* spp. phenology in Greenland waters [44] showed that little auks select foraging grounds where availability of their main prey is matched in time with their high food demands. This might be especially important in Svalbard, where some of the world's largest colonies of these birds are located [45] and which is now threatened by the new climate state due to progressing Atlantification. Little auks are the main fertilizers of ornithogenic tundra and thus play an important role in the Arctic ecosystem [46]. Therefore, studies of the phenology of their main prey are important for better understanding threats for both marine and terrestrial ecosystems and their interactions.

Consequently, two main goals emerged in this study: (1) to test the size of zooplankton both on individual (*Calanus* copepodite stages) and community (abundance-weighted mean size) levels; and (2) to investigate dynamics of complete copepodite structure of *Calanus* in the context of food demand of the little auk during the summer chick rearing period in two different regions and summer seasons on the West Spitsbergen Shelf. We hypothesized that both the population dynamics of *Calanus* and the taxonomic composition of zooplankton vary significantly between more Atlantic and Arctic water dominated regions, which in turn provide different feeding conditions for planktivores.

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

#### *2.1. Study Area*

Research was conducted on the west coast of Spitsbergen in two fjords, representing different hydrographic regimes: Kongsfjorden, located in the north and Hornsund in the south (Figure 1). Hornsund is situated on the south-western tip of Spitsbergen and is influenced by the coastal Sørkapp Current, which carries less-saline, cold Arctic water [47,48]. Kongsfjorden is exposed to advection from the warmer Atlantic water of the West Spitsbergen Current (WSC) [48–50]. The two different currents are separated by a density gradient that forms the large frontal system along the West Spitsbergen Shelf (WSS) [51]. Thus, Kongsfjorden receives twice as much Atlantic water than Hornsund, and

this results in 1 ◦C higher water temperatures and 0.5 higher salinities compared to Hornsund [48]. However, in recent years, gradual warming has been observed in both fjords [48,52].

**Figure 1.** Zooplankton sampling stations in 2015 (**A**) and 2016 (**B**) from two fjords: Kongsfjorden (upper panel) and Hornsund (lower panel) located along the coast of Spitsbergen and sampled during three study periods each year (marked by various shapes of points). The sections of laser optical plankton counter (LOPC) surveys (**C**) are marked with yellow and pink lines. Map of the study area (**D**) with current patterns in the Spitsbergen region (simplified from Sakshaug et al. 2009).

#### *2.2. Sampling Protocol*

Zooplankton samples from the little auk foraging grounds were collected during two consecutive summer seasons (hereafter study season) of 2015 and 2016 (Figure 1 and Figure S1) in two study locations of cold Hornsund, and warm Kongsfjorden, both inside the fjords and in the adjacent areas located on the WSS (Figure 1). In Hornsund, samples were collected three times during each season, (hereafter study periods): three, 13, and seven samples in 2015 (hereafter respectively H1, H2, H3) and five, 15, and seven in 2016 (respectively H1 , H2 , H3 ). In Kongsfjorden nine, five, and six samples were collected in 2015 (hereafter respectively K1, K2, K3) and nine, eight, and nine in 2016 (K1 , K1.5 , K2 ). For exact sampling dates see Figure S1. In total, during 12 zooplankton sampling campaigns, 96 zooplankton samples were collected: 50 in the Hornsund area, and 46 in the Kongsfjorden area. Most sampling was performed directly from the R/V Oceania (IO PAN). Two sampling campaigns in the Kongsfjorden area (K1 and K2 ) were performed on board of the R/V Lance (Norwegian Polar Institute), one on board of a Zodiac boat (K3 in 2015) and one in Hornsund on board the S/Y Magnus Zaremba (H3 in 2016). Most of the samples (91) were collected using a WP2-type net (0.25 m<sup>2</sup> opening area) fitted with filtering gauze of 180 μm mesh size, which is best suited to catch all copepodite stages of *Calanus*. Five samples (K2) were collected using a WP2-type net with 60 μm mesh size. All zooplankton samples were collected from the mid-surface water layer (i.e., upper 50 m), which was arbitrarily chosen taking into consideration the little auk maximal diving depth of 50 m [53]. After collection, all samples were preserved in 4% formaldehyde solution in borax-buffered seawater and transported to laboratory for analysis.

#### *2.3. Laboratory Analyses*

Detailed laboratory analyses of each sample were performed according to a standard procedure, following the instructions given by Kwasniewski et al. [14]. First, zooplankters larger than 0.5 cm were removed from the sample, identified, and counted. Then, 2 mL subsamples were taken from each sample with a macropipette according to the sub-sampling method described by Harris et al. [54]. Subsamples were taken until at least 400 individuals were counted from a single sample. All organisms were identified to the lowest possible taxonomic level, typically species. Special focus was put on *Calanus* and their copepodite stages, which were identified according to criteria described by Kwasniewski et al. [55]. The abbreviations (CI-AF) refer to six successive copepodite stages of

*Calanus*, i.e., CI, CII, CIII, CIV and CV as the first five copepodite stages, and AF to adult females. After microscope analysis, the number of individuals in each sample was converted into abundance (ind. m<sup>−</sup>3), basing on the volume of filtered water.

Due to the similar function of both species according to the trait-based approach [20] and difficulties in distinguishing them properly, for the purposes of this work *C. glacialis* and *C. finmarchicus* have been merged into one group, hereafter *Calanus*. However, the measurements of the prosome length of all copepodite stages of *Calanus* were performed for at least 30 individuals from each copepodite stage in the subsamples. To compare the prosome length of all the *Calanus* copepodite stages, 12,800 individuals were measured.

#### *2.4. Temperature Measurements*

Temperature data for comparison with the measurements of body size of *Calanus* and zooplankton community from net samples were obtained from 52 stations, (35 in Hornsund and 17 in Kongsfjorden) in two years of study by vertical profiles using a conductivity–temperature–depth sensor (CTD)). Measurements were first all binned into 1 m depth intervals, and then averaged over the upper 50 m water column for each station.

#### *2.5. LOPC Measurements*

To measure the size and concentrations of *Calanus*-type particles in different temperatures, the continuous oscillatory profiles, from the surface to 50 m depth, were performed with a conductivity–temperature–depth sensor (CTD; SBE 911plus, Seabird Electronics Inc., Washington, DC, USA) and a laser optical plankton counter (LOPC; Brooke Ocean Technology Ltd., Dartmouth, NS, Canada). The LOPC is an in-situ sensor that provides data on the abundance and size structure of a plankton community by measuring each particle passing through a sampling tunnel of a 49 cm2 cross section. As the particle passes the sensor, the portion of blocked light is measured and recorded as a digital size and converted to the equivalent spherical diameter (ESD). The ESD is the diameter of a sphere that would represent the same cross-sectional area as the particle being measured with the use of a semi-empirical formula based on calibration with spheres of known diameters [56–58]. A *Calanus*-type group of particles, that involved only older life stages (app. CV), was selected basing on size (0.9–2.5 mm ESD) and transparency (attenuance index > 0.4) [59]. Then, the measurements (size + concentration of *Calanus*-type) were grouped according to a few seawater temperature ranges (<4 ◦C, 4–6 ◦C, and >6 ◦C) to calculate the 1d kernel density (R function geom\_density) estimates of the dominating size in particular water temperatures.

#### *2.6. Data Analyses*

Multivariate nonparametric permutational ANOVA (PERMANOVA) [60] was used to test differences in: *Calanus* prosome length of individual copepodite stages; *Calanus* copepodite structure (CI-AF, stage index, and zooplankton species composition among fixed factors of regions and study periods. Prior to the analyses, abundance data were square-root-transformed [61]. The distribution of centroids representing particular samples was illustrated with a non-metric multi-dimensional scaling (nMDS) using Bray–Curtis similarities ordinations. The calculation of the Pseudo-F and *p* values was based on 999 permutations of the residuals under a reduced model [62]. To assess the magnitude of the spatial variation at each gradient, the estimated components of variation (ECV) as a percentage of the total variation were used. The relationship between seawater temperature and the *Calanus* prosome length and the mean size of zooplankton organisms was tested with the use of Pearson linear correlation. The *Calanus* stage index was calculated as a weighted mean on the basis of relative abundance of particular life stages, with each stage given values from 1 (CI) to 6 (AF) [63], where AF stage was represented by adult females. The mean size of zooplankton organisms in the sample was calculated as a weighted mean based on a local database with detailed measurements made by the Svalbard fiords in the scope of the Dwarf project (Polish-Norwegian Research project

no. Pol-Nor/201992/93/2014). To determine zooplankton community structure, relative abundance of species/taxa was used. For each study period, a median number of individuals of a given taxa was calculated. Taxa constituting more than 5% of total zooplankton abundance were distinguished; the rest were grouped into an "others" category. To compare the prosome length frequency distribution of a copepodite stage CV of *C. finmarchicus* and *C. glacialis,* the number of all CVs in the sample (ind. m−3) was divided by the total number of measured individuals in the sample. To determine the relative abundance of measured individuals, each measurement of a given individual was multiplied by a factor: the abundance of CV (ind. m<sup>−</sup>3)/number of measured individuals. Measured individuals were classified into size classes every 0.05 mm.

#### **3. Results**

#### *3.1. Size Response of Calanus and Mesozooplankton to Di*ff*erent Seawater Temperatures*

The prosome length of *Calanus* copepodite stages based on comprehensive morphometrical analysis differed significantly between the two fjords in the case of all life stages (CI–CV), except for adults (AF) (Table 1). In general, the median prosome length of *Calanus* individuals was larger in Hornsund than in Kongsfjorden (Figure 2).

**Table 1.** Results of one-factor multivariate PERMANOVA for the prosome length of *Calanus* copepodite stages in Hornsund and Kongsfjorden. Bold means *p* < 0.05, df are degrees of freedom, MS represents means of squares, <sup>√</sup>ECV are square root of the estimated components of variance.


The prosome length of individual copepodite stages of *Calanus* significantly correlated with temperature in both fjords for CI-CV (Pearson correlation coefficient, R = −0.37, *p* < 0.001 (CI); R = −0.48, *p* < 0.001 (CII); R = −0.39, *p* < 0.001 (CIII); R = −0.38, *p* < 0.001 (CIV); R = −0.37, *p* < 0.001 (CV), Figure 3), and did not correlate for AF (Pearson correlation coefficient, R = 0.01, *p* = 0.88, Figure 3). Generally, a relatively large variation in the body size of each *Calanus* copepodite stage was observed for each temperature recorded, except for the highest values (>7 ◦C), at which only relatively small individuals were observed in Kongsfjorden.

**Figure 2.** Prosome length of *Calanus* copepodite stages: CI–adult females (AF) in Kongsfjorden (red) and Hornsund (blue). Horizontal black lines in the boxes show the median, box represents percentiles, whiskers indicate ranges, dots represent values outside the range, red arrows show statistically significant size differences between two investigated fjords, N indicates number of measured individuals.

**Figure 3.** The relation between seawater temperature and prosome length of *Calanus* copepodite stages CI–AF. Dots represent individuals measured, blue in Hornsund and red in Kongsfjorden. Trendline is marked in black, N indicates the number of measured individuals, R indicates Pearson correlation coefficient.

Additionally, a significant negative correlation between seawater temperature and mean zooplankton body size was observed (Pearson's correlation R = −0.46; *p* = 0.001; Figure S2). Larger animals (>1200 μm) were observed at lower seawater temperatures (i.e., below 4 ◦C) in Hornsund, while smaller organisms were dominating mainly in warmer (4–7.5 ◦C) Kongsfjorden.

Moreover, according to the measurements of the automatic LOPC method of the *Calanus*-type particles, individual size of older life stages of *Calanus* (CIV and CV) also differed in relation to the seawater temperature (Figure 4). In Hornsund in both years the largest *Calanus*-type particles (c.a. 1.5 mm ESD) were observed in low temperatures (<4 ◦C), while the highest densities of smallest individuals (1.0–1.2 mm ESD) were noted at higher temperatures (>6 ◦C). In Kongsfjorden, the size range of older *Calanus* differed between two years, with larger size fraction dominating in 2015 (Figure 4).

**Figure 4.** Size distribution of Calanus-type particles expressed as equivalent spherical diameter (ESD, mm) recorded by laser optical plankton counter (LOPC) in different temperature ranges in Kongsfjorden and Hornsund in 2015 and 2016.

#### *3.2. Zooplankton Community Structure*

The zooplankton community structure in Kongsfjorden was significantly different between the two years (PERMANOVA, MS = 6127.4, Pseudo-F = 16.13, *p* = 0.001) and among study periods in 2015 (PERMANOVA, MS = 1776.3, Pseudo-F = 5.22, *p* = 0.003) and in 2016 (PERMANOVA, MS = 2133.0, Pseudo-F = 4.49, *p* = 0.001). Because samples in K2 period of 2015 were collected with different net mesh size, they were excluded from the current analysis (see M and M for details). In 2015, *Calanus* was the most abundant taxon in K1, when it reached 47% of the total zooplankton abundance. Its percentage decreased over time reaching 25% in K3 (Figure 5). In 2016, the *Calanus* percentage was relatively low, and decreased from 12% in K1 and K1.5 to 5% in K2 (Figure 5). The second numerous taxon, *Oithona similis* was more significant in 2015 (30% in K1 and 46% in K3) than in 2016 (20% throughout three study periods). The share of the Copepoda nauplii was higher (maximum in K1 ) in 2016 than in 2015. Similarly as in Hornsund, the percentage of Bivalvia veligers were especially high in 2016 (17–26%) and very low in 2015 (≤5% in K3). Another important zooplankton component in Kongsfjorden was *Limacina helicina* with 16% to 31% contribution to the overall zooplankton abundances in 2016. The percentage of *Pseudocalanus* spp. was rather constant over time.

**Figure 5.** Zooplankton community structure in Kongsfjorden and Hornsund.

The zooplankton community structure in Hornsund also differed between the two years (PERMANOVA, MS = 3245.0, Pseudo-F = 5.83, *p* = 0.001) and between study periods (PERMANOVA, MS = 5222.6, Pseudo-F = 9.39, *p* = 0.001) of both years (PERMANOVA, MS = 1661.3, Pseudo-F = 2.98, *p* = 0.003). In both years, *Calanus* dominated the zooplankton community, constituting approx. 60% of total zooplankton abundance in all three study periods of 2015 and approximately 40% in 2016. The second most numerous taxon was *O. similis* with at least 20% of total zooplankton abundance throughout all the periods of both study seasons (Figure 5). Copepoda nauplii composed a significant percentage of total zooplankton abundance, especially in H1 and H1 (20%). The contribution of *Pseudocalanus* spp. was in general the most important in H2 and H2 study periods, while Bivalvia veligers comprised about 10% of all taxa in H2 and H3 in summer 2016.

#### *3.3. Calanus Copepodite Structure*

The copepodite structure of *Calanus* in Kongsfjorden clearly differed between the years (PERMANOVA, MS = 6079.7, Pseudo-F = 15.19, *p* = 0.001) and between the study periods only in 2015 (PERMANOVA post hoc, MS = 1782.3, Pseudo-F = 4.45, *p* = 0.013). *Calanus* copepodite stage structures were much more similar among study periods within one year than among corresponding study periods of the two years (Figure 6), which was indicated by higher estimated components of variance (ECV) for years other than for study periods (19.5% vs. 10.8%). The *Calanus* copepodite stage index in Kongsfjorden differed between years (PERMANOVA, MS = 698.1, Pseudo-F = 46.40, *p* = 0.001), but not between study periods (PERMANOVA, MS = 21.9, Pseudo-F = 1.46, *p* = 0.257). In general, in 2015 the median values of *Calanus* copepodite stage index were higher than in 2016, when they remained relatively low (Figure 7).

**Figure 6.** The non-metric multi-dimensional scaling (nMDS) of the *Calanus* copepodite structure in two Spitsbergen fjords: Kongsfjorden and Hornsund. Vectors indicate the direction of best correlating variables determined as a percentage of each copepodite stage. Their lengths correspond with the strength of the correlation.

**Figure 7.** *Calanus* stage index (Y axis) together with copepodite stage structure (different colours) in Kongsfjorden and Hornsund in particular study periods (X axis) in 2015 and 2016.

The copepodite structure of *Calanus* in Hornsund differed first of all between study periods (PERMANOVA, MS = 3724.2, Pseudo-F = 31.16, *p* = 0.001), but also slightly between years (PERMANOVA, MS = 475.1, Pseudo-F = 3.98, *p* = 0.026; Figure 6). The higher share of ECV was obtained for the factor of periods than for a factor of year (15.9% vs. 4.4%). The *Calanus* copepodite stage index in Hornsund also differed between years (PERMANOVA, MS = 49.5, Pseudo-F = 7.62, *p* = 0.014) and study periods (PERMANOVA, MS = 746.9, Pseudo-F = 57.56, *p* = 0.001). Differences in the copepodite stage index had a stepwise character in both years, with a slightly higher median values in all investigated periods observed in 2016 (Figure 7).

Large differences in *Calanus* copepodite structure in Kongsfjorden were only observed in 2015, when the proportion of early copepodite stages (CI–CIII) decreased from over 50% in K1 to about 10% in K3 (Figure 7). The CIV was the dominating life stage in K2, then it was outnumbered by CV in K3. The proportion of AF was very low with a slightly higher value, but not exceeding 2% in K1. In 2016,

the composition of *Calanus* community was dominated by young copepodite stages and was very similar among study periods (Figure 7). The CI–CIII constituted approximately 60% of all copepodite stages with about 20% share for each stage. *Calanus* population in Hornsund in both years changed over time (Figure 7). Although the highest importance of the youngest copepodites of *Calanus* (CI and CII) were observed early in both years, (H1 and H1 ), their contribution was almost two times higher in 2015 than in 2016. Thus, the correlation of these stages with the ordination coordinates was stronger concerning samples from 2015 (Figure 6). The highest proportion of CIII was observed in 2015 in H2 (31%) and in 2016 in both H1 and H2 (25%). The highest proportion of CIV was observed in the second and the third study periods in both years and was slightly higher in 2016. Peak occurrence of CV took place in the third period of both years (H3 and H3 ). The proportion of AF was relatively low throughout the study, but slightly higher in the first study periods of both years and correlated positively with the youngest life stages.

The *Calanus* copepodite structure differed significantly between Hornsund and Kongsfjorden in corresponding periods only in 2016 (PERMANOVA, MS = 5558.3, Pseudo-F = 18.95, *p* = 0.001). In 2015, the copepodite structure in both fjords was relatively similar in corresponding periods H3 and K3 (PERMANOVA, MS = 642.8, Pseudo-F = 3.35, *p* = 0.054) with the predominance of late stages (CIV and CV) in both regions. In turn, in 2016 in the comparable H2 and K1.5' periods, the copepodite structure differed significantly, because of the high percentage of early stages CI–CIII recorded in Kongsfjorden, and the domination of CIV in Hornsund. Also the copepodite stage index for *Calanus* differed between two fjords only in 2016 (PERMANOVA, MS = 258.5, Pseudo-F = 34.61, *p* = 0.001), while in 2015 in corresponding study periods the stage index in H3 and K3 was similar (Figure 7). In analogous periods of 2016 (H2 and K1.5 ) a higher median value of stage index was observed in Hornsund (Figure 7).

#### **4. Discussion**

The phenology of *Calanus* is a critical factor for little auks reproduction success during their breeding season, as was shown recently in a spatial perspective study in Greenland [44] and in time perspective study in Svalbard waters [64]. Therefore it is now of vital importance to study the match in time and space between the availability of older life stages of *Calanus* and little auks, because as was shown by this study, the development rate and the age structure of *Calanus* may differ significantly depending on the region, water temperature and time in the season. The issue is alarming not only because temperature warming has been shown to accelerate development of *Calanus* [65], but also because the altered phenology of many species is becoming an increasingly important problem for trophic interactions [31,64,66] and thus entire food webs. To date, disturbance in interactions between predators and prey (match/mismatch) have been observed in many groups of organisms, e.g., between fish and plankton [67,68], insects and plants [69], birds and insects [70–72] shorebirds and arthropods [73] or seabirds and zooplankton [33,74,75]. The high variability in *Calanus* development, smaller body size and lower proportion in their concentration in relation to smaller zooplankton taxa in warmer Kongsfjorden, observed in our study, are in line with predicted scenarios of pelagic system modifications in the future Arctic towards faster development and prevailing role of smaller organism size [15,20,25]. Even though shortening life cycles and body size reduction of *Calanus* are expected not to have negative consequences for top predators [20], it will probably be important to the little auk, which is dependent on the availability of large, energy-rich prey [35].

Little auk preferentially search for large, lipid-rich copepods to cover the high energy costs incurred during foraging trips and feeding underwater [76–78]. Because the lipid content is strongly related to the body size, in this work *Calanus* size rather than species affiliation was utilized as the main qualitative trait. This approach is suggested as more appropriate for understanding the process of energy transfer to higher trophic levels on a larger scale [20]. Combining two species into one *Calanus* category was also supported by the recently discussed and clearly demonstrated problem of misidentification of *C. glacialis* and *C. finmarchicus* [17,19,20]. In order to gain a broader view on the *Calanus* population characteristics, in this study at the first step individuals of both species were

identified and separated in accordance with traditional morphological classification [55], by measuring the prosome length. Results have shown that neither the specimens classified as *C. glacialis* nor *C. finmarchicus* were following a normal size distribution (Figure 8), which was in opposition to what was demonstrated for this stage for both species by molecular methods [20] (Figure 1). A right-skewed size distribution of *C. glacialis* both in Hornsund and in Kongsfjorden observed in our study is most probably caused by the fact that larger individuals classified as *C. finmarchicus* are in fact smaller individuals of *C. glacialis* [20]. This confirms that the size criterion is no more a reliable tool to accurately classify an individual for a given species [79], thus justifying our approach to combine them into one group.

**Figure 8.** Prosome length frequency distribution of copepodite stage CV of *Calanus finmarchicus* and *C. glacialis* classified according to morphological identification after Kwasniewski et al. [55] in Hornsund (n = 1798) and Kongsfjorden (n = 799).

The body size of individual copepodite stages (CI–CV) of *Calanus* differed clearly between the two investigated fjords, with smaller prosome length observed in warmer Kongsfjorden than in colder Hornsund. First of all, such a difference may be explained by expected differences in proportions of *Calanus* species (*C. glacialis* vs. *C. finmarchicus*), which is also of importance for little auks [66]. But the change in individual size within particular species has also to be taken into consideration. Such observations are especially important in the light of recent studies demonstrating a great range of size plasticity of *Calanus* [17,19,20] and considering the fact that temperature is a key factor determining body size in copepods [21,80]. In this work, the smallest prosome length of *Calanus* CI-CV was observed at highest temperatures (>7 ◦C), whereas the largest individuals were observed in seawater temperature <4 ◦C. A similar trend was observed by LOPC measurements of *Calanus*-type particles selected for the older life stages, with the predominance of the largest individuals observed at temperatures below 4 ◦C and the smallest ones observed at temperatures above 6 ◦C in Hornsund. Such a clear differentiation in size modes can be explained by the co-occurrence of two water masses carrying two different *Calanus* species [14] and/or that the younger (CIV) copepodite stage was dominating in warmer, close to surface waters, while older (CV) life stage prevailed in colder conditions. Although the fraction of older *Calanus*-type particles was characterized by similar sizes in all temperature ranges in Kongsfjorden in 2015, the dominating size of this fraction was evidently shifted towards smaller sizes in the warmer season in 2016. This could be caused by different relative roles between CIV and CV life stages, or by the methodological bias, because the location of transects in both years slightly differed. In 2015 transects were performed closer to the interior part of the fjord, while in 2016 transects were mainly located in open shelf waters. Smaller individuals representing mainly *C. finmarchicus* can dominate in open waters with prevailing Atlantic water masses [19], while inside the fjord usually both species co-exist [17]. Hence the larger size of *Calanus*-type particles observed in 2015 could have resulted from a higher share of the larger species (*C. glacialis*). Despite differences in species composition, the significant relation between the body size of *Calanus* and seawater temperature observed in this study agrees with the assumptions of the temperature-size-rule (TSR) [81,82], which states that ectotherms grow

slower, but mature at a larger body size in colder environments. The smaller size of *Calanus* in warmer temperatures observed in this study may be explained by the fact that organisms tend to be smaller in response to warming [83–85] and progressive reduction of *Calanus* body size is predicted with increasing seawater temperatures [86]. The body size of *C. glacialis* was found to vary considerably along its geographical range [19,20,83,87]. In experimental studies on *C. finmarchicus*, its prosome length also significantly decreased with increasing sea temperatures. Moreover, the largest individuals of *C. glacialis* were recorded after development in waters with temperatures not exceeding 3 ◦C [16].

In addition to individual size of *Calanus*, the proportion of *Calanus* in the overall community is very important for visual planktivores such as little auks, which need their prey to occur not only in a high concentration but also as easily visible [37]. In general, *Calanus* species are the key element of zooplankton communities in Svalbard waters, especially in terms of biomass [1,88,89], however their proportion in total zooplankton abundance is highly variable in time, space and under different hydrographic conditions [63]. In this study the proportion of *Calanus* in total zooplankton abundance was higher in Hornsund than in Kongsfjorden in both studied years, which confirmed more favourable foraging conditions for little auks in Hornsund than Kongsfjorden [77,90]. A similar predominance of *Calanus* in Hornsund was also recorded in 2007 [91]. Likewise, Trudnowska et al. [92] found higher proportions of *Calanus* in zooplankton communities in the colder Hornsund than in the Atlantic-influenced Magdalenefjorden. The lower proportion of *Calanus* in Kongsfjorden could result from the strong advection of the Atlantic Water carrying high concentrations of small copepods (e.g., *Oithona similis*) [27,93,94]. Therefore, the Atlantic water masses are typically avoided by little auks, due to high proportions of small copepods, which hinder detection of preferred prey [14,37]. The abundance of *O. similis*, which was really high in Kongsfjorden in this study, has been increasing gradually in Spitsbergen fjords since 2006 [91] as a consequence of the progressive Atlantification of these waters [27]. The increasing importance of small copepods in the zooplankton composition [89] is one of the most spectacular examples of the progressing warming that have already been documented [23].

Studies of *Calanus* development are challenging because its reproductive strategies are highly variable in time and space due to corresponding changes in environmental conditions and food supply [6,20,86,95]. In this study the development of the *Calanus* population and in consequence also its stage index, followed similar trends in Hornsund in both years. A similar gradual development of *Calanus* population, reflected by a dominance of young stages CI–CIII during the first study period and older CIV–CV stages during the third period was observed in Rijpfjorden [6]. Such similarity in the *Calanus* age structure observed in Rijpfjorden and Hornsund indicates a coincident timing of reproductive events and its synchrony with ice algae bloom in April and phytoplankton bloom in July [96]. In addition, the presence of early copepodite stages in all the studied periods (this study) might suggest continuous reproduction, or at least, the presence of more than one generation of *Calanus*, which is likely in high latitudes according to several new studies [20,65,97]. In Kongsfjorden, the trend of a gradual population development was observed in this work only in 2015. This observation coincided with the seasonal dynamic of the *Calanus* population structure emphasized for this fjord by a year-round investigation of Lischka and Hagen [94] with a higher contribution of early copepodite stages in July and a more advanced population in August. In turn in 2016, *Calanus* age structure was very similar in all three studied periods in Kongsfjorden and persisted relatively young according to low stage index. To some extent this might be caused by shorter time intervals between sampling periods in this year, but most probably it was caused by different advection impacts, according to different sea surface temperatures in the two years investigated (Figure 9). The events of advection are often associated with a transport of younger populations [55,93,95], which could explain the higher contribution of early stages in 2016 in Kongsfjorden in K1 (7–8 ◦C SST, Figure 9). This fact was also confirmed by a multiyear study conducted in the WSC region, where the copepodite structure of *C. finmarchicus* was younger during 'warmer' than 'colder' summers [27].

**Figure 9.** Sea surface temperature isoline and sea ice conditions (colour surfaces) in Svalbard waters in 2015 during K1 (**a**) and in 2016 during K1 (**b**) study periods. Star icon represents Kongsfjorden study area. Data source: Norwegian Meteorological Institute.

#### **5. Conclusions**

This study confirms the difficulty in proper *Calanus finmarchicus* and *Calanus glacialis* recognition, which, together with the similar functions of these species proved in recent studies, was an important argument for aggregating them into a single *Calanus* group. This study also provides evidence that the development rate and size structure of *Calanus* is highly variable in time, space and in relation to seawater temperature. Furthermore, according to the observations of this research, the copepodite structure is much more dynamically affected by Atlantification in Kongsfjorden than in Hornsund. Additionally, seawater temperature was confirmed to correlate negatively with both the mean size of mesozooplankton organisms and the body length of *Calanus* copepodites (CI–CV). These results support the hypothesis about shortening the life span and associated reduction of the body size of *Calanus* along with climate warming. The accelerated development of *Calanus* can cause a significant shift in time in availability of its fifth copepodite stage in the foraging grounds of their key predator, the little auks. These findings confirm the hypothesis of the possible mismatch in timing between the availability of *Calanus* CV and the little auks highest food demands and therefore highlight the necessity to continue further seasonal studies of *Calanus* phenology in Svalbard waters.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/7/1405/s1, Figure S1: Zooplankton sampling dates in two fjords (Kongsfjorden and Hornsund) along the coast of Spitsbergen. Codes for sampling in Kongsfjorden: K1, K2, K3 in 2015 and K1 , K1.5 , K2 in 2016. Codes for sampling in Hornsund: H1, H2, H3 in 2015 and H1 , H2 , H3 in 2016, Figure S2: Relationship between seawater temperature and the weighted mean size of zooplankton in communities sampled in Hornsund (blue) and Kongsfjorden (red). Each dot represents mean individual zooplankton size in a particular sample.

**Author Contributions:** K.B.-S. and K.B. were responsible for the research design. K.B. sampled zooplankton material and analysed the data. E.T. analysed LOPC data and prepared related figure. K.B. prepared drafted the text and figures and performed statistical analyses. All authors participated in discussions and editing.

**Funding:** This study received support from Polish Ministry of Science and Higher Education decision No. 3605/SEAPOP/2016/2 for international, co-funding project SEAPOP II for years 2016–2020. K.B. has been supported as a PhD student by the Centre for Polar Studies, Leading National Research Centre, Poland. E.T. was financed by Polish National Science Centre (ecoPlast No. 2017/27/B/NZ8/00631).

**Acknowledgments:** We thank Agnieszka Promi ´nska for CTD data. We acknowledge all people who were involved in extensive zooplankton sampling campaigns: Anette Wold, Anna Maria Kubiszyn, Justyna Wawrzynek, Maciej Jan Ejsmond, Wojciech Moskal, Michał Procajło and Mateusz Orma ´nczyk.

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

#### **References**


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

### *Article* **The Copepod** *Acartia tonsa* **Dana in a Microtidal Mediterranean Lagoon: History of a Successful Invasion**

#### **Elisa Camatti \*, Marco Pansera and Alessandro Bergamasco**

Consiglio Nazionale delle Ricerche, Istituto di Scienze Marine (CNR ISMAR), Arsenale Tesa 104, Castello 2737/F, 30122 Venezia, Italy; marco.pansera@ve.ismar.cnr.it (M.P.); alessandro.bergamasco@ismar.cnr.it (A.B.) **\*** Correspondence: elisa.camatti@ismar.cnr.it; Tel.: +39-041-2407-978

Received: 13 May 2019; Accepted: 5 June 2019; Published: 8 June 2019

**Abstract:** The Lagoon of Venice has been recognized as a hot spot for the introduction of nonindigenous species. Several anthropogenic factors as well as environmental stressors concurred to make this ecosystem ideal for invasion. Given the zooplankton ecological relevance related to the role in the marine trophic network, changes in the community have implications for environmental management and ecosystem services. This work aims to depict the relevant steps of the history of invasion of the copepod *Acartia tonsa* in the Venice lagoon, providing a recent picture of its distribution, mainly compared to congeneric residents. In this work, four datasets of mesozooplankton were examined. The four datasets covered a period from 1975 to 2017 and were used to investigate temporal trends as well as the changes in coexistence patterns among the *Acartia* species before and after *A. tonsa* settlement. Spatial distribution of *A. tonsa* was found to be significantly associated with temperature, phytoplankton, particulate organic carbon (POC), chlorophyll *a*, and counter gradient of salinity, confirming that *A. tonsa* is an opportunistic tolerant species. As for previously dominant species, *Paracartia latisetosa* almost disappeared, and *Acartia margalefi* was not completely excluded. In 2014–2017, *A. tonsa* was found to be the dominant *Acartia* species in the lagoon.

**Keywords:** *Acartia tonsa*; Lagoon of Venice; nonindigenous species; zooplankton distribution; coexistence patterns; niche overlaps; long-term ecological research

#### **1. Introduction**

Alien invasive species, also called non-native or nonindigenous species (NIS), together with marine pollution, overexploitation of living resources, and physical alteration of habitats, represent the main threats to the world's oceans on local, regional, and global scales [1]. The planktonic crustacean *Acartia tonsa* Dana (1849) (calanoid copepod) is an NIS recently introduced in the Mediterranean Sea [2]. Several studies have shown that *Acartia* NIS are colonizing coastal areas and estuaries by propagation or introduction (e.g., [3,4]). The ability of Acartiidae to cross geographic barriers relies mainly in their capability of producing resting stages [5]. These NIS are modifying the status of native species, which are subject to competitive pressure [6]. *A. tonsa* is widely distributed in estuarine environments along the Atlantic coasts of North and South America [7] and the Pacific coast of North America [8] where it is the most abundant species. It appeared in the European coasts in the first half of the 20th century, possibly transferred by ship ballast waters [9,10]. In the Mediterranean Sea, *A. tonsa* was first reported in 1985 in the Etang de Berre, a eutrophic lagoon near Marseilles, France [2]. Only after 1985 was the presence of *A. tonsa* confirmed in several Italian transitional waters such as in a lagoon of the Po River Delta (northern Adriatic Sea) [11], in the Lagoon of Lesina [12], and in the Lagoon of Venice [13]. In estuarine ecosystems, *A. tonsa* generally reaches relatively high abundances and becomes often the dominant zooplankton species during summer [14], provided that high concentrations of particulate

organic carbon and particulate organic matter are available. In fact, the life cycle of this copepod is strictly dependent on the quantity of the available food—the larger the trophic supplies are, the more accelerated its growth rate is [15–17].

The Lagoon of Venice, a large Mediterranean lagoon located in the northwest coast of the Adriatic Sea, presents marked habitat heterogeneity, and the classification of its habitats is still a matter of debate [18,19]. The lagoon is considered the main hotspot for invasive species in the whole Italian coast with the presence of more than 60 NIS, including 29 invertebrates and 34 macrophyte species among which *A. tonsa* also appears [20–25]. The lagoon is also part of the Long-Term Ecological Research (LTER) network (http://www.lteritalia.it), a global network of research sites located in a wide array of ecosystems. LTER research is a fundamental tool for monitoring environmental changes over time. Zooplankton communities exhibit an intrinsically high variability in transitional environments, and the elucidation of coexistence patterns is a critical question in ecology as well as in accounting for differences in abundances among species. In addition to zooplankton ecological relevance related to the role in the marine trophic network, changes in the community have implications for environmental management and ecosystem services.

This work describes relevant steps in the history of invasion and establishment of *A. tonsa* in the Lagoon of Venice with respect to local stressors. This work provides a recent picture of its distribution along gradients of environmental parameters and identifying the relevant biogeochemical quantities assisting or limiting the successful colonization of the species, mainly compared to congeneric residents in the lagoon (*Acartia margalefi* Alcaraz (1976), *Paracartia latisetosa* Krichagin (1873), and *Acartia clausi* Giesbrecht (1889)). The study will focus on abundance, distribution, and niche interactions of the Acartiidae community to address the question of niche separation and to investigate multidimensional niche breadths under a temporal and trophic gradient.

The work aims to identify which ecological factors are most important for *A. tonsa* population structure and organization and to provide a possible key to disentangle the roles of *Acartia* lagoon dominant species based on their niche characteristics. Identification or exclusion of possible overlaps among realized niches in the space of considered environmental variables will aim to highlight the specialization of each species and the relevance of coexistence mechanisms in the structuration of the community. Combining spatial and temporal frames in the three different periods taken into consideration, along a gradient from the mainland to the inlets, determining species' relative abundance distributions (RADs) in a given habitat and time will support the comparison of the reciprocal position of the species of interest within the mesozooplankton community and also how they eventually interacted in a competitive way.

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

#### *2.1. Study Area*

The Lagoon of Venice is a large Mediterranean lagoon with an area of approximately 500 km<sup>2</sup> stretching along the northwest coast of the Adriatic Sea (Figure 1). Three inlets connect the lagoon to the sea, and its general circulation results from the superposition of tide, wind, and topographic control [26,27]. The effective renewal rate of water is on the order of a few days for the areas closest to the inlets and up to a month for the innermost areas [26,28]. Its depth is very shallow, some 1.3 m on average and more than 10 m in the deepest channels connecting the three inlets with the city of Venice and with the industrial area of Marghera. The Lagoon of Venice is surrounded by urban and industrial areas. Moreover, fishery, mariculture, and tourist activities are very developed. Because of the complex interplay of the variety of simultaneously occurring stressors, the lagoon experiences high variability in most of the environmental parameters, showing high habitat heterogeneity [29,30].

**Figure 1.** Study area (Lagoon of Venice, Italy) and location of sampling stations during 1975–1980 (grey triangles, including Lido), 1997–2002 (black circles), 2003–2004 (grey squares), and 2014–2017 (black triangles, including San Giuliano, Marghera, Fusina, and Lido).

#### *2.2. Zooplankton and Environmental Variables Datasets*

Four zooplankton datasets (DSs) from surveys carried out in the Lagoon of Venice were considered in this study: the first dataset refers to the situation before the settlement of *A. tonsa* (DS1: period 1975–1980, monthly sampling) and the three after the first record *A. tonsa* in the lagoon (DS2: period 1997–2002, monthly sampling; DS3: period 2003–2004, monthly sampling; and DS4: period 2014–2017, seasonal sampling).

The first dataset (DS1), concerning data prior to the establishment of *A. tonsa*, was used for RAD analysis to compare the different relative species' abundance distributions as a function of the presence of *A. tonsa* along an expected environmental gradient from the sea to the mainland (Lido inlet, Crevan, Dese; Figure 1). The second dataset (DS2), the most evenly distributed in terms of spatial and temporal coverage, was used to investigate patterns and trends of *A. tonsa* almost immediately after its first record in an intermediate area of the lagoon, called Palude della Rosa. In this case, zooplankton samplings were collected from May 1997 to April 2002 at five stations located in the northern part of the lagoon (San Giuliano, Marghera, Fusina, Lido, and Palude della Rosa; Figure 1) characterized by a complex interplay of freshwater and marine inputs [26,28] and by anthropogenic pressure [31]. San Giuliano collects urban waste water from the town of Mestre, where phytoplankton blooms often develop [32,33]. Marghera is influenced by industrial pollution. Fusina is affected by thermal pollution from a thermo-electric power plant. Lido, in the northernmost inlet of the lagoon, is characterized by substantial lagoon shelf exchanges. Palude della Rosa is a typical lagoon environment, influenced both by freshwater and, to a lesser extent, by shelf water transported this far by tides.

The recent spatial distribution of *A. tonsa* populations in the lagoon, and its potential diverse response in areas subjected to different natural and anthropogenic influences, were investigated through the analysis of samples derived from monthly samplings, from April 2003 to March 2004, at six stations inside the lagoon and one station in the Adriatic Sea site just outside the Lido inlet (third data set DS3, Figure 1). Thus, with respect to the 1997–2002 stations, this dataset provided more coverage and representativeness of different areas in the lagoon with different natural and anthropogenic stressors, including a station on the outer shelf. Station 1M (Figure 1) was chosen to verify the presence of the species in the coastal sea, whereas the other stations were selected based on the following rationale: station 2B was influenced by freshwaters; station 7B collected urban waste; station 9B was affected by thermal pollution; station 11B was located in an intermediate area between northern and central lagoon basins; station 14B was in a site with seagrass meadow; and station 17B was located in the southern lagoon basin (Figure 1). Station 1M was not sampled in September 2003 and February 2004, and station 17B was not sampled in April 2003. All samples reported abundances of phytoplankton and mesozooplankton organisms as well as physical and chemical parameters related to water quality, in particular: temperature, salinity, dissolved oxygen, chlorophyll *a,* and particulate organic carbon (POC).

Samples related to the most recent period 2014–2017 (DS4) were used to confirm the presence and distribution of *A. tonsa* ten years after the previous survey along a gradient from the inner lagoon to the outer shelf. In this case, eight stations were sampled. The Inside station (Figure 1) was located in the inner channel inside the industrial harbour of Marghera, Port1 was in the lee of the cruise ships docks, 7M was at the yachting mooring docks, and PTF was the CNR Acqua Alta Oceanographic tower (http://www.ismar.cnr.it/infrastructures/piattaforma-acqua-alta), situated few miles offshore the lagoon in the Northern Adriatic. The remaining four stations corresponded to Marghera, Fusina, San Giuliano, and Lido sites of the period 1997–2002. All the considered 65 samples reported abundances of mesozooplankton organisms (66 present taxa) as well as physical and chemical parameters related to water quality, in particular: temperature, salinity, dissolved oxygen, chlorophyll *a*, nutrients, and Secchi disk.

Sampling procedures and analyses were performed in the same way for all datasets; thus, data from different years and projects were homogeneous and comparable with each other. Zooplankton samples were collected by a plankton horizontal sampler (200 μm mesh net size) and preserved with borax buffered formaldehyde for microscopic analysis. Taxonomic and quantitative zooplankton determinations were performed using a Zeiss stereomicroscope at the lowest possible taxonomic level (species for copepods and cladocerans). Each sample was poured into a beaker with 200 cm<sup>3</sup> of filtered seawater to allow a thorough mixing for random distribution of the organisms. At least four aliquots of the samples were analysed while the entire sample was checked against the identification of rare species [34,35]. Phytoplankton samples resulting from the DS3 were collected in 250 cm3 dark glass bottles and fixed with 10 cm<sup>3</sup> of 20% hexamethylentetramine buffered formaldehyde [36]. Counts were taken according to Utermöhl's (1958) [37] method, using an inverted microscope (Zeiss Axiovert 35) after 2 to 10 cm<sup>3</sup> subsamples had settled for 24 h [38]. Temperature, salinity, and dissolved oxygen were measured in situ using a CTD (Conductivity Temperature Depth) Idronaut Ocean Seven 316 Multiprobe. Particulate organic carbon (POC) was determined using a CHN (Carbon Hydrogen Nitrogen) analyser [39], whereas surface chlorophyll *a* (Chl *a*) was analysed with spectrophotometric methods [40].

#### *2.3. Statistical Analyses*

In order to characterize and explain the *A. tonsa* population dynamics and trends after its arrival in the lagoon, a seasonal Kendall test (SK) [41] was performed over the data set 1997–2002 (DS2). The test performed the Mann–Kendall (MK) trend test for individual seasons of the year, where season was defined by the user (here, corresponding to monthly values). It then combined the individual

results into an overall test to depict whether the dependent variable changed in a consistent direction (monotonic trend) over time or not.

To ascertain if species abundances and environmental parameters were associated, the nonparametric Kendall rank correlation test was carried out [42] on the 2003–2004 dataset (DS3). A correspondence analysis (CA) [43] was carried out on this dataset to highlight the possible relationships and groupings between stations (representative of different habitats) and zooplankton species (proxy of different water masses). CA is the most suitable statistical technique to analyse enumerative data [44]. As it is based on the chi squared metric, this algorithm automatically weights both low and high frequencies. The dataset was organized in a species⁄station matrix. Statistical analyses were performed using MATLAB® software.

Species RADs were used to support the comparison of the whole zooplankton community throughout the four decades and the different habitats covered by the available datasets. In particular, three types of habitats were identified along a gradient from the mainland to the sea: the innermost shallow zones characterized by low water circulation (inner), the inlet areas at the interface between the lagoon and the Adriatic Sea with strong marine characters (inlet), and the transition areas with less direct links to the sea (intermediate). For these three habitats, all available samples of each of the three decades (decade D1: 1975–1985, from data set DS1; decade D3: 1995–2005, from data set DS3; recent period D4: 2005–2017, from DS4) were integrated to obtain 3 × 3 snapshots (before invasion, early settlement, and more mature condition in each habitat) of the RADs of the overall zooplankton community.

Relationships between environmental variables and zooplankton community composition in the current decade (D4)—with particular focus on *A. tonsa* and the congeneric *A. margalefi*, *P. latisetosa,* and *A. clausi*—were studied using an outlying mean index (OMI) analysis through the 'ade4' package in R, extended to the overall zooplankton taxa (dataset DS4). OMI analysis [45] is a multivariate technique to perform niche analyses of species assemblages and explore the relationships between environmental gradients and community structure. The focus is on a specialization criterion called OMI index (i.e., species marginality), which measures the distance between the average habitat conditions used by each species and the average habitat conditions of the sampling area. To characterize the realized niche of each species, the analysis extracts two other terms: the tolerance index (Tol), which measures the habitat breadth of the species, and the residual tolerance (RTol), which represents the variance in the species niche not explained by the measured environmental variables.

#### **3. Results**

During the decade before the first record of *A. tonsa* in lagoon (DS1), the mesozooplankton community was composed of more than 40 taxa, with a weak decreasing gradient of overall richness from the inlet (46 taxa) to the intermediate areas (43 taxa) and to the inner zone (41 taxa). Among them, copepods were represented by 24 species, and Acartiidae largely dominated the community with *A. clausi* being the most abundant taxon in every habitat. The species *A. tonsa* was not present in any sample of DS1, even in the innermost areas, where the three dominant species (*A. clausi, A. margalefi, P. latisetosa*) represented more than 80% of the overall community.

Conversely, over the analysed 1997–2002 period (DS2), *A. tonsa* reached higher abundances at the inner lagoon stations San Giuliano (Figure 2) and Palude della Rosa (Figure 2) and constituted about 90% of the lagoon mesozooplankton community (Figure 3). Seasonal cycles (Figure 2) showed that the species was, generally, nearly absent in the cold season while annual peaks in abundances were reached in summer. The population started to rapidly increase in May, when suitable water temperatures (>15 ◦C) assisted the growth, and decreased in fall. Annual maxima were reached in different months, depending mainly on the station features. San Giuliano and Palude della Rosa, the stations with larger abundances, had annual peaks in July.

**Figure 2.** Box plot of monthly *A. tonsa* abundance from the dataset 1997–2002 for the five sampled stations (San Giuliano, Marghera, Fusina, Lido, and Palude della Rosa). The cross (μ) indicates the average, 25–75% indicates the interquartile range, and 10–90% indicates the interdecile range.

**Figure 3.** Relative abundance, station by station, of *A. tonsa* with respect to the other *Acartia* species and the remaining zooplankton during the 1997–2002 period.

SK analysis showed a statistically significant increase of the species only at the San Giuliano station (Table 1) along with total zooplankton abundance. The positive trend of total zooplankton appeared to be due to *A. tonsa* exclusively. At Palude della Rosa, the second station in order of abundance, the trend of *A. tonsa* and total zooplankton abundance was also positive, although not statistically significant (Table 1). In the other stations, where *A. tonsa* abundance was generally lower and reached the minimum at Lido station, the trend was found to be negative and not significant (Table 1).


**Table 1.** Results of the seasonal Kendall test (1997–2002 period).

Despite the positive, significant trend of *A. tonsa* at San Giuliano, the estimate of Sen's slope [46] was 30 ind (number of individuals)/m3/year. This negligible increase over time (Figure 4) may be suggestive of a steady state of the population (i.e., a mature stage of colonization). Alongside the significantly increasing trend of *A. tonsa* abundance at the San Giuliano station, SK highlighted an opposite trend for the other representative species of the *Acartia* genus historically relevant in the lagoon: *A. margalefi* and *A. clausi* (Table 1). In particular, *A. margalefi* significantly (in a statistical sense) decreased at all stations, providing additional evidence of the ongoing decline in abundance of this species.

**Figure 4.** Time series of monthly *A. tonsa* abundance (bars) along with a 12-month running mean and median and Sen's trend line at the San Giuliano station during 1997–2002.

Analyses related to the April 2003 to March 2004 period (DS3) allowed characterization of the actual distribution of *A. tonsa* over an environmental gradient in the Lagoon of Venice. The presence of a decreasing gradient of salinity can be noted, shifting from stations near inlets to stations well inside the lagoon (Table 2). Higher mean salinity values were found at station 14B, which was located near an inlet and, therefore, directly influenced by tidal exchanges, and the marine coastal station 1M (Table 2). Lower salinity values were observed at the innermost stations (2B, 7B, and 9B). The presence of *A. tonsa* at station 1M was negligible and exclusively was due to ebb tidal transport, as this species is usually not resident in marine waters. Indicators of trophic water quality—particulate organic carbon (POC), Chl *a* concentration, and phytoplankton abundance—showed an opposite gradient with respect to salinity, with the larger values at stations 2B, 7B, and 9B (Table 2). The species recorded the largest abundances in those stations (2B and 9B) where POC, Chl *a*, and phytoplankton values were high and where also the mean annual mesozooplankton abundances were highest (6386 ± 10,327 and 4892 ± 11,350 ind m−<sup>3</sup> respectively; Table 2). Lowest abundance at the site 14B corresponded to the lowest concentrations of POC and Chl *a* as well as low phytoplankton abundance (Table 2).

Percentage contributions of the most abundant taxa to total zooplankton community in each station are reported in Table 3. A few species, such as the copepods *A. tonsa*, *A. clausi*, *Paracalanus parvus* Claus (1863), and *Centropages ponticus* Karavaev (1895), were always present in the zooplankton community over the whole investigated area with a percentage contribution greater than 76%. *A. tonsa* was the most abundant in all lagoon stations with the exception of 14B and 1M stations, which were dominated by *A. clausi*. The percent contributions of the codominants *P. parvus* and *C. ponticus* were generally lower than 5% except for the 14% contribution of *P. parvus* at site 1M. The cladoceran *Penilia avirostris* Dana (1849) was present only at the marine coastal station 1M. The other taxa (i.e., decapods larvae) contributed with low abundances to zooplankton composition.


**Table 2.** Average and standard deviation values of the hydro-chemical data and biological concentrations (phytoplankton, *A. tonsa*, and main mesozooplankton

With respect to the entire zooplankton community, total abundances increased in late spring and summer at all stations except 1M where, instead, spring maxima were observed in May (17 <sup>×</sup> <sup>10</sup><sup>3</sup> ind m<sup>−</sup>3) and June (8 <sup>×</sup> 10<sup>3</sup> ind m−3) mainly because of the copepod *A. clausi* (Figure 5). Similar to the datasets 1997–2002 (DS2), the warm season corresponded to the maxima in abundances of *A. tonsa*. In general, the main mesozooplankton groups showed seasonal fluctuations with maxima during warm periods; copepods dominated throughout the year. *A. margalefi* showed very low abundances (maximum mean value of 10 ind m−<sup>3</sup> in June), while *P. latisetosa* was never found. The rest of the zooplankton community was mainly represented by the meroplankton, especially by decapod, gastropod, and cirriped larvae during the spring and summer periods (Figure 5).

**Figure 5.** Monthly abundance of *A. tonsa* with respect to other *Acartia* species and remaining zooplankton at coastal station 1M (**top panel**) and at lagoon stations (**bottom panel**) for the 2003–2004 period.

A match between abundances and distribution of *A. tonsa* and environmental parameters emerged from the Kendall rank correlation test for DS3 (Table 4), which showed that *A. tonsa* was the only species of dominant copepods positively and significantly correlated with temperature, phytoplankton, POC and Chl a concentrations, as well as *A*. *margalefi* species. However, the positive correlation of *A. margalefi* with POC was not significant. With respect to the same parameters, the congener *A. clausi* showed the opposite correlation and was not statistically significant.


**Table 3.** Mean percentage contributions of the most abundant species/taxa at the seven stations for the 2003–2004 period.

**Table 4.** Results of the Kendall rank correlation test at seven stations for the 2003–2004 period. Statistically significant values are marked in bold.


The CA was suggestive of a separation of the stations into one main group, which included lagoon stations 2B, 9B, 17B, 11B, and 7B, and the other two isolated stations located (i) in an area filled by seagrass in the central basin close to an inlet (14B), which was heavily influenced by tidal exchanges, and (ii) the coastal station 1M (Figure 6). The main group, associated mainly with inner lagoon sites, included all stations where *A. tonsa* dominated and where all trophic parameters favored its abundance. These stations had common species compositions because their species percentage distributions were similar with respect to all other stations. In fact, the same species percentage distribution considerably differed at the inlet station 14B and coastal station 1M. Indeed, Figure 7 shows that station 14B was characterized by a more heterogeneous community at the taxa level, and that higher percentages of marine taxa prevailed at station 1M such as appendicularians and cladocerans *P. avirostris*, *Evadne* spp., and *Podon polyphemoides* Leuckart (1859). The strictly neritic species *A. clausi* was found in both stations 1M and 14B, sites with high salinity and low inorganic nutrient waters, confirming the preference

for coastal waters. Station 7B, the most distant in CA space from the other lagoon stations, was characterized mainly by the presence of *A. margalefi*, fish larvae, and copepod nauplii.

**Figure 6.** Results of the correspondence analysis (CA) performed on seven stations and total annual abundances (May 2003–April 2004). The first letter of each species' name corresponds to the actual position on the CA space. For enhanced readability, only species with relative abundance >5% are shown. Stations B2 and B9 as well as Polychaeta larvae, Decapoda larvae, and Cirripeda nauplii species have been moved from their respective original positions, indicated by the black lines.

**Figure 7.** Relative abundance, station by station, of A. tonsa with respect to the other Acartia species and the remaining zooplankton during 2003–2004 period.

Data related to the last period 2014–2017 (DS4), though only seasonal, confirmed the presence and the dominance of three principal *Acartia* species inside the lagoon: *A. margalefi*, *A. clausi* and *A. tonsa*. A different distribution of these species along physical and trophic gradients was also found: *A. margalefi* coexisted with *A. tonsa* mainly in the intermediate lagoon areas (Figure 8). Both virtually disappeared at stations with marine features (Lido and PTF) where *A. clausi* instead dominated. *P. latisetosa* was found with very limited abundance in the inner lagoon stations (Inside, Marghera, and Fusina) only in May 2014 at 4 ind m−<sup>3</sup> or less.

**Figure 8.** Seasonal relative abundance of *A. tonsa* with respect to other *Acartia* species and remaining zooplankton at stations for the 2014–2017 period. Data related to 7M, Port1, and inside refer only to May 2014, October 2014, and February 2015 sampling dates.

RDA analysis showed that the relative rank distribution of the most representative zooplankton species in the community changed over time and in space by comparing the periods before and after the arrival of *A. tonsa*. First of all, Figure 9 shows that much of the RAD shape was steeply dominated by the arrival of *A. tonsa*, especially in inner areas where the species reached maximum abundance values and exceeded those of the congeneric ones. *A. clausi* maintained the most representative relative abundance of the *Acartia* genus throughout the period prior to the arrival of *A. tonsa* and in all three considered areas (Figure 9). It maintained the same relative rank also in the following periods but only in the inlet areas. In the years following its appearance, *A. tonsa* replaced *A. clausi* in the intermediate and inner areas. *P. latisetosa* reached an important distribution rank only in the internal areas and in the pre-*A. tonsa* period, after which it could always be found on the low-rank curve tail both in spatial and temporal terms. *A. margalefi*, also with an important relative abundance, before the arrival of *A. tonsa*, especially in the most typically lagoon areas (intermediate and inner), fell in rank with the appearance of *A. tonsa* (D3), but then (D4) recovered in terms of presence and abundance mainly in the intermediate areas (Figure 9).

**Figure 9.** Rank abundance distribution (RAD) and occurrence of the zooplankton community of the Venice lagoon during the period of study (decade D1: 1975–1985, from data set DS1; decade D3: 1995–2005, from data set DS3; and recent period D4: 2005–2017, from data set DS4). The samples were lumped according to the three habitat types along a gradient from the mainland to the sea (see text). In every frame, the upper curve gives the ranked occurrence of each species in the lumped N samples; the *Acartia* species's relative ranks are highlighted along the lower RAD curve.

With an OMI analysis, we computed and tested niche parameters (Table 5) to describe marginality, tolerance, and, thus, the variability of responses of *Acartia* species to environmental variables as well as their possible niche overlap. Figure 10 shows the representation of the statistically significant species' realized niche position of *Acartia* species on the first factorial plane of the OMI analysis whose origin represents the average habitat conditions of the sampling area. The two first axes of the OMI analysis accounted for 81% of the marginality (69% for the first axis). As a consequence, subsequent graphs used only these two axes. Seventeen out of 66 taxa showed a significant deviation (p < 0.05) of

their niche from the origin (whereas among the Acartiidae, only *A. clausi* was significant, p < 0.01) suggesting a significant influence of the environmental conditions for a relevant part of the community (Monte-Carlo krandtest, perm = 999). Furthermore, the between-site analysis confirmed that the environmental gradient of sites (sea, inlet, intermediate, and inner) could be discriminated through the inhabiting zooplankton community (Monte-Carlo rtest, perm = 99, p < 0.01). In our case, the tendency of low values of marginality for *Acartia* species indicated that there was no significant difference between the average overall environmental conditions and those where the species were preferentially found. The niches of the three typical lagoon species, *P. latisetosa*, *A. margalefi,* and *A. tonsa*, gravitated around less salty and more trophic environments and also presented an evident overlap. Furthermore, *A. tonsa* and *A. clausi*, the latter having more neritic coastal characteristics, were equidistant from the average environmental conditions. Therefore, this showed a clear affinity, or the opposite, depending on the environment taken into consideration: coastal marine sites for *A. clausi* (sea and inlet in the chart) and lagoon for *A. tonsa* (intermediate and inner). In particular, *A. tonsa* had a low marginality and a high residual tolerance to environmental conditions, whereas *A. clausi* showed an opposite trend but was still relative low OMI values (Table 5). As for the other two typical lagoon species, *A. margalefi* had the lowest value of OMI, while *P. latisetosa* had the lowest RTol in comparison to the rest of the *Acartia* species (Table 5). This meant that *A. clausi* was greatly influenced by environmental conditions, whereas *P. latisetosa* seemed to poorly depend on unconsidered environmental parameters (e.g., POC). Instead, *A. tonsa* showed wider potential ecological tolerance. In general, most common habitat conditions covered by the sampling units corresponded to the ubiquitous or generalist species. In contrast, specialists, which deviated from these general habitat conditions, demonstrated high OMI values, as in the case of P. *latisetosa* with respect to the other two lagoon congeners. Tol values showed a high correlation and dependence of *P. latisetosa* on environmental variables, while *A. tonsa*, with its high RTol, showed a greater adaptability to the variations of the ecosystem in which it gravitated. *A. margalefi* had the lowest OMI and Tol values.

**Table 5.** Niche parameters of the 20 most abundant zooplankton taxa in the Lagoon of Venice during 2014–2017 (outlying mean index (OMI) analysis). The inertia of each species, OMI, the tolerance index (Tol), and the residual tolerance index (RTol) are indicated. The last column (p) represents the percentage of random permutations (out of 1000) that yielded a higher value than the observed OMI (significant cases are in bold, p < 0.05).


<sup>1</sup> *P. latisetosa* does not fall within the 20 higher ranks.

**Figure 10.** OMI analysis. Realized niche position of *Acartia* species in the Lagoon (2014–2017 dataset). The site scores are projected on the same first two axes of the OMI analysis, and the habitat type is highlighted (Labels: acla = *A. clausi*; aton = *A. tonsa*; plat = *P. latisetosa*; and amar = *A. margalefi*).

#### **4. Discussion**

The occurrence of NIS is increasing in marine and estuarine systems. Among the 955 new taxa reported for the Mediterranean Sea, 42 are planktonic copepods [47]. Their invasive behaviour has been recognized as one of the major threats to the conservation of the biodiversity and the functioning of marine ecosystems [48–50]. *A. tonsa* dated its first appearance in the Venice lagoon in 1992. Since then, it never disappeared, and it is now present throughout the year and dominant with the exception of colder months. The species is widespread throughout the lagoon with significant abundance except for the areas closer to the tidal inlets, where the trophic condition is lowered and the hydro-chemical characteristics become less favourable to the species. In the Venice lagoon, *A. tonsa* is currently considered a stabilized species. The most plausible hypotheses about the introduction of this species are that it was brought into the outer port area thanks to its ability to produce resting eggs, then via ballast water from ships [51–53], or released by aquaculture, fisheries, or pet industries [54].

Our study showed how the trophic conditions of the lagoon system were (and still are) the main factors that influenced its adaptive success. Unlike the other main species present within the lagoon zooplankton communities, *A. tonsa* seems strongly dependent on trophic conditions, as it is positively correlated with nutrients, Chl *a*, POC, and phytoplankton concentrations. The favourable environmental conditions that the species found in 1992, fundamentally characterized by high habitat trophic levels [55], are congruent with our current findings and to the ecological characteristics of the species [8,15,56]. This allowed the settlement of *A. tonsa* in the Lagoon of Venice.

Much of the available knowledge about zooplankton communities of the Venice lagoon, mainly copepods, derives from studies along the physical and trophic gradient of the northern lagoon. In this area, *Acartia* genus usually represents a quantitative, important component of the zooplankton communities. It shows spatial and seasonal segregation patterns associated with the hydrological conditions, the seasonal variability, and the trophic status and pollution of the water [31,57,58].

Comparison of the whole zooplankton community throughout the four decades and the different habitats highlighted that in the Venice lagoon, before the arrival of *A. tonsa*, the genus *Acartia* dominated the copepod community with *P. latisetosa* and *A. margalefi* usually found in the innermost and intermediate parts of the lagoon, respectively. In particular, in the more internal areas, *A. margalefi* and *P. latisetosa* ranked second and third, respectively, in terms of presence and abundance, though *P. latisetosa* greatly diminished in significance in the intermediate areas. The dominant species *A. clausi* ranked always first, mainly in the coastal area. This picture instead has been profoundly modified with the passage from the first analysed decade (i.e., before the invasion of *A. tonsa*) to the following ones. In fact, the current scenario shows, again, the *Acartia* species among the dominant species of the entire lagoon zooplankton community, and it is always organized along a spatial gradient determined by salinity and trophic conditions but with modified rankings in the different considered areas. *A. clausi* no longer dominates all the considered sub-basins, except for the inlet area where it is, again, the most representative species both in space and in time. *P. latisetosa*, whose presence and dominance were already weak in the intermediate areas, further lost importance in the inner areas and virtually disappeared except in spring 2014. Very little abundance of the species was observed only in the internal channels of the lagoon (Inside, Marghera, and Fusina stations). We argue this could either be due to its segregation in more confined and not yet investigated areas or to the occupation of its original habitat by *A. tonsa*. In the presence of a large salinity and trophic range, the *Acartia* species coexist in time with a strong spatial segregation that shows *A. clausi* limited to the inlet and coastal areas, *A. margalefi* in the intermediate area, and partly overlapping with *A. tonsa* that dominates in the innermost areas.

During the second decade (1995–2005), coinciding with the first phase of stabilization of the species after its first occurrence dated 1992, *A. tonsa* seems to clearly supplant the *A. clausi* and *A. margalefi* congeners in the intermediate and mostly internal areas, whereas occupation of *A. clausi* near the lagoon mouths failed. It is interesting to note that, in more recent times (2005–2017), *A. margalefi* appears to be recovering its original ranking, returning to coexist with *A. tonsa* in the intermediate areas. As reported in other Adriatic lagoons [59,60], the arrival of *A. tonsa* may have caused the disappearance or drastic decrease of congeneric species, and the current coexistence of *A. tonsa* and *A. margalefi* is probably favored by a greater resilience of the latter as well as a return to more favorable environmental conditions always for the latter (decline in nutrient concentration and increase in ecological status) [61].

Therefore, among the main factors that seem to have favoured the adaptive success of *A. tonsa* compared to its congeners, we can ascribe the deterioration of the environmental conditions during the 1970s and the 1980s to when the lagoon of Venice was affected by abnormal inputs of pollutants and trophic loadings, mainly phosphorus and nitrogen compounds, which induced massive macroalgae growth (*Ulva rigida* C. Agardh, 1823) whose decomposition in summer led to frequent dystrophic crises [62–65]. Later, during the 1990s, the Lagoon of Venice experienced significant changes in both primary production and trophic conditions. In particular, *Ulva* coverage rapidly declined [64], mainly because of climatic changes and the increase of sediment resuspension and sedimentation [66]. As a consequence, since the end of 1980s and through the 1990s, the presence patterns of *Acartia* species also changed the mesozooplankton composition of the lagoon.

Within this scenario, and knowing that *A. tonsa* can tolerate low dissolved oxygen concentrations and adapt well to hypoxia [67,68], we hypothesized that the hypoxia and hypertrophy in the inner lagoon area could not have hindered the settlement of *A. tonsa* whose adaptive success almost seemed to have benefited. Our results seem to support the hypothesis that, after a first important collapse of *A. margalefi* coinciding with the deterioration of environmental conditions and with the advent of *A. tonsa*, the first remained segregated in the intermediate areas of the lagoon, leaving *A. tonsa* to the innermost ones. Given the niche overlap that emerged from the OMI analysis, it is quite clear that competition for the same resources may be one of the main factors affecting the coexistence of the different *Acartia* in the same area, and that the observed decline in abundance of *P. latisetosa* and *A. margalefi* could be linked to changes in environmental conditions during the 1980s that favored *A. tonsa* development. *Acartia* genus is known to have both herbivorous and omnivorous feeding habits, depending on the environmental resources [69], although little knowledge exists about the nutritional requirements of *A. margalefi* [70]. The available studies report that *P. latisetosa* probably has an omnivorous or detritivorous diet just like *A. tonsa* [71]. Our analyses suggest for these two

congeners a stricter dependence on the considered trophic variables resulting in a stronger competition for food with *A. tonsa*.

Therefore, change in distribution and size of the populations of *A. margalefi* and *P. latisetosa* occurred from the second half of the 1980s to the present. This probably depended on the success of recruitment that, in turn, could have been influenced by the system's trophic conditions but also, consequently, by the presence of effective competitors such as *A. tonsa*, which are also known to be an active predator of early larval forms. Competition with other copepods, especially congeneric, appears to be the main documented impact of *A. tonsa* (e.g., [4,51,52,72,73]) besides being known that *A. tonsa* species may coexist with indigenous Acartiidae (e.g., [6]) or exclude them by competition [5], as suggested in some Mediterranean systems [57,59,60], where it seems that *A. margalefi* and *A. clausi* are not able to settle down in large populations in areas colonized by *A. tonsa* [59,74].

This may indicate that large-scale changes, as well as a biotic interactions with the species and mainly with the "new" species *A. tonsa,* might be responsible for the evolution in the lagoon of Venice of different dominant copepod associations with respect to the past. The triggering phenomenon seems to have been the anthropic action on such a delicate and dynamic ecosystem, whose impact *A. tonsa* took advantage of.

Currently *A. tonsa* dominates the intermediate and internal areas of the lagoon, and it reaches the highest population abundance to the detriment of the indigenous congener species. The process, initially favored by the negative effect of changes in environmental conditions, worsened because a clear niche overlap emerged, which was confirmed by our study, in particular between *A. tonsa* and *A. margalefi* with *A. tonsa* benefiting from the point of view of trophic resources. The extent of its area of distribution in the lagoon confirms its euryhaline species characteristics and its excellent adaptability to instantaneous variations in habitat conditions. At present, the lagoon is oligo-mesotrophic over a great part of its surface [61] and dominated by areas classified as polyhaline (18 to 30 salinity) that correspond to the optimal development of *A. tonsa* (salinity ranging from 15 to 22), where the lowest salinities limit its egg production [75].

The adaptive success of *A. tonsa* with respect to its congener in the Lagoon of Venice corroborate the findings that *A. tonsa* is an opportunistic, tolerant species that can take advantage of eutrophic/impacted ecosystems [4,9,31,76,77] and invade them.

Restriction of *A. tonsa* distribution offshore seems mainly influenced by food availability and less by the salinity. In the Lagoon of Venice, *A. tonsa* can sustain the high energy requirements. Moreover, competition with true marine copepods is reduced because of the habitat selection that takes place along the lagoon–sea gradient. It is also probable that the spatial width and heterogeneity of the Venetian lagoon favored and allowed the gradual recovery of the community of *A. margalefi* in recent years, which we observed in the present work. Succumbing to the highly competitive level of the alien *A. tonsa*, *A. margalefi* seems to have rediscovered a niche that does not completely overlap that of *A. tonsa* right along the saline and trophic gradient that exists in the lagoon.

The reports of new NIS in the Venice lagoon—the copepods *Pseudiaptomus marinus* Sato (1913), *Oithona davisae* Ferrari F.D. and Orsi (1984), and ctenophore *Mnemiopsis leidyi* A. Agassiz (1865) [53,78]—are also very recent. The latter, reported for the first time in the lagoon in 2016, is still present throughout the year, and its strong predator characteristics on the planktonic component [78] could revolutionize again the structure of the *Acartia* genus and in general of the zooplankton communities in a few years. Copepods are important components of the pelagic food web, as they are themselves food resources for numerous bentho-pelagic invertebrates and planktivorous juveniles of fishes. The advent of new species capable of competing for the same resources, as *M. leidyi*, could trigger a process at the base of the trophic levels that would affect the highest ones. So, the story of an invasion could therefore not have ended here.

#### **5. Conclusions**

The Lagoon of Venice has been recognized as a hot spot for the introduction of nonindigenous species [25]. Several anthropogenic factors (industrial and urban pollution, mariculture, shipping, and tourism) as well as environmental stressors (e.g., warming) concurred to make this ecosystem ideal for invasion as in other estuaries and coastal areas of the Mediterranean. In this work, four datasets of mesozooplankton were examined, with particular emphasis on Acartiidae and the alien species *A. tonsa*. The first dataset, dated 1975–1980, was used as a term of comparison before *A. tonsa* settlement. The second dataset (during 1997–2002) was used to investigate the seasonal cycle and temporal trends. The third dataset (2003–2004) was used to underpin the spatial variability of abundance and species composition, and the last dataset (the recent period until 2017) was used to confirm the presence of *A. tonsa* several years later. The selected trend test on *A. tonsa* abundance did not point out significant overall increases in the period 1997–2002, but it did point out increases in one single station (the most eutrophic one). This may be suggestive that the population has reached a mature stage of colonization in the lagoon. The annual cycle showed maxima in the warm season (mostly July) and negligible abundances in the cold season. Spatial distribution of *A. tonsa* was found to be significantly, positively associated with temperature, phytoplankton, POC, chlorophyll *a*, and counter gradient of salinity, which confirmed that *A. tonsa* was an opportunistic tolerant species that could take advantage of eutrophic ecosystems. This would seem to have occurred at the expense of the previously dominant species *A. margalefi* and *P. latisetosa,* which clearly declined in abundance over the last years. While *P. latisetosa* almost disappeared, *A. margalefi* was not completely excluded. Stations inside the lagoon showed similar species compositions remarkably different from the station in the outer shelf, which was more representative of coastal conditions and dominated by *A. clausi*. In 2014–2017, *A. tonsa* was found to be still one of the dominant *Acartia* species in the lagoon.

Zooplankton is known to be particularly sensitive to environmental changes, whether resulting from natural or anthropogenic forcing [79–82]. The interaction between anthropogenic activity, climate change, and plankton communities, focusing on systematic changes in plankton community structure, abundance, distribution, and phenology over recent decades, is a key global issue as well as the potential socioeconomic impacts of these changes [83]. In environments of relentless evolution such as the Lagoon of Venice, monitoring the ecological dynamics of species is of the utmost importance. In particular, continuous modifications in zooplankton assemblages in response to anthropogenic and environmental stressors must be considered. Ongoing plankton monitoring in the LTER site of the Lagoon of Venice will act as sentinel research to identify future changes in this complex, heavily impacted ecosystem and its related food web.

**Author Contributions:** E.C. contributed substantially to the study's conceptualization, data acquisition, analysis, and to the original draft preparation; M.P. contributed to recent data acquisition and analysis and to the review and editing; and A.B. contributed substantially to statistical analyses and to the review of the manuscript. All authors approved the final submitted manuscript.

**Funding:** This work was supported by the LIFE-WATERS Program, EC-sponsored, for the two-year period 1997–1999. The data set 2002–2003 was collected within the framework of the "Integrazione delle Conoscenze del Sistema Ecologico Lagunare (ICSEL)" project, promoted by the Venice Water Authority (Magistrato alle Acque) through Servizio Ambiente—Consorzio Venezia Nuova, while the dataset 2014–2017 derived from the ADRIATIC IPA BALMAS (Ballast Water Management System for Adriatic Sea Protection) Project and from the LTER.

**Acknowledgments:** The authors acknowledge Thetis S.p.A. for providing the environmental data during 2003–2004 period and the CNR colleague F. Acri for his help in carrying out biochemical analyses. We wish to thank F. Bernardi Aubry, M. Bastianini, S. Finotto, A. Schröder, S. Pasqual, L. Dametto, and M. Casula for the fieldwork technical support during sampling and the crew of the M/B "Litus" G. Zennaro, D. Penzo, and M. Penzo.

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

#### **References**

1. IMO (International Maritime Organization). Global Ballast Water Management Project: The Problem. Available online: http://globallast.imo.org/index.asp?page=problem.htm&menu=true (accessed on 7 May 2019).


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

### *Article* **E**ff**ects of Asian Dust and Phosphorus Input on Abundance and Trophic Structure of Protists in the Southern Yellow Sea**

**Xi Chen 1, Guang-Xing Liu 1,2, Xiao Huang 1,3, Hong-Ju Chen 1,2, Chao Zhang 1,2 and Yang-Guo Zhao 1,2,\***


Received: 16 April 2019; Accepted: 3 June 2019; Published: 7 June 2019

**Abstract:** To reveal the effects of Asian dust and phosphorus (P) input on the structure and function of micro-food web in the Yellow Sea, an experiment was conducted onboard the southern Yellow Sea where P was deficient. The response of the abundance and trophic structure of planktonic protists to different concentrations of dust and P were studied. The results showed that the sand-dust deposition presented variable effects on different sizes of protists during incubation periods. At the initial stage of incubation with dust, the amount of all sizes of autotrophic protists, especially 10–20 μm, were improved; on the contrary, the heterotrophic and mixotrophic protists were inhibited. At the late period, the increase of autotrophic protists was restricted, while the 2–5 μm heterotrophic and mixotrophic protists obviously increased. Similarly, adding P demonstrated the obviously positive effect on the 10–20 μm autotrophic protists at the initial period, and then the growth was restricted at the late period. These results were consistent with that of sand-dust deposition. Hence, it could be presumed that the positive effect of sand-dust deposition on autotrophic protists in the southern Yellow Sea was achieved by the release of P from dust. P in the early stage of sand-dust deposition promotes the growth of large-size autotrophic protists, which may accelerate red tides in eutrophic ocean. The stimulation of small-size heterotrophic protists at the late period of sand-dust deposition contributed to the material cycle and food transmission in the ocean. Therefore, the effects of sand-dust deposition on the abundance and trophic structure of different sizes of planktonic protists could change the structure of micro-food web in the southern Yellow Sea and further affected the ecological function of planktonic protists.

**Keywords:** Yellow Sea; sand-dust deposition; protists; trophic structure

#### **1. Introduction**

The Yellow Sea, as one of the four marginal seas in China, is an area where land, ocean and atmosphere interact more intensely with concentrated human activities and marine economic development. It owns abundant natural resources and developed coastal economy. Previous studies found that the Yellow Sea was intermittently restricted by phosphorus (P) [1,2], which leads to a constant changing process of protists trophic structure in this ocean area, and then altered the ecological functions of protists, such as nutrients utilization and transformation [3,4]. Sand-dust deposition is an important pathway for transporting land-based nutrients and pollutants to the ocean and providing nutrients for the marine planktonic protists. Researchers found that sand-dust deposition held a significant positive correlation with chlorophyll *a* and primary productivity [5,6]. Asian dust is an important part of global dust, and the Yellow Sea, located in the downwind zone of the Asian dust source area, is the greatest probability of being affected by Asian dust in China's offshore waters. Some studies showed that Asian dust deposited into the Yellow Sea [6–8] through a long-distance transportation in sand-dust weather, and obvious "fertilization" phenomenon was observed to affect the primary productivity. Thus, Asian sand-dust deposition is an important factor that affects the primary productivity of the Yellow Sea. The nutrients carried by sand-dust [9–12] could alleviate the phosphorus deficiency in the Yellow Sea, thereby affected the abundance and trophic structure of marine biota. However, the effect of sand-dust on the growth of nanoplanktonic protists in the P-limited Yellow Sea has not been reported.

Marine nanoplankton with size ranging 2–20 μm, is a key component of marine micro-food web and they play an irreplaceable role in maintaining primary productivity and material cycle [13–15]. It is remarkable that different size groups of nanoplankton (e.g., 2–4, 4–5, 5–7, 7–10 and >10 μm) revealed different roles in the food web due to the variations of species and proportion [16]. As such, three size groups, i.e., 2–5, 5–10, and 10–20 μm of nanoplankton were proposed in this study. According to the mechanisms of energy and nutrient acquisition, nanoplanktonic protists are divided into autotrophic, heterotrophic (i.e., protozoa) and mixotrophic protists [2,17–19]. Autotrophic nanoplanktonic protists are the key contributor to marine primary productivity [17], Heterotrophic nanoplanktonic protists affect the community structure and function of nanoplankton by preying on bacteria, cyanobacteria or smaller protists [18,19], and being preyed by medium-sized zooplankton [20,21], and then flow to a higher trophic level from the bottom of the food web [22,23]. Mixotrophic nanoplanktonic protists own both the autotrophic and heterotrophic mode, and improve the utilization of nutrients by changing nutritional habits [2,3,17]. Hence, the material conversion and energy flow of the marine plankton ecosystem depends on different trophic marine nanoplankton protists. Therefore, it is of great significance to investigate the effects of dust and phosphorus input on the abundance and trophic structure of planktonic protists.

However, there are several blind spots that need to be further revealed in the relationship between the sand-dust and planktonic protists. Firstly, the effect of sand-dust input on nanoplanktonic protists is unknown. Though the sand-dust can provide nutrients during the transport process, it also adsorbs substantial heavy metals such as copper, cadmium, plumbum, and other land-based pollutants, which present a strong toxic effect on plankton [24,25]. Secondly, whether the sand-dust, as a supplement of nutrients, can effectively compensate for P limitation in the Yellow Sea is still unclear. Thirdly, it is difficult to accurately distinguish the ecological functions of different species due to the complexity of marine nanoplanktonic protists.

Hence, in this study, to make a better understanding of the effects of sand-dust and phosphorus on the abundance and trophic structure of different sizes of planktonic protists, an experiment was conducted onboard in the southern Yellow Sea where P was deficient., the corresponding changes of nanoplanktonic protists were following investigated according to the size groups (i.e., 2–5, 5–10, and 10–20 μm) and trophic types (i.e., autotrophic, heterotrophic and mixotrophic protists) [3,18,22,26–30]. The results could provide a scientific supplement for revealing the effects of sand-dust deposition on the trophic structure and ecological function of marine planktonic protists. Meanwhile, this study supply a supplement for further analyzing the effects of sand-dust deposition on marine ecosystems.

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

#### *2.1. Sand-Dust Preparation*

In this study, sand-dust samples were collected from the surface soil of Hunshandake Sandy Land (42◦22 28 N, 112◦58 34 E) in May, 2011. The samples were stored at −20 ◦C until use. The large stone in sand-dust samples was eliminated by a sieve with 20 μm of pore in diameter and then aged artificially as the previous method [31]. To analyze the concentrations of nutrients, DOC and trace metals in the sand-dust, 30 min of ultra-sonication treatment for the sand-dust was conducted under the temperature of about 0 ◦C. The contents of nutrients were analyzed by an ion chromatograph (ICS-1100, Dionex Corporation, Bannockburn, IL, USA). Heavy metals were analyzed by an inductively coupled plasma mass spectrometry (ICP-MS-7500c, Agilent Technologies, Palo Alto, CA, USA). Dissolved organic carbon (DOC) was determined by the high temperature combustion oxidation method (TOC-V, Shimadzu Corporation, Kyoto, Japan). The contents of nutrients, DOC and dissolved trace elements in sand-dust are shown in Table 1.

**Table 1.** Concentrations of nutrients, Dissolved organic carbon (DOC), dissolved trace metals in the sand-dust.


*2.2. Onboard Culturing Experiment*

In November 2014, an onboard artificial sand-dust deposition experiment was carried out on the 'Dongfanghong NO.2' scientific research ship at H03 station (36◦06 00 N, 121◦39 00 E) in the southern Yellow Sea. The location of sampling station H03 is shown in Figure 1. Surface seawater was collected by using a shipborne CTD water collector (USA). The collected seawater was immediately filtered by a silk membrane with pore of 20 μm, to remove micro, small, medium and large plankton. The filtered seawater was filled into 1.5 L sterile polyethylene terephthalate (PET) culture bottles and afterwards the bottles were placed in a water tank. The tank was tightly fixed on the deck without shading. The surface seawater from the station H03 was fed into the outer layer tank to maintain the inner bottle temperature and ensure it was the same as the surface water.

**Figure 1.** Incubation Station.

The culture experiment was comprised of five groups, i.e., control group, low dust group (LD), high dust group (HD), low phosphorus group (LP), and high phosphorus group (HP). The amounts of added dust or P (as NaH2PO4) referred to the observation data on the Yellow Sea as the previous report [32]. The amount of added dust and P in each experiment group is shown in Table 2. Dust or

NaH2PO4 were added only one time immediately after the establishment of the culture system. All culture experiments were conducted triplicate. Samples were taken every day to determine the protists, DOC and inorganic nutrients during the five days' culture period.



#### *2.3. Sample Collection and Measurement*

#### 2.3.1. Nanoplanktonic Protists Samples

Nanoplanktonic protists samples were collected at 8 a.m. each day during the incubation period. Before sampling, PET culture bottles were slowly inverted for three times to uniformly mix nanoplankton in the incubation system. There was 10 mL seawater obtained from the culture bottle and then transferred into the sterile freezing tube. Immediately, 10% (*w*/*v*) paraformaldehyde was added into the tube to fix protists (for paraformaldehyde, the final concentration was 0.5% (*w*/*v*)). After slightly mixing, the tubes were quickly frozen in a liquid nitrogen tank and then transferred to the refrigerator at −80 ◦C for storage. Three parallel samples were taken at one time.

The abundances of different sizes of nanoplanktonic protists were counted by a flow cytometry (BD C6 plus) at different wavelength of fluorescence, as explained by Christaki et al. (2011) [16] and Zubkov et al. (2007) [33]. The sample was stained by SYBR Green I, then determined by flow cytometry with excitation wavelength 488 nm. The trophic mode of the protists was determined by green fluorescence (FL1, 530 ± 20 nm) and red fluorescence (FL3, >630 nm), the protists which showed high FL1 value and low FL3 value were heterotrophic, the protists with high FL1 value and high FL3 value were autotrophic, the rest of protists were mixotrophic. The size of protists were measured by calibration beads (2, 5 and 10 μm).

#### 2.3.2. Dissolved Organic Carbon (DOC)

The seawater samples were filtered with a high-temperature treated GF/F filter (Whatman, Maidstone, UK) and stored in a high-temperature treated glass bottle at 4 ◦C in the dark. The samples were determined by using a total organic carbon meter (TOC-VCPN) [34].

#### 2.3.3. Inorganic Nutrients

The seawater samples were filtered by acid treated 0.45 μm acetic acid fibre membrane and frozen at −20 ◦C after adding the fixative. The concentrations of different nutrients were determined by QuAAtro nutrient automatic analyzer. NH4 <sup>+</sup> was oxidized by sodium hypochlorite with indigo-phenol blue (660 nm), NO3 − reduced by copper-cadmium column with naphthalene ethylenediamine hydrochloride (550 nm), NO2 − with naphthalene ethylenediamine hydrochloride (550 nm) and PO4 <sup>3</sup><sup>−</sup> with phosphomolybdenum blue (880 nm) [32,34].

#### *2.4. Data Processing and Analysis*

Multivariate ANOVA in SPSS Inc. (Chicago, IL, USA) was used to analyze the significant difference of different groups in the culture system.

#### **3. Results and Analysis**

At H03 Station, the concentration of DOC is 79.68 μmol/L, the concentration of DIN is 2.71 μmol/L, and the concentration of PO4 <sup>3</sup><sup>−</sup> is 0.13 μmol/L. The N:P ration is 21 (>16) showed that H03 station was P-limited.


Changes of 10–20μm protists under sand-dust and phosphorus addition are shown in Figure 2.

**Figure 2.** Variations in abundance of 10–20 μm protists.

The changes of 10–20 μm of protists in different concentrations of sand-dust and phosphorus supplementation groups were different from the control. The maximum values of autotrophic and heterotrophic protists in control group were 20 and 119 cells mL<sup>−</sup>1, respectively, appeared on the 2nd day of incubation. The maximum values of the mixotrophic protists were 223 cells mL−<sup>1</sup> on the first day of incubation. The peak of HP treated autotrophic protists was 43 cells mL−1, appeared on the 2nd day of incubation, which was significantly higher than that of control (*p* < 0.05). The peaks of HP treated heterotrophic and mixotrophic protists appeared on the 1st day of incubation, 393 and 571 cells mL<sup>−</sup>1, respectively, which were 3.30 and 2.56 times higher than that of control (*p* < 0.05). The high value lasted for one day and then declined rapidly. The maximum values of LP treatment for autotrophic, heterotrophic and mixotrophic protists were 26, 167, and 336 cells mL−<sup>1</sup> on the first day of

incubation, respectively, which were 1.30, 1.40 and 1.51 times of the control (*p* < 0.05). At the end of incubation, the abundance of mixotrophic protists was 2.98 times higher than that of control (*p* < 0.05).

In summary, the abundances of 10–20 μm protists in P supplying groups were higher than the control, and the maximum value of HP treated autotrophic protists was the largest. It showed that the addition of P will promote the growth of all planktonic protists, especially when the ocean area was P-limited. In the early stage of incubation, the higher concentration of P led to a stronger growth promoting effect on autotrophic protists; they increased and lasted for a longer time. At the end of incubation, the abundance of all kinds of protists in P-supplemented groups was lower than that in the early stage, and some were even lower than that at the control group. This might be caused by the more rapid consumption of nutrients at the early stage of incubation, and the nutritional deficiency that appeared at the late stage. Therefore, one-time addition of P could only promote the growth of protists in a short time. If nutrients are not supplied continuously, their growth will be limited, due to lack of nutrients in the later period. Compared with that, the abundance of LP treated mixotrophic protists was higher than that at control group at the late stage, which indicated that the P utilization efficiency of mixotrophic protists in this condition might be higher.

The differences were observed between the three trophic types of protists at the sand-dust addition groups. The abundance of autotrophic protists at HD and LD groups increased first then gradually decreased. The maximum values were 23 and 17 cells mL−<sup>1</sup> on the 2nd day after culture. The maximum value of HD treatment group was 1.15 times higher than the control (<0.05). Although the peak value of HD treatment group was lower than that of P treated group, the value on the last day of incubation was significantly higher than that of P treated group (<0.05). The results showed that the addition of high concentration sand-dust could promote the growth of 10–20 μm autotrophic protists, which might be caused by the continuously supplement P from the sand-dust. The abundances of heterotrophic and mixotrophic protists decreased first and then increased at the sand-dust adding groups. At the initial stage of incubation, they were significantly lower than that of the control and P adding groups (*p* < 0.05). At the late stage of incubation, the abundance of heterotrophic protists in the two sand-dust groups was higher than that at control and P adding groups. The abundance of mixotrophic protists was higher than that at control and HP groups, but lower than that at LP group. The abundances of heterotrophic and mixotrophic protists in HD treatment group were always higher than that of LD treatment group. Hence, the early sand-dust deposition presented an inhibiting effect on the growth of 10–20 μm heterotrophic and mixotrophic protists, which was ascribed to the toxic effect of the heavy metals and other harmful substances dissolved from the sand-dust [5]. The inhibition effect of sand-dust was much greater than the promotion of P dissolution. At the late stage of incubation, the sand-dust stimulated the growth of heterotrophic and mixotrophic protists, especially heterotrophic protists. This indicated that the heterotrophic protists gradually tolerated the harmful substances from sand-dust. Furthermore, the nutrients dissolved from the sand-dust including trace elements such as Fe [8,10,11] (Table 1), supplemented the nutrients requirement at the late stage of incubation. The effect of dust was more significant for 10–20 μm heterotrophic and mixotrophic protists at the late stage of culture.

#### 3.1.2. Response of 5–10 μm of Protists to Sand-Dust and Phosphorus Addition

Changes in the abundance of 5–10 μm of nanoprotists in different trophic modes under the stress of sand and phosphorus addition are shown in Figure 3.

**Figure 3.** Variation of protists abundance in 5–10 μm size.

The variation trend of the abundance for 5–10 μm protists during incubation is basically the same as that of 10–20 μm protists, while the variation range is quite different. For the control group, the peak values of autotrophic, heterotrophic and mixotrophic protists appeared on the 2nd day of incubation with the cell concentrations of 116, 615 and 930 cells mL−1, respectively. Compared with that, at HP treatment group, the peak values of autotrophic, heterotrophic and mixotrophic protists were 395, 2245 and 1732 cells mL<sup>−</sup>1, which were 3.41, 3.85 and 1.86 times higher than those of control (*p* < 0.05). The LP treatment groups followed the peaks of HP treatment. The maximum concentrations of autotrophic, heterotrophic and mixotrophic protists were 176, 980, and 989 cells mL−1, respectively, higher than those at control groups (*p* < 0.05), but they were significantly lower than those at HP treatment group (*p* < 0.05). The abundance decrease of protists in P-supplemented groups at the late stage of incubation might be related to the nutrient deficiency due to their rapid growth at the early stage. As such, it can be referred that the growth of 5–10 μm planktonic protists in this ocean is also limited by P. Addition of P can promote the reproduction of all of the planktonic protists. The higher the concentration of P, the stronger the promotion effect on reproduction of planktonic protists. Compared with responding results of 10–20 μm protists to P addition, it showed that the response of 5–10 μm protists to P addition is faster, but the duration is shorter.

Similarly, there were also differences between three trophic types of 5–10 μm protists in the sand-dust supplementation groups during the cultivation period. The abundance of autotrophic protists in HD and LD treatment groups increased first and then decreased. Same as P-adding group, the largest value of HD and LD groups appeared on the 1st day of incubation, with the 150 and 123 cells mL<sup>−</sup>1, respectively. The peak value of autotrophic protists abundance in HD group was 1.29 times higher than that at control group all the time. Nevertheless it was lower than that of P-adding group in the early stage of incubation (*p* < 0.05), and higher than that of P-adding group in the final stage. It showed that the addition of high concentration of sand-dust stimulated the growth of 5–10 μm autotrophic protists for a longer time, which might be ascribed to the release of P from sand-dust [6]. The heterotrophic and mixotrophic protists abundance of dust addition groups decreased at the early stage of incubation, which was significantly lower than that at control group and P-adding groups (*p* < 0.05). The higher concentration of dust resulted in a sharper decline, which was confirmed by the fact that heterotrophic and mixotrophic protists abundance of HD treatment groups were always lower than of LD group during all the culture period. Interestingly, at the late stage of incubation, the abundance of two trophic modes protists was higher than that at P-adding groups and control group. Hence, the growth of heterotrophic and mixotrophic protists was significantly inhibited by the sand-dust deposition in the early stage. The higher concentration of dust presented the stronger inhibition effect, which was related to the toxic effect of heavy metals and other harmful substances dissolved from the dust [5]. At the late stage, the growth of protists was promoted, especially for heterotrophs. This might be related to the increase of tolerance for protists to harmful substances and the supplementation of nutrients dissolved from the dust.

#### 3.1.3. Response of 2–5 μm of Protists to Sand-Dust and Phosphorus Addition

Figure 4 shows that the variation trend of 2–5 μm protists abundance is similar to that of above mentioned two sizes of protists, but the variation ranges are quite different. For the control group, the maximum values of autotrophic, heterotrophic and mixotrophic protists appeared on the 2nd day of culture, which were 4071, 9491 and 11,965 cells mL<sup>−</sup>1, respectively. However, the peak values of P supplementation groups mainly emerged on the 1st day of culture. The highest values of HP treated autotrophic, heterotrophic and mixotrophic protists were 9303, 18,867 and 23,501 cells mL<sup>−</sup>1, which were significantly higher than those of control (*p* < 0.05). The maximum values of LP treated three trophic modes protists were 6245, 11618 and 15,240 cells mL<sup>−</sup>1, respectively, higher than those of control (*p* < 0.05), but lower than those at HP group (*p* < 0.05). It can be seen that the growth of 2–5 μm planktonic protists in this ocean area was also restricted by P. The addition of P could promote the growth of planktonic protists of different trophic modes. Comparably, the response of 2–5 μm planktonic protists to P addition was faster than that of 10–20 μm protists, and the duration lasts longer on HP treated heterotrophic and mixotrophic protists.

The changes of 2–5 μm protists abundance of different trophic modes in the sand-dust addition groups were different. The highest abundance of HD and LD treated autotrophic protists appeared on the 1st day of culture, same as that of P-adding groups, with the maximum values of 4860 and 4255 cells mL<sup>−</sup>1, respectively. The peak value of HD group was significantly higher than that at control group, but lower than that at P-adding group (*p* < 0.05). At the late stage of incubation, HD group was higher than those at control and P-adding groups, indicating that the addition of high concentration sand-dust benefits to the growth of 2–5 μm autotrophic protists. This result might be related to the dissolution of P from sand-dust [6]. For the heterotrophic and mixotrophic protists, the concentration decreased at the early stage of incubation, which was significantly lower than that at control group and P-adding groups (*p* < 0.05). At the end of cultivation, the abundance of heterotrophic and mixotrophic protists increased, and got higher than those at control and the P-adding groups, especially heterotrophic protists in the LD group was 2.98 times higher than at control group. The result showed that the growth of 2–5 μm heterotrophic and mixotrophic protists, especially heterotrophic protists, was also inhibited by the sand-dust at the initial of incubation, then been promoted in the late stage.

**Figure 4.** Variation of protists abundance in 2–5 μm size.

#### *3.2. E*ff*ects of Sand-Dust and Phosphorus Addition on the Composition of Protists*

3.2.1. Effects of Dust and Phosphorus Dosage on Trophic Structure of 10–20 μm Protists

The changes of trophic structure of 10–20 μm protists under the addition of dust and phosphorus are shown in Figure 5. The mixotrophic protists predominated in each group, accounting for 63.33% on average, followed by heterotrophic protists, accounting for 30.68% on average. On the 1st day of culture, the peak value of heterotrophic protists was appeared at P-supplemented groups and the growth rate was higher than that at other groups. On the 1st day of HP addition culture, the proportion of heterotrophic protists was higher than that at control group with 10.23%, and the proportion of mixotrophic protists was lower than that of the control group with 10.30% (*p* < 0.05), the proportion of autotrophic protists did not change obviously. At the end of the culture, the proportion of autotrophic protists at HP group increased by 8.02% (*p* < 0.05) because its high growth rate maintained for a long time, and the proportion of mixotrophic protists decreased by 7.41% (*p* < 0.05). The proportions of heterotrophic protists had no significant difference among different experimental groups (*p* > 0.05). The results showed that the autotrophic protists over competed heterotrophic protists and mixotrophic protists under the P dosage condition and thus the proportion of autotrophic protists increased significantly.

**Figure 5.** Variation in the trophic structure of protists in 10–20 μm size.

The trophic structure of 10–20 μm protists in sand-dust addition groups was different from that of P addition groups. Compared with the control group, the proportion of autotrophic protists increased at the initial stage of cultivation, while the growth of heterotrophic and mixotrophic protists was inhibited, and the higher the concentration of sand-dust lead to the stronger effect. On the 1st day, the proportion of autotrophic protists at HD group increased by 8.58% (*p* < 0.05), while the proportion of heterotrophic protists and mixotrophic decreased compared with that of control group (*p* < 0.05). However, at the end of culture, the growth rate of autotrophic protists decreased due to the significant increase of the heterotrophic growth rate at sand-dust group, which made the proportion of heterotrophic protists in sand-dust group higher than that at control group by 17.86% (*p* < 0.05), and significantly higher than that of P addition, and the proportion of autotrophic and mixotrophic protists decreased by 7.80% and 10.06%, respectively, compared with the control group. The results showed that the sand-dust could stimulate the growth of autotrophic protists at early stage, but inhibit the heterotrophic protists significantly. At the late stage, sand-dust presented a stronger promoting effect on the heterotrophic protists.

#### 3.2.2. Effects of Dust and Phosphorus Dosage on Trophic Structure of 5–10 μm Protists

The response of trophic structure of 5–10 μm protists to dust and phosphorus dosage is shown in Figure 6. It showed that mixotrophic protists and heterotrophic protists were the dominant species, accounting for 45.74% and 44.88% of the total abundance. The relative abundance of autotrophic protists was relatively small, accounting for 9.38% on average. Compared with 10–20 μm groups, the proportion of 5–10 μm mixotrophic protists decreased by 17.59% (*p* < 0.05), while the proportion of heterotrophic protists increased by 14.20% (*p* < 0.05). On the first day of incubation, the growth of heterotrophic protists was higher, the proportion of mixotrophic protists reduced. At HP treatment group, heterotrophic protists increased by 12.16% (*p* < 0.05), mixotrophic protists decreased by 14.34% (*p* < 0.05), and autotrophic protists decreased slightly (*p* > 0.05), compared with control group. At the end of culture, at HD addition group, the proportion of heterotrophic protists increased by 20.47% (*p* < 0.05), the proportion of mixotrophic protists decreased by 20.58% (*p* < 0.05), and the proportion of heterotrophic protists presented little change (*p* > 0.05), compared with the control groups, respectively. Hence, the 5–10 μm heterotrophic protists were more sensitive to P addition in the early stage of culture and more competitive than the mixed and autotrophic protists. However, at the end of culture, the autotrophic protists became more competitive and their proportion increased significantly, indicating that P dosage had a stronger and longer effect on the growth of autotrophic protists.

**Figure 6.** Variation in the trophic structure of protists in 5–10 μm size.

The trophic structure of 5–10 μm protists under sand-dust addition was different from that of P addition. The change of autotrophic protists in sand-dust adding group was the same as that of the P adding group. On the first day of culture, the autotrophic protists increased significantly and the higher concentration sand-dust led to the higher scale increase. At the initial stage of culture, the proportion of autotrophic protists at LD and HD groups increased to 16.10% and 21.65% respectively, which was 9.25% and 14.79% higher than that at control group (*p* < 0.05), but the growth of heterotrophic and mixotrophic protists were significantly inhibited (*p* < 0.05). At the end of culture, the proportion of heterotrophic protists at LD and HD groups increased to 53.50% and 53.76% respectively, 37.71% and 37.97% higher than that at control group (*p* < 0.05). These increases were significantly higher than that at P-supplemented groups (*p* < 0.05). Meanwhile, the proportion of autotrophic and mixotrophic protists decreased, especially for the HD and LD treatments, they were 34.91% and 27.91% lower than those at control group (*p* < 0.05). The results showed that the early stage of sand-dust dosage had obvious promoting effect on the growth of autotrophic protists, but inhibited the heterotrophic protists. At the late stage, it presented stronger promoting effect on the heterotrophic protists.

#### 3.2.3. Effects of Dust and Phosphorus Dosage on Trophic Structure of 2–5 μm Protists

Changes in the composition and trophic structure of 2–5 μm protists under sand-dust and phosphorus addition stress are shown in Figure 7. The results showed that the mixotrophic protists dominated in 2–5 μm protists, accounting for 45.11% of the total abundance, the heterotrophic followed by 37.00%, and the autotrophic was about 17.89%. Compared with the 10–20 μm and 5–10 μm protists, the autotrophic increased by 11.89% and 8.50%, respectively, meaning which reflected that the trophic structure of 2–5 μm protists were more balanced than other sizes. On the 1st day, the protists trophic structure in P adding groups did not change obviously compared with control group, but at the end of culture, the proportion of HP and LP treated mixotrophic protists increased by 10.22% and 11.43% respectively. While the autotrophic protists decreased by 12.55% and 13.41%, respectively. The results showed that P addition presented little effect on the trophic structure of protists at the early stage, but decreased the proportion of the autotrophic and increased the proportion of the mixotrophic at the late stage.

The trophic structure of 2–5 μm protists under sand-dust adding changed largely. At the early stage, the effect of sand-dust addition on the autotrophic presented the same trend as that of P addition compared with the control. Both P and sand-dust adding improved the proportion of the

autotrophic protists. The relative abundance of autotrophic protists at LD and HD treated groups were 28.15% and 29.26%, respectively, which increased by 8.23% and 9.15% comparing with the control. However, the growth of heterotrophic and mixotrophic protists was inhibited. At the end of incubation, the proportion of heterotrophic protists increased to 52.04% and 45.94% at LD and HD treatment groups, respectively, 18.80% and 12.70% higher than that at control group. The proportion of the autotrophic protists at LD and HD treatment groups decreased to 4.40% and 14.49%, respectively, which was 19.05% and 8.97% lower than that at control group. The results showed that the sand-dust deposition had a positive effect on the growth of 2–5 μm autotrophic protists at the early stage, but increased the proportion of the heterotrophic protists at the late stage.

**Figure 7.** Variation in the trophic structure of protists in 2–5 μm size.

#### **4. Discussion**

According to the changes of relative abundance and trophic structure of protists by the rank sum test of Kruskal-wallis and Nemenyi-Wilcoxon-Wilcox, it showed that P promoted the growth of protists with the different particle sizes and the trophic modes at the early stage of culture. While the sand-dust presented different effects on the growth of different trophic modes. It stimulated the autotrophic but inhibited the heterotrophic and mixotrophic protists. The order of relative abundance of autotrophic protists with different particle sizes was 10–20 μm > 5–10 μm > 2–5 μm. The order of decreasing the proportion of heterotrophic protists with different particle sizes was 5–10 μm > 10–20 μm > 2–5 μm; and the order of inhibiting the proportion for mixotrophic protists with different particle sizes was 10–20 μm > 5 μm > 2–5 μm. The early stage of sand-dust deposition obviously inhibited the heterotrophic and mixotrophic protists by the dissolution of toxic and harmful substances in sand-dust. Mixotrophic and heterotrophic protists are the primary predator of autotrophic protists; Pearce et al. (2011) found 42%~82% primary production was consumed by mixotrophic and heterotrophic protists [35]. The decrease of mixotrophic and heterotrophic protists could accelerate the growth of autotrophic protists. In this study, sand-dust addition promoted the growth of autotrophic protists at the early stage, and the decrease of mixotrophic and heterotrophic protists accelerated this process. The further study confirmed that sand-dust storm could significantly induce the occurrence of red tide in the southern Yellow Sea [18,30]. Besides, the harmful inhibiting effect of sand-dust on heterotrophic protists community would weaken the marine matter cycle and food transfer efficiency.

At the late stage of culture, P presented a prohibitive effect on the growth of protists, which was related to the rapid growth of protists at the earlier stage, which resulted in nutritional deficiency and growth restriction at the late stage. In the late stage of sand-dust deposition, the autotrophic protists were limited due to the excessive growth of autotrophic protists in the earlier stage. At the late stage, sand-dust demonstrated stimulating effect for heterotrophic and mixotrophic protists, and the order was heterotrophic > mixotrophic protists. As for the heterotrophic protists with different particle sizes, the order was 5–10 μm > 2–5 μm > 10–20 μm, and for mixotrophic protists, the order was 10–20 μm > 5–10 μm > 2–5 μm. In the late stage of sand-dust deposition, the toxicity of harmful substances to heterotrophic protists gradually weakened. These protists slowly adapted to the sand-dust environment and dissolution of organic substances and other nutrients in sand-dust provided abundant nutrients for the heterotrophic protists [17]. As such, their abundance maintained stable and their proportion gradually increased. The effects of low sand-dust group and high sand-dust group on the proportion of heterotrophic and mixotrophic protists were basically coincident. However, the effect of high sand-dust group was stronger and the nutritional supplementation lasted longer in the later period.

For the trophic structure, at the early stage of incubation, P addition could promote all kinds of phytoplankton protists, especially the 10–20 μm autotrophic protists, which were consistent with that of autotrophic protists at the early stage of sand-dust deposition. At the late stage of incubation, both P addition and sand-dust deposition restricted all sizes of autotrophic protists. Therefore, it is speculated that the effect of sand-dust deposition on the autotrophic protists in this ocean were ascribed to the dissolution of P from sand-dust. The early promotion of sand-dust was to supplement the P in the sea area, while at the late stage, the inhibition for the autotrophic protists was mainly related to the early rapid growth of autotrophic protists by consuming up the P and other nutrients. The limiting capacity of the sand-dust group at the late stage is less than that at P-adding groups. This reflected that the dust could continuously supply the P and other nutrients.

The abundant P at the early stage of sand-dust deposition promoted the rapid growth of large-size autotrophic protists. This phenomenon will further accelerate the occurrence of red tides in eutrophic sea areas. Sand-dust deposition in the late stage stimulated the small-size heterotrophic protists and accelerated the material cycle efficiency and food transfer capacity in the sea. Therefore, the influence of sand-dust on the structure of different particle size and trophic protists will change the structure of the micro-food web in the ocean. It will also change ability of material transformation in the water body and ultimately, and affect the ecological function of protists in the transforming matter and producing food in the sea.

#### **5. Conclusions**

(1) Sand-dust deposition affected the trophic structure of different particle sizes of planktonic protists in the southern Yellow Sea. This could lead to change the structure of micro-food webs in the sea, and affect the ecological functions of micro-food webs in material transformation and food production.

(2) The growth of planktonic protists of all trophic modes in this ocean was restricted by P. The early addition of P could promote the growth of planktonic protists of all trophic modes in the southern Yellow Sea. The effect on the10–20 μm autotrophic protists was most obviously, while the late addition of P mainly restricted the growth of different sizes of protists.

(3) The effect of initial sand-dust deposition on autotrophic protists was the same as that of P, it inhibited heterotrophic and mixotrophic protists. The positive effect of sand-dust deposition on heterotrophic and mixotrophic protists was strong at the late stage, and it improved the abundance of small-sized heterotrophic protists.

(4) The positive effect of sand-dust deposition on autotrophic protists in the Yellow Sea might be related to the dissolution of P from the sand-dust. The promotion of small-size heterotrophic protists in the late stage of sand-dust deposition could accelerate the material circulation efficiency and food transformation in the sea.

**Author Contributions:** Conceptualization, X.C. and Y.-G.Z.; methodology, X.C., G.-X.L. and H.-J.C.; investigation, X.C., H.-J.C., and C.Z.; formal analysis, X.C. and C.Z.; project administration; G.-X.L. and H.-J.C.; resources, G.-X.L. and H.-J.C.; software, X.C.; supervision, G.-X.L.; writing—original draft preparation, X.C., and Y.-G.Z.; writing—review and editing, X.C., Y.-G.Z. and X.H.; funding acquisition, G.-X.L.; validation, G.-X.L.; visualization, X.C and X.H.

**Funding:** Financial support for this project was provided by the National Natural Science Foundation of China (NSFC; 41210008) and the Major State Basic Research Development Program of China (973 Program; 2014CB953701). This study was conducted at the Ocean University of China.

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

#### **References**


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

### *Article* **Stable Isotope Analysis and Persistent Organic Pollutants in Crustacean Zooplankton: The Role of Size and Seasonality**

#### **Roberta Piscia 1,\*, Michela Mazzoni 1,2, Roberta Bettinetti 1,2, Rossana Caroni 1, Davide Cicala <sup>1</sup> and Marina Manca <sup>1</sup>**


Received: 12 June 2019; Accepted: 14 July 2019; Published: 18 July 2019

**Abstract:** Zooplankton is crucial for the transfer of matter, energy, and pollutants through aquatic food webs. Primary and secondary consumers contribute to the abundance and standing stock biomass, which both vary seasonally. By means of taxa- and size-specific carbon and nitrogen stable isotope analysis, the path of pollutants through zooplankton is traced and seasonal changes are addressed, in an effort to understand pollutant dynamics in the pelagic food web. We analyzed zooplankton plurennial changes in concentration of polychlorinated biphenyls (PCBs) and dichlorodiphenyltrichloroethane and its relatives (DDTs) and in taxa-specific <sup>δ</sup>15N signatures in two size fractions, <sup>≥</sup>450 <sup>μ</sup>m and ≥850 μm, representative of the major part of zooplankton standing stock biomass and of the fraction to which fish predation is mainly directed, respectively. Our work is aimed at verifying: (1) A link between nitrogen isotopic signatures and pollutant concentrations; (2) the predominance of size versus seasonality for concentration of pollutants; and (3) the contribution of secondary versus primary consumers to carbon and nitrogen isotopic signatures. We found a prevalence of seasonality versus size in pollutant concentrations and isotopic signatures. The taxa-specific δ15N results correlated to pollutant concentrations, by means of taxa contribution to standing stock biomass and δ15N isotopic signatures. This is a step forward to understanding the taxa-specific role in pollutant transfer to planktivores and of zooplankton enrichment in PCBs and DDTs.

**Keywords:** stable isotope analysis; persistent organic pollutants; crustacean zooplankton; freshwater; size fractions; seasonality

#### **1. Introduction**

In lakes, crustacean zooplankton are low-order consumers and represent an important link between the base of the pelagic food web and the organisms at higher trophic levels, which may have economic or conservation value. Although considered a homogeneous compartment, zooplankton is composed of organisms that differ substantially from each other not only in their taxonomy, but also in body size, metabolic rate, and ecological roles [1]. While heterogeneity of the zooplankton community has been fairly well-documented in basic ecological studies, it is often overlooked in ecotoxicological ones, particularly in models, which are usually focused mainly on sources and top-levels of pollution patterns, the latter being directly related to human health. In particular, it is a crucial step in the flow of pollutants through aquatic ecosystems [2–5]. Both these pathways are involved, because zooplankton is directly related to primary producers and is able to process organic matter to feed on phytoplankton and to actively contribute to the microbial loop [6,7]. Zooplankton also represents a food source for all fish types at the larval stage and for zooplanktivorous fish at the adult stage [8,9]. Predation by fish and vertebrates is, in general, visual, therefore depending on prey size and visibility [10]. Both primary and secondary consumers may contribute to standing stock biomass in each size fraction considered. The contribution of primary and secondary consumers to zooplankton standing stock biomass generally varies with the season in deep, temperate, subalpine lakes, such as Lake Maggiore. In deep, temperate lakes, earlier spring warming and conditions of thermal water stratification influenced by climate change promote an early spring peak in phytoplankton biomass, followed by a peak in zooplankton biomass, mainly due to primary consumers, after about one month. The subsequent slow decline of zooplankton primary consumers is concomitant with an increasing contribution of zooplankton secondary consumers to the total zooplankton biomass in late summer. This pattern of change in zooplankton community composition also includes a change in organism size; in fact, a prevalence of larger sizes is observed when the contribution of predatory zooplankters increases.

Lake Maggiore has a long history of severe pollution brought about by persistent organic pollutants, such as dichlorodiphenyltrichloroethane and its metabolites. The International Commission for the Protection of the Italian–Swiss Common Waters (CIPAIS) supports monitoring activities in this lake, with a particular focus on priority substances, such as DDx and other organic chemicals of industrial origin, such as polychlorinated biphenyls (PCBs) and polybrominated diphenyl ethers (PBDEs), known for their harmful effects on human health and on environment [11,12]. Changes in physical and chemical conditions may cause their partial release into the dissolved phase, with increased bioavailability for organisms [13]. Because of their physico-chemical properties, persistent organic pollutants (POPs) can be detected everywhere, including worldwide remote regions far away from emission sources, such as polar or Alpine ecosystems [14–17]. In this study, we analyzed the seasonality of zooplankton in relation to persistent organic pollutant concentrations (i.e., DDTs and PCBs) in two different size fractions, one representative of most zooplanktonic biomass (≥450 μm) and the other representative of fish food source (≥850 μm). The former is representative of the major part of the crustacean zooplankton, while the latter represents the size fraction to which fish predation is mainly directed [10,18]. Despite its importance, zooplankton body size is usually neglected in studies on trophic pollutant transfer. Both seasonality and size are expected to drive zooplankton pollutant concentration over the season and a proper biomass estimation is essential to determine whether there is a prevalence of either of the two drivers. Furthermore, in order to study the modifications in the food web, we used carbon and nitrogen stable isotope analysis (13C, 15N; SIA). By this method, food sources (e.g., littoral versus pelagic, deep versus shallow) can be clearly identified in a system, as can the time-specific contribution of organisms to food webs. The method was originally applied to marine environments and to fish in particular, not least for fingerprinting their origin. Only recently it has been increasingly applied in freshwater environments, not only regarding fish or other organisms of direct human consumption, but also regarding those which sustain their production and growth, thus allowing for a more reliable reconstruction of the food web [19]. SIA is applied to investigate ecosystems in the making, but also to predict, by means of organism nitrogen isotopic signature, timeand space-specific trophic position, from which, among others, biomagnification depends [20–22]. The trophic role and food sources of organisms can be evaluated by carbon and nitrogen stable isotope analysis [23–25]. Isotopic signatures provide information on the resources exploited by a consumer within an environment and its trophic positions in the food web [24,26]. In particular, the isotopic carbon signature is almost constant between food source and consumer, whereas a nitrogen isotopic enrichment is observed with an increase of the trophic status [27]. In the present work, we related the trophic role of zooplankton with persistent organic pollutant (POPs) concentrations in Lake Maggiore (northern Italy). Previous long term studies in this deep, large Italian lake suggested that zooplankton seasonal dynamics and the different primary and secondary consumer contributions to zooplankton biomass are reflected by different contributions of the two important native predatory cladocerans (*Bythotrephes* and *Leptodora*), of diaptomids and cyclopoid copepods, and of the cladoceran *Daphnia*.

Our purpose is to investigate, in particular, how DDTs and PCBs vary with the season and in the two different size fractions, the former representing the bulk of crustacean zooplankton biomass and the latter representing the size fraction on which fish predation is mainly oriented. We hypothesize that variations are related to changes in δ15N signature of the two fractions and that the latter results from the taxa-specific contribution to zooplankton biomass. We expect to find a correlation between taxa-specific signatures and taxa contribution to measured zooplankton standing stock biomass.

#### **2. Study Site**

Lake Maggiore is the second deepest (Depthmax = 370 m) and largest (area = 212 km2, volume 37.5 km3) subalpine lake located in the northwestern part of Italy (Figure 1). By its nature oligotrophic, the lake underwent eutrophication during the 1960s and 1970s. The eutrophication reversal and the return to oligotrophic conditions were due to a reduction in total phosphorus discharge into the lake after wastewater treatment facilities were opened and thanks to the ban of phosphorus-containing detergents.

**Figure 1.** Map of Lake Maggiore. Sampling station location is shown with the red dot (from Google Earth).

Since 1996, contamination by dichlorodiphenyltrichloroethane and its relatives (DDTs), produced by a farm plant located in the lake catchment basin was detected [28]. Polychlorinated biphenyls (PCBs) in a steady-state condition were also reported [29–32].

Lake Maggiore is one of the best-studied lakes in Italy, with the first ecological research, including studies on zooplankton population dynamics, dating back to the early 1900s [33,34]. Long-term studies on zooplankton seasonal dynamics have been conducted over the last ten years by C and N stable isotope analyses, aimed at identifying trophic interactions and taxa-specific roles in the transfer of matter and pollutants within the pelagic food web.

#### **3. Materials and Methods**

Zooplankton sampling was performed seasonally at Ghiffa station (45◦57' N; 8◦38 E; Figure 1) in correspondence with the lake's deepest area (370 m). We collected two samples every day using two zooplankton nets that were 58 cm in diameter (450 μm and 850 μm in mesh size), which were hauled from 0 to 50 m deep several times in order to obtain a sufficient amount of organisms needed to perform the analyses (min = 2, max = 15). The total volume filtered was approximately 26 m3 per zooplankton collection. Samples of the size fraction ≥450 μm were collected from January 2012 to January 2016, while samplings of the size fraction ≥850 μm started later (May 2013). All samples were collected in duplicate. Live organisms of one of the two samples were concentrated in a laboratory setting in approximately 1 L of lake water, frozen at <sup>−</sup><sup>20</sup> ◦C, and subsequently used for carbon (13C) and nitrogen (15N) stable isotope analysis (SIA). One third of the second zooplankton sample was immediately fixed and preserved in 95% ethanol for microscopic analysis at 6.3×. We calculated the biomass of the different taxa by way of length–weight regression equations [35,36]. On the basis of the taxa, zooplankton organisms were grouped into primary consumers (herbivores, i.e., *Daphnia*, *Eubosmina*, *Diaphanosoma,* and copepod diaptomids) and secondary consumers (predators, i.e., *Bythotrephes*, *Leptodora,* and copepod cyclopoids). The remaining portion of the second sample (two-thirds) was filtered on a 1.2 μm pore glass–fiber-filter (GF/C, 4.7 cm of diameter) and then frozen at −20 ◦C. As a consequence of using large net mesh sizes, rotifers and early developmental stages of crustacean zooplankton (i.e., the smallest body size organisms) and large phytoplankton colonies were not included in our samples. However, they are of only marginal importance for the purpose of our study, in which large zooplankton was regarded as an important link in the transfer of pollutants from primary producers to fish owing to selective predation on visible prey, particularly large-bodied zooplankton [37]. SIA was performed on pooled samples of the two size fractions (≥450 μm and ≥850 μm) and on the already-listed crustacean taxa present in sufficient amounts in the sample. Under the dissecting microscope, we pooled individuals of each taxon in order to reach the minimum dry weight (DW) of 2 mg per sample (approximately 70–700 individuals depending on the biomass). Samples were oven-dried for 24 h at 60 ◦C, before being homogenized and transferred into tin capsules of 5 × 9 mm in size. The isotopic compositions of organic carbon and nitrogen were determined by Ján Veizer Stable Isotope Laboratory (Ottawa University, Ontario, Canada) following methods already described in Visconti and Manca [38]. Isotope ratios were expressed as parts per thousand (%) difference from a standard reference of PeeDee Belemnite for carbon and atmospheric N2 for nitrogen, according to the equation:

$$^{13}\text{C}\_{\prime}\text{ }^{15}\text{N}=\left[\left(\frac{\text{R}\_{\text{sample}}}{\text{R}\_{\text{standard}}}\right)-1\right]\times 1000\tag{1}$$

where R is the isotopic ratio 13C/ 12C and 15N/ 14N.

For determination of organic compounds (OCs), the materials and methods described in Mazzoni et al. [1] were followed. Quantitative DDT homologue analyses were performed by the external standard method using a solution containing pp'DDT, pp'DDE, and pp'DDD as the reference standard, prepared from single pure compounds (Pestanal, Sigma–Aldrich, Steinheim, Germany) in iso-octane (Carlo Erba, pesticide analysis grade, Val de Reuil, France). Arochlor 1260 (Alltech, IL, USA), while the addition of PCB 28, 52, and 118, was used for PCB quantification. The analyzed congeners consisted of PCB 28, 52, 101, 118, 138, 149, 153, 170, and 180. The detection limit for each OC was 1 ng·g−<sup>1</sup> lipid weight (l.w.).

We used the software SigmaPlot 11.0 (version 11, Systat Software Inc., San Jose, CA, USA) to perform statistical analyses. After verifying that the datasets were normally distributed, we applied two-way ANOVA in order to identify which factor, size or seasonality, drove changes throughout the year.

#### **4. Results**

A strong seasonality characterized stable isotope data and zooplankton biomass composition. In Figures 2 and 3, we reported data of taxa/groups biomass percentage composition, δ13C% and <sup>δ</sup>15N% values and total crustacean biomass of the two size fractions analyzed (≥<sup>450</sup> <sup>μ</sup>m and <sup>≥</sup><sup>850</sup> <sup>μ</sup>m). A remarkably higher isotopic carbon (or less depleted) signature in summer and a higher nitrogen signature in autumn and winter was observed in the samples of both size fractions. The most depleted carbon signatures in the fraction ≥850 μm occurred with high relative abundance of primary consumers

to total zooplankton biomass, concomitant to the spring peak or abundance of the cladoceran filter feeder *Daphnia*. In the ≥450 μm fraction, where most depleted carbon signatures were recorded, cyclopoids and diaptomids were also important contributors to the total biomass, together with Cladocera filter feeders.

**Figure 2.** (**A**) Biomass percentage composition of planktonic crustacean taxa (bubbles) of Lake Maggiore and their δ13C% values (bubble centers) and δ13C% values of pooled samples (black line); (**B**) biomass percentage composition of planktonic crustacean taxa (bubbles) of Lake Maggiore and their δ15N% values (bubble centers) and δ15N% values of pooled samples (black line); (**C**) total biomass of Lake Maggiore planktonic crustacean taxa. All data refer to the fraction ≥450 μm.

Zooplankton standing stock biomass largely changed with the seasons. In samples belonging to the ≥450 μm size fraction, a spring peak was followed by a second phase of biomass increase in autumn.

In the ≥850 μm size fraction, peak biomass values were generally recorded in August. The unusual increase in May 2015 was related to the detection of the peak population growth of *Daphnia*. As already highlighted, *Daphnia* can be an important contributor to standing stock biomass, due to its peak biomass in spring being related to the development of young, smaller specimens.

We observed a clear seasonality of <sup>δ</sup>13C isotopic signatures for the <sup>≥</sup>450 <sup>μ</sup>m size fraction during the studied years. The least depleted values in δ13C% were regularly measured in August, while more depleted values characterized the spring samples (May), when the contribution to the peak biomass of primary consumers, particularly *Daphnia*, was at its peak. Limited seasonal variation regarding the contribution to the total biomass was characteristic of the cyclopoid copepods, whose contribution to carbon and nitrogen isotopic signatures showed similar values during the year. In the >450 μm size fraction, cyclopoids and calanoids were related to depleted values of carbon isotope signatures and to the highest values of nitrogen signatures. The biomass contributions of different zooplankton taxa, both primary and secondary consumers, drove patterns of both isotopic signatures; therefore, the measured values, e.g., in January 2012, integrated isotopic signatures of the two taxa present, i.e., cyclopoid copepods and, to a lesser extent, *Daphnia*.

**Figure 3.** (**A**) Biomass percentage composition of planktonic crustacean taxa (bubbles) of Lake Maggiore and their δ13C% values (bubble centers) and δ13C% values of pooled samples (black line); (**B**) biomass percentage composition of planktonic crustacean taxa (bubbles) of Lake Maggiore and their δ15N% values (bubble centers) and δ15N% values of pooled samples (black line); (**C**) total biomass of Lake Maggiore planktonic crustacean taxa. All data refer to the fraction ≥850 μm.

A marked seasonality was also characteristic of the <sup>δ</sup>13C and <sup>δ</sup>15N signatures in the <sup>≥</sup>850 <sup>μ</sup>m size fraction. The least carbon-depleted values in August were related to the high contribution of secondary consumers, while more carbon-depleted values in spring and winter were related to a higher contribution of primary consumers. Values tended to be most enriched in winter, with secondary consumers more enriched than primary ones. The signatures of <sup>δ</sup>15N% measured in the <sup>≥</sup>850 <sup>μ</sup>m size fraction were related to the relative biomass contribution of primary and secondary consumers. Among predatory cladocerans, *Bythotrephes* was generally much more δ15N enriched than *Leptodora*. For example, in August 2013, the nitrogen enrichment of *Bythotrephes* was higher in comparison to *Leptodora* and *Daphnia*. Differences in δ13C signatures suggest that *Leptodora* was living in deeper waters than *Bythotrephes* and *Daphnia*. Changes in isotopic signature and contribution to total biomass were mainly related to cladoceran primary and secondary consumers, thus enabling the investigation of infra-Cladocera food relationships.

As hypothesized, we found that carbon and nitrogen isotopic signatures of the two pooled size fractions could be reconstructed from taxa-specific signatures weighted on contribution to total biomass, by applying the following equation:

$$\delta^{15} \text{N}\_{\text{\textquotedblleft}} \delta^{13} \text{C}^{\text{\textquotedblleft}} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblleft} = \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblright} \text{\textquotedblleft} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblright} + ... + \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblright} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedblleft} \text{\textquotedbl$$

where *i* = δ15N% and δ13C% are taxa-specific signatures.

Relationships between measured and estimated isotopic signatures were significant (δ13C%: R2 = 0.904, N = 28, F = 244.682, *P* < 0.001; δ15N%: R<sup>2</sup> = 0.942, N = 28, F = 438.413, *P* < 0.001; Figure 4). Such relationships are the result of the accuracy of estimation of both biomass and taxa-specific isotopic signatures. In fact, taxa-specific biomass estimation is based on length–weight regression equations applied to a consistent number of specimens (at least 25 specimens per taxa) and taxa-specific isotopic signatures analyses are based on a very high number of specimens (cfr. Section 3), which is necessary

for Isotope Ratio Mass Spectrometry (IRMS) analysis (at least 1 mg of dry weight). The accuracy of standing stock biomass determination from length–weight regression of taxa-specific analysis was also confirmed by the comparison with the direct biomass dry weight estimation, which was shown to be statistically significant (R<sup>2</sup> = 0.814, N = 8, F = 26.251 *P* = 0.002).

**Figure 4.** Relationship between isotopic fingerprint (carbon and nitrogen) of pooled zooplankton samples belonging to ≥450 μm and ≥850 μm size fractions determined by instrumental analysis (measured) and reconstructed values (estimated) following Equation (2).

The difference between the smaller and the larger size fraction was mainly determined by cladoceran taxa composition. In the case of DDTs, the main effects cannot be properly interpreted because a weak, but significant, interaction between size and seasonality was determined when two-way ANOVA was applied (N = 28, F = 3.149, *P* < 0.048), so that the size of a factor's effect depends upon the level of the other factor. Changes in PCBs and δ15N values were shown to be strongly driven by seasonality (PCB: N = 28, F = 14.312, *P* < 0.001; δ15N: N = 28, F = 29.169, *P* < 0.001) and for both datasets there was no statistically significant interaction between factor size and seasonality (*P* = 0.181 and *P* = 0.596, respectively). Results of the two-way ANOVA are reported in Table 1.

A regression analysis was performed with δ15N data and pollutant concentrations (Figure 5) with ≥450 μm and ≥850 μm size fractions. Changes in both DDT and PCB concentrations on a logarithmic scale were related to zooplankton δ15N% isotopic signatures (DDTs: R2 = 0.541, N = 26, F = 28,340, *P* < 0.001; PCBs: R<sup>2</sup> = 0.502, N = 27, F = 25.159, *P* < 0.001). The data in the upper right side of the graph

refer to autumn and winter samples, while those in the lower side of the graph are related to spring and summer samples.


**Table 1.** Results of the two-way ANOVA tests performed on DDTsPCBs concentrations and δ15N% fingerprint of zooplankton samples of the two size fractions (≥450 μm and ≥850 μm).

**Figure 5.** Relationship between nitrogen isotopic fingerprint of pooled zooplankton samples from Lake Maggiore, belonging to ≥450 μm and ≥850 μm size fractions, and persistent organic pollutant (POPs) concentrations, expressed as logarithmic scale.

POP concentrations in the ≥850 μm size fraction were generally lower than in the ≥450 μm size fraction (Figure A1).

#### **5. Discussion**

The pelagic zooplankton is composed of a large variety of organisms differing in taxonomic composition and body size. Together, protists, monogonont rotifers, Cladocera, and Copepoda contribute to the zooplankton community, the latter with both individuals at their adult stage and a number of developmental stages that are different in size and trophic role. Regarding net zooplankton, as a result of the mesh size, only rotifers and crustaceans were reliably sampled. The former can be numerically dominant during the spring phase, along with naupliar and early copepodite stages of copepods. Their contribution to zooplankton biomass is little, however, with respect to Cladocera and adult or sub-adult copepods [38–40]. Altogether, the latter are ≥450 μm in size.

Zooplanktivorous fish, however, tend to select larger (≥850 μm) and more visible prey, such as *Bythotrephes*, large, ovigerous *Daphnia*, and *Leptodora*. Therefore, these are a component of the zooplankton population, which is active in transferring POPs to fish [10,41–44].

By choosing to focus on the two size fractions of crustacean zooplankton, we aimed to investigate components which are directly involved in the transfer of pollutants in the trophic chain. In addition, investigating the role of micro-zooplankton is almost impossible, given that POP concentrations and δ13C and δ15N analyses require high weight amounts of zooplankton material (approximately 1 μg dry weight and 1 mg dry weight, respectively). Further studies are required to elucidate the role of micro-zooplankton as carrier of pollutants.

As expected, zooplankton biomass composition largely varied with the season in both size fractions analyzed, because in deep subalpine lakes, the abundance and composition of zooplankton populations vary along the season [45–47]. A spring peak abundance of phytoplankton population triggers the growth of primary zooplankton consumers, such as the large filter-feeder *Daphnia*, which, in Lake Maggiore, usually reaches its peak population density in May [47]. The increase of *Daphnia* and other primary consumer cladocerans (*Eubosmina*, *Diaphanosoma*) in spring is followed by an increase in predatory cladocerans (*Bythotrephes* and *Leptodora*), able to selectively exploit primary consumers and contribute to their decline [48,49]. In turn, the decrease in zooplankton primary consumers promotes a second phytoplankton peak in summer, which leads to a second phase of increase in primary consumers in autumn.

In this study, seasonality was also observed in isotopic carbon and nitrogen zooplankton signatures. In pelagic zooplankton, δ13C% tended toward less depleted values in summer, likely mirroring phytoplankton isotopic signatures [38,39,50]. During thermal stratification, the high growth rates of phytoplankton cells can lead to the consumption of dissolved atmospheric CO2, causing a shift in carbon exploitation sources toward the uptake of the HCO3 − ion [51–54]. The unselective filter feeder *Daphnia* feeds on seston particles, which fit the distance between their intersetular filtering combs (i.e., 0.2–50 μm) [55]. Since Lake Maggiore seston is mainly composed of phytoplankton cells and carbon fractionation is negligible, *Daphnia* isotopic signature specimens matched variations of the phytoplankton carbon fingerprint [45,50]. The pelagic zooplankton nitrogen isotopic signature also varied with the season, with the maximum enrichment in winter. Isotopic nitrogen tended to increase along the trophic food chain and stepwise enrichment varied in different environments and with the structure of the trophic web [45,56,57]. Isotopic nitrogen enrichment is diet-dependent [58]. In eutrophic lakes, where high values of δ15N% are usually found, zooplankton organisms do not feed only on phytoplankton, but also exploit alternative organic particulate matter, such as bacteria, protozoa, exuviae, and fecal pellets of zooplankton specimens, which are more δ15N-enriched than phytoplankton algae. We can hypothesize that, in winter, when the abundance of phytoplankton cells is low, alternative, δ15N-enriched food sources (e.g., bacteria, protozoa, exuviae) became more important, leading to the observed enrichment in isotopic nitrogen content. In addition, under food shortage conditions, the zooplankton might be able to exploit δ15N%-enriched lipids which accumulate in high food concentration conditions [59,60].

Distinguishing between primary and secondary consumers is crucial for the reconstruction of patterns of matter and energy through the pelagic food web. As primary and secondary consumers largely overlap each other in body size, such distinction cannot be based on size fraction analyses of zooplankton. We found that the isotopic signature of pooled zooplankton samples of the two size analyzed fractions was tightly related to taxa-specific contribution to total biomass, thus allowing for the identification of taxa-specific contribution to measured δ15N% isotopic signatures. The latter is, in turn, correlated to POP concentrations, therefore allowing for a better definition of primary and secondary zooplankton consumer taxa in the transfer of toxicants through the pelagic food web. The significant relationship between pooled zooplankton samples and taxa-specific contributions to the total biomass in Lake Maggiore results from the fact that other zooplankton particles are excluded when samples are filtered through 450 μm nets. In Lake Maggiore, detritus is usually <126 μm and the contribution of phytoplankton colonies >76 μm is largely negligible. Therefore, any other zooplankton contribution to the total biomass in the samples of the two size fractions analyzed was certainly negligible. This is not the case in eutrophic or mesotrophic lakes, where zooplankton is smaller and detritus and phytoplankton colonies can be found in zooplankton samples in non-negligible amounts [1,61].

We hypothesized that both size and seasonality can influence POP concentrations in zooplankton populations and that the two parameters are related to the δ15N% isotopic signature. Our results highlighted a predominance of seasonality versus size. The ≥450 μm fraction was composed of primary (*Daphnia*, *Eubosmina*, *Diaphanosoma* and diaptomids) and secondary consumers (cyclopoids, *Bythotrephes* and *Leptodora*) throughout the year, while in the ≥850 μm fraction, the biomass of secondary consumers is dominant in August. As secondary consumers are expected to be higher in POP concentrations, we would have expected to record the highest values of POP concentrations in summer; on the contrary, we observed high POP concentrations in winter, when zooplankton nitrogen enrichment was at its maximum. In fact, the δ15N% and POP results directly correlated, with both tending to increase from spring to winter. Ecotoxicological studies on marine trophic chains demonstrated a similar relationship along the food chain [62,63]. In Lake Maggiore, the seasonal pattern of POP concentrations was very similar terms of values between years of the study period [28], suggesting a steady-state condition. In these conditions, the major δ15N-enrichment of carnivorous species was probably hidden by seasonal variations in δ15N%, which was dependent on quality of organic particles, i.e., the bacteria, protozoa, or organic particles deriving from dead organisms should have higher POP concentrations than phytoplankton cells.

The predominance of seasonality versus size was confirmed also by DDT concentrations in the smallest size fraction. Theoretically, the big predators, *Bythotrephes* and *Leptodora*, and the largest *Daphnia* specimens of the ≥850 μm size fraction should be characterized by higher values of POP concentrations, the former due to this specimen belonging to an upper trophic level and the latter because adults of a larger size can accumulate more pollutants in their tissues. However, our results likely reflect the fact that copepods rich in lipids (substance where POPs are mainly stored) were not present in the larger size fraction. However, because zooplanktivorous fish mainly select large cladocerans [43], the ≥850 μm fraction should be more representative of fish food sources.

#### **6. Conclusions**

Zooplankton play a crucial activity in matter and energy transfer in the food web. In temperate lakes, biomass, taxa, and size composition largely vary with the season. The contribution of primary and secondary zooplankton consumers also varies seasonally, along with the δ15N% isotopic signature. By investigating seasonal changes in biomass and taxa-specific contributions and in carbon and nitrogen isotopic signatures along with DDT and PCB concentrations, we aimed at quantifying the relative importance of size versus seasonality. We found that the δ15N% of two size fractions, representative of bulk crustacean zooplankton biomass and of the fraction on which fish mainly prey, respectively, significantly correlated to POP concentration, as both varied with the season. The bulk zooplankton nitrogen isotopic signature resulted from a biomass-weighted taxa-specific isotopic signature, thus enabling us to distinguish between the contributions of secondary and primary

consumers to measured nitrogen enrichment, from which the concentration of pollutants depends. DDT concentration was higher in the >450 μm fraction, while no difference was found in PCB concentration in the two size fractions. The difference between the two was mainly related to copepod adults, which were entirely lacking in the larger, >850 μm size fraction, the one on which fish selectively prey. Overall, seasonality was largely predominant over size for the dynamics of DDTs and PCBs, therefore suggesting that, in temperate lakes, more than one time spot studies are required; indeed, multi-year studies are ideal for verifying, among others, the condition of a steady state which must be fulfilled before applying models regarding the transfer of pollutants through an ecosystem.

The oligotrophic, subalpine Lake Maggiore was ideal for investigating the role of zooplankton in pollutant transfer, as size spectra distinction allows for the neglect of background noise caused by micro-zooplankton, protists, and phytoplankton algae. This is not the case for more productive lakes, in which such size-based distinction cannot be applied in the same way.

**Author Contributions:** Formal analysis, M.M.M. and R.P.; SIA and zooplankton investigation, M.M.M. and R.P.; POPs investigations, M.M. and R.B.; writing—original draft preparation, M.M.M., R.P., R.C., and D.C.; Review, R.B. and M.M.

**Funding:** This research was funded by CIPAIS (Commissione Internazionale per la Protezione delle Acque Italo-Svizzere).

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

#### **Appendix A**

**Figure A1.** Concentration of POPs in two crustacean zooplankton size fractions (≥450 μm and ≥850 μm) of Lake Maggiore.

#### **References**


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