The non-indigenous Oithona davisae in a Mediterranean transitional environment: coexistence patterns with competing species

The Venice lagoon (VL) has been recognized as a hot spot of introduction of non-indigenous species (NIS), due to several anthropogenic factors and environmental stressors that combined may facilitate NIS invasions. In the last decades an increasing number of zooplankton NIS have been observed in the VL. This work aims to provide a picture of the annual cycle and distribution of the recently recorded non-indigenous copepod Oithona davisae, considering the coexistence patterns with the congeneric resident Oithona nana. Therefore, zooplankton samplings were carried out monthly from August 2016 to July 2017 at five Long-Term Ecological Research LTER stations in the VL. Oithona davisae showed a persistent occurrence throughout the year with the highest abundances in the warm season and in the inner areas, while the congeneric O. nana, showing a different distribution pattern, resulted more abundant near the inlets of the Lagoon, where O. davisae reached the minimum density. Oithona davisae seems to find local conditions that promote its settlement and distribution, especially in the inner and more trophic lagoon sites. In other European coastal embayments or transitional waters, O. davisae occupied the niche left by the indigenous O. nana or can replace this congeneric species through competitive exclusion mechanisms. Our data indicate that, for now, such species replacement has not occurred in the VL. One of the causes is the extreme variety of habitats and niches offered by this environment allowing a balanced coexistence with O. nana and in general with the resident copepod community.

Research Network (LTER-Italy, www. lteri talia. it). It is a shallow coastal water body located in the north-western part of the Adriatic Sea, with a surface area of about 500 km 2 ( Fig. 1) and an average depth of 1.5 m. It is connected with the Adriatic Sea through three inlets allowing water and sediment exchange with the sea, named Lido, Malamocco and Chioggia. The semi-diurnal microtidal regime of the VL has with a mean range of 0.40 m during neap tides and about 0.80 m during spring tides. The residence times ranged from few days, in the proximity of the three inlets, to over 60 days in the inner lagoon areas 28 , while the amount of seawater that is exchanged during each tidal cycle is about one third of the total volume of the lagoon 29 . The brackish nature of the VL results from the freshwater inputs from the drainage basin, which is about 35 m 3 s −130 . The lagoon is a heterogeneous and complex system, characterized by different physical and environmental gradients, and a mosaic of landforms and habitats that are the result of natural and man-induced drivers 31 ; in the entire lagoon basin, typical salinity and temperature range of 9-37 and 5-28 °C, respectively 32,33 . The VL is affected by several anthropogenic pressures like fisheries, shellfish culture, domestic and industrial discharges and maritime transport 34,35 . From a biological point of view, planktonic copepods represent the dominant fraction of zooplankton assemblages across the VL, albeit with varying species composition along the salinity gradient 37 . In the VL, resident zooplankton species are particularly able to adapt to highly variable environmental conditions. In the inner areas, characterized by salt marshes and channels, the copepod A. tonsa, a typical species in transitional waters, is dominant in the warm period 36 . Towards the inlets, with the increasing of marine features, neritic and marine species such as Acartia clausi, Paracalanus parvus and Centropages spp. begin to dominate 38 . In some periods, especially during the warm season and in areas close to the inlets, the presence of cladocerans is relevant (Evadne spp., Podon spp. and Penilia avirostris), while in the internal lagoon basin, meroplanktonic forms (such as larval stages of Decapoda, Bivalvia and Polychaeta) represent an important fraction of the mesozooplankton community 37,38 . Within oithonids, only O. nana represents a dominant component of the summer zooplankton lagoon assemblage, while O. similis and O. plumifera, the other two oithonids detected in the Lagoon, were generally found in areas with more marine characteristics and always with low abundances 39 . Zooplankton sampling was carried out in the context of the LTER monitoring activities, which studied since few decades five stations located in the central and the northern part of the Venice Lagoon ( Fig. 1) and representing a gradient of environmental conditions and different anthropogenic impacts. Station St.1 (S. Giuliano) is affected by urban discharges, St.2 (Marghera) is located in the industrial port area influenced by industrial pollutants and St.3 (Fusina) receives warm water from a thermoelectric power plant. These three stations, considered the inners stations, are characterized by high water residence times of 18-25 days 28 . The inlet station St.4 (Lido), located in the northernmost inlet of the Lagoon, is influenced by Adriatic coastal waters with a residence time of few hours and finally, St.5 (Palude della Rosa) is a typical lagoon environment influenced both by fresh water and by sea water inputs and characterized by a residence time of 6-12 days 28 .
Sampling methods and laboratory analysis. To 31 . This particularly small mesh size was adopted, as it is known that both O. davisae and O. nana, being relatively small species, may pass through the classically used nets with a 200 μm mesh 40 . The characterization of the rest of the mesozooplankton community was based instead on seasonal samplings, from May 2014 to February 2017 (Supplementary  Table S4), performed by the ongoing long-term zooplankton monitoring activities. Those zooplankton community samples were collected by horizontal hauls of HydroBios Apstein net (mouth diameter 0.4 m, mesh size 200 μm) and of WP2 standard net (mouth diameter 0,57 m, mesh size 200 µm) equipped with a flowmeter HydroBios to calculate the filtered water volumes (7 to 50 m 3 ).
A multiparametric CTD probe (SBE 19plus) was used to measure environmental parameters like temperature, salinity, oxygen and chlorophyll a at the sampling sites, while the water transparency was measured with the Secchi disk. The zooplankton samples were fixed in a solution of 4% sodium tetraborate buffered formalin and ambient water. In the laboratory, copepods were identified and counted, according to the International Council for the Exploration of the Sea (ICES) protocol 41 to species level, where possible, or to higher taxonomic levels. As males of the genus Oithona could not be easily identified at species level, only copepodites and females were considered regarding the two species O. davisae and O. nana. While adult males and earlier stages (nauplii) were not counted.
Statistical analysis/data processing. Spatial and temporal patterns of the environmental factors based on Euclidean distances were assessed using repeated-measure permutational analysis of variance (PER-MANOVA) with the sampling months as fixed factor and the stations as random factor (PRIMER 6 + and PER-MANOVA software package; PRIMER-E, Ltd., UK). In order to visualize the similarities between the samples in terms of environmental conditions, a PCoA (Principal coordinates analysis) based on Euclidean distances and the vectors of environmental parameters was performed, and the occurrence of the two oithonids, O. nana and O. davisae, (as Bray-Curtis-Similarity matrix) was superimposed. Moreover, the Spearman correlation of the environmental and community similarity matrices were calculated with the software package PRIMER6+, utilising the package RELATE 42 in order to estimate the relationship between the environmental factors and the abundances of the two Oithona species.
The Pearson correlations to estimate the relationship between the abundances of the two Oithona species and the environmental parameters were calculated with the Software package Statgraphics 43  www.nature.com/scientificreports/ between and within groups of abundance data (months, stations), ANOVA LSD test was employed after the checking for the homogeneity of the variances.

Rank abundance distributions-RAD. The characterisation of abundance and distribution of most
Oithona species calls for the use of sampling nets with a mesh size smaller than 200 μm 12,40 . Parallel comparison of efficiency of different sampling devices in catching Oithona species has been only seldom faced and investigated 15 and no conclusive guidelines have been issued to scale the catches or merge the data collected with different approaches in single tows. Nevertheless, for the scope of the present study, it was necessary to detail the presence/absence and relative abundances of to the two species (O. davisae and O. nana) within the picture of the overall copepod community in each station and possibly in each seasonal/environmental condition.
Given the quite even distribution of the samplings events during the year (monthly with the 80 μm net, and almost seasonally during 2014-2017 period with the 200 μm net), we approached the issue by assuming as "typical" the yearly trend of Oithona species coming from the samplings with 80 μm net ( Supplementary Fig. S1). The linear interpolation of the expected abundances in specific calendar days made possible their merging with the specimen counts of the 200 μm samplings.
In this way we were able to construct the synthetic Rank Abundance Distributions (RAD) of the overall copepod community in each station by adding together all the catches collected with both nets (80 μm and 200 μm) and then calculating the relative abundances and ranks of every species. In particular we obtained, for each station, a pool of copepod communities that (1) integrates the species abundances of standard 200 μm net and the 80 μm net counts for O. nana and O. davisae, and (2) corresponds to several seasonal conditions during the year. These pools were also used to investigate the mechanisms of niche apportionment.
Niche apportionment models. To give insights into the mechanisms driving the structuration of the community, captured by the relative abundance distributions of species, the Niche Apportionment Models (NAM) were used.
The theoretical framework regarding the NAMs has been developed in early 90 s by Tokeshi 44 who built on Sugihara's niche-hierarchy models (1980), to support in the practical problems of fitting stochastic models to real data for the interpretation of species abundance patterns. The developed framework gives logical coherence to a range of niche apportionment models that can be categorized through the sequential breakage process of total niche. To describe how species break up the resource pool, determining the distribution of abundance of different species 44,45 , the theory proposed by Tokeshi [45][46][47] takes into consideration six different models which are (sorted by decreasing evenness): Dominance Preemption (DP), Random Assortment (RA), Random Fraction (RF), Power Fraction (PF), MacArthur Fraction (MF) and Dominance Decay (DD); the abundance of a species is assumed to directly correspond to the amount of niche/resource apportioned to that species.
In these models of sequential breakage of the total niche, the dominance pattern depends on the probability p associated to each niche fraction. The resulting Species Abundance Distribution (SAD) vary from uneven patterns like the geometric series to extremely even ones similar to the broken-stick model 48 .
The power of this analysis ultimately relies on the robustness of the Species Abundance Distribution within a community in catching (in the shape of the SAD itself) the underlying mechanisms that drive the relations (coexistence, competition, etc.) among the species and therefore the structuration of the community. So that fitting one of the proposed stochastic models to the SAD (i.e. to find the model that fits better than others to the survey data, and this can be done with classical statistical methods) means that the probability of that specific sequential niche breakage (caught by the modelled scheme) is higher than others.
In order to highlight the mechanisms regulating the resource allocation within a community of related species, several authors suggest applying the NAMs only on that specific group of taxonomically related species or trophic guild [45][46][47] . As copepods represent about 80% of the total zooplankton community in the VL 38 , we applied the NAMs to this dominant planktonic group in order to elucidate the mechanisms regulating the resource allocation in different areas and to understand if those mechanisms could explain the different distribution of O. davisae and O. nana. This is the first attempt to apply this approach to copepods communities. In fact, to date, these models have been used to explain and describe changes in relative abundance distributions of terrestrial and aquatic communities of both fresh and marine waters, like larval chironomid, stream fish assemblages, salt meadow vegetation, dragonfly community, terrestrial arthropods, marine macroalgae and phytoplankton community 44,46,[49][50][51] . The analyses have been performed using the R-package 'nicheApport' (available in GitHub, 2017).

Results
Environmental characteristics. The investigated area is characterized by a high variability of abiotic conditions typical of transitional environments. The environmental parameters differ significantly both spatially, between stations, as well as temporally, between the month of sampling (PERMANOVA Table 1). During summer, higher temperatures and Chl-a values and lower transparency were measured compared to other seasons. The salinity, influenced by fresh and marine water inputs, varies between stations and throughout the year. Trophic and salinity gradients are evident with the highest mean Chl-a values at the inner stations and the highest mean salinity value in the area close to the inlet (St. 4) (     (Fig. 4).
O. nana generally presented abundances higher than O. davisae in about 87% of the samples. Moreover, in 65% of the samples its relative abundance, within these two oithonids, was greater than 90%; while O. davisae reached relative abundances of more than 90% only at St.1 in February, from April to July and in September (Fig. 4).
After one year of observation, it emerged that the two species seem to be distributed along opposite physicochemical gradients. While O. nana resulted to be more abundant near the inlets, where O. davisae reaches its minimum, O. davisae is more abundant in the inner areas, where O. nana occurs with low abundances (Fig. 4). When calculating the Pearson correlation of the two species to every single parameter, it is evident, that the presence of O. davisae is significantly positively correlated to Chl-a (r = 0.358, p = 0.005), while O. nana is weakly though significantly positively correlated to salinity (r = 0.295, p = 0.0231) and significantly to transparency (r = 0.43, p = 0.0007) (Table 1, Fig. 2). In general, our data show weak correlations with all the environmental parameters considered, except for chlorophyll a in which the correlation is medium and significant; with respect to temperature, the correlation is at the limits of significance (for O. davisae) and there is a positive correlation signal which supports the seasonal characteristics of the species. The correlation with chlorophyll is however a clear signal of a larger scale trophic dynamics, which supports a trophic network, also of prey and predators, within which O. davisae is placed.     (Fig. 5).  (Table 3). This reflects, for St.1, an increasing trend in the evenness of the community from spring to winter: in spring, the best fitting model was the Dominance Preemption, characterized by a steeper RAD curve (Fig. 6). This model predicts that the first species that invades a resource space always occupies more than half of the available niche letting the remaining portion to be exploited by a second species in the same way and so on. Therefore, this species will be dominant on all the remaining species, resulting in a strong dominance structure. In summer and autumn, the best fitting model for both stations was the Random Assortment model (Fig. 6), which refers to a situation where the niche space is randomly divided among simultaneously or sequentially colonizing species. Thus, the abundances of different species vary independently of each other or, in other words, the species appear to be ecologically identical according to the neutral theory of biodiversity. This could arise both as a result of a mismatch between abundance and niche size or as dynamic non-hierarchical niche partitioning in a highly variable environment 45 . During winter, the best fit was for the Random Fraction model, with a RAD curve with a flatter slope than the previous ones, indicating a higher degree of evenness within the copepod community in this season (Fig. 6). This model is based on the assumption that sequentially colonizing species in an assemblage have the same probability of niche occupancy so that each new species carves out a random portion of the available niche space 44,46,52 . In this case, the coexistence mechanism during winter, characterized by a more diverse community, could rely on the development of independent population dynamics with random exponential growths under a variety of environmental limiting factors, as in Preston's lognormal distributions. At St.4, the best fitting was for Random Assortment model all year round, but also a weak signal for Random Fraction was detected in spring, autumn and winter ( Table 3). None of the schemes MacArthur Fraction and Dominance Decay exhibited a detectable fit.

Discussion
The increasing anthropogenic activities, including maritime transport and aquaculture, enhance the spread of NIS at a global level 53,54 , potentially increasing the local biodiversity or endangering local communities. Climate change and the resulting increase in seawater temperature, especially in Mediterranean, as a marginal sea, and particularly in Northern Adriatic 55,56 , favour the spread of thermophilic invasive species. Alien species that occupy trophic niches, similar or close to those of native species, can alter the balance of the latter and lead to a significant reduction in their number 57,58 . The invasive copepod O. davisae, known to have potentially a negative impact on the local copepod community, is often reported among the dominant species within the copepod community, both in native and introduced areas 17,23,59 . In the past, the genus Oithona has been probably often underestimated since the coastal species are less than 1 mm in body length and the use of conventional mesozooplankton nets (200 μm mesh size) fails to quantitatively catch these small cyclopoid species, resulting therefore unsuitable to characterize their real abundance 12,40 . With respect to previous works [e.g. 9 ], the present study adopts a net equipped with a smaller mesh size, allowing to detect the presence of the species with greater precision both in terms of quantity and distribution area. The present study reports the first description of the spatial and temporal distribution of the non-indigenous cyclopoid copepod O. davisae in the VL, after its first observation in 2014 9 , including details regarding potential invasion patterns and the niche relationships with the autochthonous congener O. nana and the overall copepod community. In this study, O. davisae revealed patterns of presence and abundances related to the seasonality, with highest densities in the warm period and very low abundances in the cold season. These data confirm the thermophilic character that this nearshore estuarine species exhibits not only in native areas but also where it is a recentlyintroduced player 20,60 . In the VL, the O. davisae population began to increase in May, when water temperature exceeded the values of 18-20 °C, threshold at which the population densities start to actively increase 61 . The maximum abundances in the VL were reached in June and July, the same months in which 18 reported the highest densities of O. davisae in Lake Faro and Ganzirri (Strait of Messina, Italy), although with 20-fold greater than those observed in the VL. From August onward, the population exhibited a decrease in abundance. This trend  www.nature.com/scientificreports/ could be also related, as reported for other areas 17,61 , to direct jellyfish predation. In our case by the invasive comb jelly Mnemiopsis leidyi, which started to become very abundant in the VL from August 2016 (personal observation). We do not know if Mnemiopsis, which is a generalist predator, is the main predator of O. davisae in VL, but we hypothesize that the decrease of O. davisae in August and September could also be due to the massive presence of Mnemiopsis. According to available literature, in the native areas and in those where it has been introduced, O. davisae presents a second peak of abundance in late summer/early autumn 20,59-61 . Conversely, this has not been observed in the VL where they actually exhibited very low densities in the same period, also depending on the sudden drop of water temperature and low chlorophyll a values. In winter, when the water temperature further decreased down to 1.2 °C, the species resulted almost absent, or it was not detected due to very low abundances. In the cold period, it is known that pre-fertilized O. davisae females reduce motion and respiratory activity by holding the sperms in the spermatheca until spring 61,62 . Therefore, these survival strategies may explain how the species copes with these extreme conditions despite its thermophilic characteristics. However, as mentioned above, in the VL the maximum abundances of O. davisae are not comparable with those reached in the native areas 17 , or in the Black Sea 63 and in the Lakes of Faro and Ganzirri 18 . Another hypothesis regarding such relevant differences in abundance, could be related to the qualitative difference in food resources between VL and the aforementioned areas. The dominance of O. davisae, both in the Tokyo Bay and in Black Sea, was related to a shift from a diatom-based phytoplankton community towards one with prevalence of small flagellates 63,64 . Oithona davisae is, in fact, an ambush feeder that prefers motile preys like flagellates, ciliates and dinoflagellates 13,65,66 . Several authors reported these feeding preferences, both in laboratory experiments and in natural conditions, confirming that diatoms are not part of the diet of this species [67][68][69] . The second hypothesis explaining the low density of O. davisae in VL could therefore take into account the availability of suitable food. As described by 70 , the phytoplankton seasonal cycle in the VL shows a temporal distribution pattern which www.nature.com/scientificreports/ justifies the presence of O. davisae especially in certain areas and at certain times of the year: diatoms are in fact the dominant group all year round but, after their spring bloom, going toward summer, also cryptophyceans, dinoflagellates and euglenophyceans, appreciated food by O. davisae, appear in the phytoplankton community. This is consistent with our data and explains the high abundance of O. davisae only in June/July. This feeding preference may also explain the spatial distribution of O. davisae as it shows its highest densities especially at St.1 one of the most eutrophic areas in the lagoon 71,72 , where cryptophyceans, dinoflagellates and euglenophyceans were found to be abundant during the warm season 70,73 . O. davisae, as most copepods, is able to form dense populations in short time 17 and it can have a strong impact on the ecosystem structure and functioning especially following an ecological crash. In fact, in 1992, the copepod A. tonsa found the ideal conditions for its development and settlement during a critical event of hypertrophy followed by intense anoxia in the VL. In the case of the Black Sea, the successful introduction and later establishment of O. davisae has been preconditioned by the complete disappearance of O. nana, previously caused by the invasion of M. leidyi, enabling O. davisae to occupy this empty niche 74 , and later, by the reduced predation pressure of M. leidyi on zooplankton due to the predation of Beroe ovata on the ctenophore. Similarly, in Lakes Faro and Ganzirri (Sicily), the occasional remixing of anoxic bottom waters rich in hydrogen sulfide, in conjunction with high temperatures during summer (> 32 °C 75 ), caused the complete disappearance of zooplankton (Zagami, pers. comm). The temporarily underexploited niche, in the Lakes Faro and Ganzirri, could have then provided an opportunity for the successful introduction and establishment of O. davisae. In the VL, at present, there have been no documented ecological crashes and the possible role of grazing pressure on zooplankton and on O. davisae in particular is still undocumented, though under study. However, rather than a niche left possibly empty, predation is expected to be the factor controlling the spread and dominance of zooplankton species in the Lagoon.
In addition to niche employment, replacement through a competitive exclusion mechanism was also observed 18,26 . This latter phenomenon has been explained assuming that the ecological traits of O. davisae may confer competitive advantage over O. nana or other species, like A. tonsa in the Black Sea 59 . At the moment, such species replacement has not yet occurred in the VL probably because of the extreme variety of habitats and environmental niches offered by this environment 76 allowing a balanced coexistence with O. nana and in general with the resident copepod community. In fact, according to our data, the population of O. davisae has not exceeded, in terms of abundance, those of O. nana except in the internal and more eutrophic station (St.1) in the hot period.
Based on the different distribution gradients and the two target stations, one characterized by higher degrees of marine influence (St.4) and the other by greater degrees of confinement (St.1), in the present work the comparison of the synthetic RAD curves allows to define an updated and more reliable picture of the species rank of O. davisae, including its role with respect to the autochthonous congeneric O. nana and other copepod species like A. tonsa, which is dominant in the same areas where O. davisae has also currently settled 36 .
At St.1, the inner station, the copepod community showed a lower species richness and evenness compared to St.4, the more marine station close to the inlet, confirming the difference between these two areas in term of habitat features. The more restrictive environmental conditions of the inner station (St.1) trigger in fact a selection among competing species, regulated by environmental factors, limiting the species richness and favouring the most rapidly growing and tolerant ones. The inlet station (St.4) conversely, due to the marine influence and the resulting mixing, with more stable environmental conditions, maintains a more even species richness, regulated by interspecific competition. At St.1, O. davisae and O. nana were, together with A. tonsa, the three most abundant species and ranked at the second and third places respectively, while at St.4, O. nana was at the first and O. davisae at the eighteenth rank place. Also, with respect to the calanoid copepod A. tonsa, which shares the same distribution pattern of O. davisae, these species seem to coexist or perhaps O. davisae is still controlled in terms of abundances by the presence of A. tonsa. Both species are opportunistic and benefit from more unstable environmental conditions. An analysis of the interannual abundance of copepods in the Sevastopol bay showed a reduction in the abundances of the previous invader A. tonsa after the arrival of O. davisae 59 . Nonetheless, in the VL, the establishment of O. davisae does not seem to have caused a significant reduction in the abundances of A. tonsa, probably because to low abundances of O. davisae until now. From its first report in the VL to date, O. davisae has not dominated over A. tonsa, but only the continued monitoring of the evolution of the abundance and distribution dynamics of these species will allow to define their respective roles and/or dominance within the trophic lagoon network over time.
To support inferences on whether the mechanisms that regulate the resource allocation among the copepod guilds could favor or not the diffusion of O. davisae in the Lagoon of Venice, we tested the performances of several model of niche apportionment (NAMs) on the seasonal copepod abundances in the two already mentioned areas of the Lagoon (St. 1 and St.4). Niche apportionment models support the investigation of the mechanisms of apportionment of the resource pools in different habitats. When a niche contender, generally a non-indigenous species (NIS), is introduced into a new environment, the scenarios that this will trigger can be various. At St.1, the picture highlighted by the NAMs well describes the changes in the resource partitioning mechanisms. Starting from spring, the low-richness community is strongly dominated by A. tonsa and O. davisae, as major competitors in pre-emptying the resource pool and subtracting niche space to the less abundant species. Then the picture gradually evolves in Summer and Autumn, into a situation in which the available niches are unsaturated, as the high variability and the rapid changes in environmental conditions do not allow the establishment of strong competition mechanisms between species. The mechanism of apportionment randomly assigns the niche space to species irrespective to their abundances, which vary independently from each other. Finally, in Winter more stable environmental conditions and low productivity lead instead to the establishment of competition mechanisms without strong dominance of one or few species in a relatively richer and more even community. www.nature.com/scientificreports/ trophic and physical, but also by the established mechanisms that regulate the resources allocation. The NAMs of St.4 showed a relatively constant situation throughout the year, with the best fitting for Random Assortment model, but also with a weak signal for Random Fraction. In our opinion, this could be the effect of the intrinsic characteristics of the station at the inlet, where water and organisms of marine and estuarine origin periodically mix: these results in a blended community where both randomness and competition mechanisms determine the relative ranks. As a consequence, the trophic and physical local habitat features, disliked by O. davisae, may limit its distribution to the more confined and eutrophic area and this is reflected by the nice partitioning mechanisms acting on the copepod guilds at this station.
To date, O. davisae is well-established in the VL, but we do not know yet if the species has reached its growing plateau. Mechanisms of competitive exclusion with O. nana, as reported for other areas, have not been observed in the VL. In our opinion the most credited hypothesis is that there is not overlapping between the preferred habitat by the two congeneric species and the extreme habitat variability leads to coexistence. However, continuing a constant monitoring with the approach of combined net mesh size will allow to get a more realistic assessment of the mesozooplankton community, thus observing the future evolution of the populations of the two species. This updated picture of the copepod community, with a view on the arrival of a new species, also highlights and confirms that the long-term monitoring programmes of zooplankton are an essential tool for sustainable development and environmental management. www.nature.com/scientificreports/