Impact of Colonizer Copepods on Zooplankton Structure and Diversity in Contrasting Estuaries

The impact of the occurrence of the non-indigenous species (NIS) Acartia tonsa, Oithona davisae, and Pseudodiaptomus marinus and of the range-expanding copepods Acartia bifilosa and Calanipeda aquaedulcis on the structure and diversity of zooplankton communities was analyzed using 18 years (1998–2015) time series from the contrasting estuaries of Bilbao and Urdaibai (Basque coast, Bay of Biscay). Changes in the structure of communities were assessed by using multivariate analyses of taxa abundances and changes in diversity by using descriptors of alpha, beta, and gamma diversity. The most evident changes occurred at the upper reaches of the estuary of Bilbao, where an abundant and less diverse brackish community, dominated by the NIS, A. tonsa and O. davisae, succeeded a low abundance and more diverse community of neritic origin. The later establishment of C. aquaedulcis was linked to further changes in the structure of the community and a progressive increase in diversity. The seasonal pattern of diversity at the inner estuary and the beta diversity in the estuary were also significantly affected by the arrival of the NIS and C. aquaedulcis. In contrast, the original low diversity brackish community of the estuary of Urdaibai, clearly dominated by A. bifilosa, was far less affected by the arrival of the same copepod NIS, and A. bifilosa remained as the species best related to the changes in zooplankton structure and diversity in the brackish habitats of the estuary.


Introduction
Globalization and increasing trade have fostered the introduction and establishment of non-indigenous species (NIS) beyond their natural distributional range in all ecosystems, and as a consequence, there is growing concern about the potential impacts at several levels (Thomaz et al. 2015).
This phenomenon is one of the most important ecological disturbances (Mollot et al. 2017), and warnings about the serious threats posed by NIS that become invasive to biodiversity (Vitousek et al. 1997;Strayer 2012;Simberloff and Vitule 2014) have been given. Although more studies have been carried out in terrestrial habitats than in aquatic ones (Jeschke et al. 2012;Lowry et al. 2012), the impact of human activities, in general, and of invasions, in particular, are likely greater in aquatic ecosystems than in terrestrial ecosystems (Thomaz et al. 2015). This is more evident in estuaries and connected habitats where human impact has caused a decrease of > 90% of formerly important species, destroyed > 65% of the main habitats, degraded water quality, and increased species invasions (Lotze et al. 2006). The increase in commercial shipping has made estuaries one of the most vulnerable ecosystems to the arrival of new species (Frisch et al. 2006), due to the unavoidable presence of a huge variety of organisms, including copepods, in ship ballast waters (Bailey 2015), and because NIS may be more resistant to adverse environmental conditions than native Communicated by Wim J. Kimmerer ones (Katsanevakis et al. 2014). Although studies on benthic species are much more abundant than on pelagic ones (Thomsen et al. 2014), the pelagic communities are pivotal in the functioning of many aquatic systems. Indeed, zooplankton play a key role in the pelagic food webs, being the link between primary producers and higher trophic levels and also in biogeochemical cycles and being also good sentinels of environmental changes in the marine environment due to their physiology and short life cycles (Webber et al. 2005;Ji et al. 2010).
NIS can affect the recipient community in highly variable ways, from changes in the density, even local extinction of native species, to changes in the structure, diversity, and functions (Thomsen et al. 2014;Gallardo et al. 2015;Chan and Briski 2017), which can influence ecosystem stability over time because species respond differentially to temporal environmental variations (Mccann 2000;Schindler et al. 2015). In many systems, the establishment of multiple NIS has occurred, and scientists need to address the consequences of multiple invasions within a community (Olden and Poff 2003). Nevertheless, it has to be born in mind also that many of the NIS effects on the ecosystems they colonize are not easy to observe because they may take a long time to reveal, and therefore, the full impact on biodiversity may not be apparent in the short term (Gallardo et al. 2015). Conversely, several theories suggest that the diversity and community structure of an habitat could also influence the invasion success of a NIS in that habitat. One of the forefathers of these ideas was Charles Elton who developed the so called diversity-invasibility hypothesis (DIH) which predicts that it is more difficult that biologically more diverse communities are invaded by novel species (Elton 1958). He argued that more diverse ecosystems should be more resistant to invasion by NIS, that is to say, that they would show a higher ecological resistance.
Although global-scale diversity trends are easier to detect under the current situation of global biodiversity loss, localscale diversity trends are more complex (Richirt et al. 2019), because they may be strongly influenced by local ecological context (Elahi et al. 2015). One of the main tools to study the changes in long-term monitoring of large numbers of species and ecosystems is diversity indices (Magurran et al. 2010;Vackár et al. 2012). There are different diversity types depending on the focus: alpha diversity which reflects the variability in a small area of homogeneous habitat, gamma diversity which indicates variability at the regional level (or within a system with different habitats), and beta diversity which gives us the turnover of taxa from one habitat to another (Ricklefs 2010). It is relevant to note that the effect of invasions in the organization of diversity in the space can differ from local species assemblages (alpha diversity) to diversity that accumulates from compositional differences between local assemblages (beta diversity), depending on the balance of processes that cause species heterogenization or homogenization between sites (Socolar et al. 2016). These indices are easy to calculate, and they are not only relevant for making science but also to inform policy, because diversity is something our societies care about (Aslaksen et al. 2015).
Zooplankton monitoring programs are an essential tool to improve our understanding and management of effects related to NIS in estuaries. These programs have recently started to focus on the variations in community structure and functional diversity based on indicators (Chiba et al. 2018), such as the Aichi Biodiversity Target indicators at the global scale described in the Convention on Biological Diversity (2016).
Accordingly, we made use of the monitoring program of the zooplankton community carried out since 1998 in two contrasting estuaries of the Bay of Biscay, i.e., the estuaries of Bilbao and Urdaibai, to analyze variations in community structure and diversity induced by multiple new copepod occurrences. Over the study period, we detected the arrival of the copepod NIS Acartia tonsa, Oithona davisae, and Pseudodiaptomus marinus in both estuaries, as well as the previously undetected copepods Acartia bifilosa and Calanipeda aquaedulcis in the estuary of Bilbao (Albaina et al. 2009;Barroeta et al. 2020). The latter two species were already found in the first studies on zooplankton carried out in the estuary of Urdaibai and other estuaries of the Basque coast in the late 1970s and early 1980s of the past century (Villate et al. 2004), so we do not know if they were also in the estuary of Bilbao before it became highly polluted. The present work aims to determine the role of the occurrence and settlement of NIS and C. aquaedulcis and A. bifilosa on the changes in the structure and diversity of the zooplankton communities of the estuaries of Bilbao and Urdaibai from 1998 to 2015. Moreover, the contrasting community structure and diversity conditions in these estuaries before the colonization allow us to discuss about their influence in invasion success.

Study Area
The estuary of Bilbao is a small (~ 23 km) and shallow (between 0.5 and 32 m deep) estuary located in the inner Bay of Biscay (43º 23′ N, 03º 07′ W) (Fig. 1). It is a meso-macrotidal estuary where euhaline waters dominate (Villate et al. 2013); the inner part is strongly stratified with salinity below the halocline around 30, whereas the outer part is partially mixed (Intxausti et al. 2012). Tidal flushing is relatively low; therefore, water residence time is medium-high (~ 29 days), being lower in the channelized upper and middle reaches than in the outer Abra embayment, and much lower in above halocline layers (0.1-1.6 days) than in below halocline layers (0.3-11.6 days) for most of the estuary length (Uriarte et al. 2014) (Table 1). The main tributaries are Ibaizabal and Nerbioi rivers entering at the head of the estuary, although there are other small streams draining into the middle reaches. For most of the twentieth century, it was a highly polluted system due to the high industrialization of the area and raw sewage discharges, which gave rise to very high heavy metal concentrations and a permanent hypoxia/anoxia situation (Villate et al. 2013;Irabien et al. 2018). However, in 1979, a rehabilitation program was initiated with a comprehensive plan for the sanitation of the metropolitan  (26, 30, 33, 34, and 35) that were sampled area of Bilbao, and there was also an industrial decline in the area. Since then, heavy metal concentrations have decreased (Fernandez-Ortiz De Vallejuelo et al. 2010), and dissolved oxygen (Villate et al. 2013) and diversity (Borja et al. 2006) have increased. Nevertheless, it is a completely man-modified estuary, where land reclamation has been very extensive, the whole estuary except the outer Abra zone became channelized (Cearreta et al. 2004). Port facilities located in this Abra zone are an important marine transport and logistics center in the European Atlantic Arc.
The estuary of Urdaibai is a smaller (~ 12.5 km) and shallower (3 m depth on average) meso-macrotidal estuary located in the inner Bay of Biscay (43º 22′ N, 02º 43′ W) (Fig. 1). The inner part of the estuary is partially stratified, with salinities below the halocline around 26, but the outer part is well mixed. Freshwater inputs of the main river flowing into the estuary (Oka river) are low, and seawater dominates at high tide (Iriarte et al. 2015), water residence time being low . At the upper and intermediate reaches, there are salt marshes and muddy intertidal areas, whereas extensive intertidal flats, mostly sandy ones, are found at the lower reaches (Table 1). The estuary of Urdaibai is part of a biosphere reserve and is one of the most important wetlands of the Iberian Peninsula. However, until 2003, dredging was a usual activity to allow launching ships built in a small shipyard located in the middle reaches. Nevertheless, since then, ships are only partially built in this shipyard, no dredging has been carried out, and the estuary has gradually become shallower (Monge-Ganuzas et al. 2013). The main source of pollution is a small sewage treatment plant located in the inner zone, the estuary having a low to moderate anthropogenic impact (Cotano and Villate 2006). In the outer reaches of the estuary, there are also two small recreational harbours (Mundaka and Busturia) and the port of Bermeo, which is located a few kilometers away from the mouth of the estuary and is one of the main Basque fishery ports in the Bay of Biscay.

Data Source
The biological and environmental data (1998-2015 period) used in this study were obtained from an ongoing monitoring program of the zooplankton community of the estuaries of Bilbao and Urdaibai. Because both systems have small estuary basins and are subject to torrential hydrological regimes and semidiurnal tidal regimes that produce strong temporal variations in the salinity zonation, a Lagrangian type of sampling strategy (Modéran et al. 2010) was used, sampling in water masses of specific salinities instead of doing it at fixed points. Thus, samples were collected monthly at the sites of 30, 33, 34, and 35 salinities in the estuary of Bilbao (B30, B33, B34, and B35) and 26, 30, 33, and 35 salinities in the estuary of Urdaibai (U26, U30, U33, and U35) during high tide (Fig. 1). The salinity values correspond to those measured at the zooplankton sampling depth. At each salinity site, zooplankton samples were collected below the halocline by 2-3 min horizontal tows using a 200-µm mesh size net equipped with a mechanical flow meter (Hydro-Bios) and subsequently preserved in 4% buffered formaldehyde for later counting and identification to the lowest possible taxonomic level under an inverted stereomicroscope. Zooplankton samples were diluted to a specific volume (10-100 ml) in order to obtain a suitable density of individuals, and we took a number of aliquots for the identification and counting of individuals until at least 100 individuals of the most abundant taxon and 30 individuals of the second and third most abundant taxa were counted.

Data Pre-Treatment
Zooplankton missing values (around 6%) of the time series were filled with the mean of the preceding and following month values (Head and Pepin 2010;Fanjul et al. 2018). Only those species that showed an abundance > 0.01% during the period 1998-2015 in data pooled for the two estuaries were taken into consideration for the calculation of the diversity indices using the abundance sorting method adapted from Ibañez et al. (1993). Zooplankton data were grouped and analyzed at two levels: (i) the zooplankton group level that included 19 groups and (ii) the copepod species level that included 25 taxa (Table 2). For the multivariate ordination analyses, taxa densities were transformed to log 10 (x + 1) (hereinafter referred to as log (x + 1)) so that data could approach the required normal distribution. In order to test changes in diversity indices due to the new colonizing copepods, we distinguished three periods according to their colonization process. Thus, period 1 lasted from 1998 to 2002 (before the occurrence of A. tonsa and O. davisae in large numbers), period 2 from 2003 to 2009 (after the establishment of A. tonsa and O. davisae and before the occurrence of P. marinus), and period 3 from 2010 to 2015 (after the occurrence of P. marinus and a marked increase of C. aquaedulcis).

Data Analyses
Multivariate ordination methods were used to model the variability in the taxonomic structure of zooplankton communities at each estuary using the Canoco v 4.55 software. A principal component analysis (PCA) of the zooplankton taxa densities (log (x + 1)) of each estuary was conducted using month as co-variable to remove seasonal variability.
The year scores of the first two axes for each salinity site were depicted, in order to show year-to-year variation patterns in zooplankton community structure. The position of each taxon in the first two axes was depicted in order to visualize which taxa contributed most to the main interannual and spatial (salinity sites) changes in the structure of communities. For their representation in figures, zooplankton taxa and copepod species were classified into the following: NIS copepods, A. bifilosa, C. aquaedulcis and brackish copepods, meroplankton groups, neritic groups and copepods, freshwater copepods, and other categories (benthic or nektobenthic groups). The diversity of zooplankton groups and copepod species was analyzed by calculating alpha, beta, and gamma diversity indices from raw abundance data using Primer 6 software. The Shannon index (H´ = − Ʃp i ·lnp i ; being p i = relative abundance) was calculated at a monthly scale for each salinity site of each estuary to estimate alpha diversity, and pooled for all salinity sites of both estuaries to estimate gamma diversity.
To determine beta diversity, the Whittaker index was calculated at a monthly scale for each estuary (Magurran 2004). The monthly variations from 1998 to 2015 of all diversity indices are shown in contour plots made using the Surfer ® 10 software (Golden Software, LLC). In order to establish the relationships of interannual changes in diversity (alpha and gamma) with those in the density of NIS species, A. bifilosa, C. aquaedulcis, key neritic copepod species, and zooplankton groups (the five most abundant ones were selected in each case) over the study period, dispersion plots of the selected taxa deseasonalized density vs. deseasonalized diversity values were used. For this purpose, the seasonal effect, which is the major temporal pattern of variation in plankton time series (Ribera D´Alcalà et al. 2004;Benedetti et al. 2019), was removed by using standardized anomalies of both copepod species/zooplankton group densities (log (x + 1)) and Shannon index values. Standardized anomalies were calculated as the difference between each value and the mean value for each month for the study period divided by the standard deviation (Taboada et al. 2019). To visualize covariation patterns between diversity and density of key taxa, linear regression fits were calculated and plotted. In order to eliminate relationships likely to be spurious due to data scarcity, linear regression lines were only shown when the number of cases was higher than 5% of the total in the data series.

Community Structure
The first component (axis 1) of the PCA of zooplankton taxa densities from the estuary of Bilbao accounted for 41.4% of the total zooplankton variability. The plot of year scores along axis 1 (Fig. 2a) revealed similar trends of interannual variation at B30, B33, and B34 sites, but the magnitude of the change increased with decreasing salinity and was highest at the innermost site (B30). At B34, however, a slight return trend towards the initial community structure was observed in the year 2009. At the outermost salinity region (B35), however, no clear trend of community change was observed. The ordination of zooplankton taxa along axis 1 (Fig. 3) emphasized the major role of the NIS A. tonsa and O. davisae as responsible for the zooplankton community changes during the study period together with the NIS P. marinus, the brackish copepods C. aquaedulcis and A. bifilosa, and bivalve larvae. All those taxa increased in abundance over time, particularly in the lowest salinity site. Several neritic copepods, such as Temora stylifera, Oncaea media, Oithona plumifera and Acartia clausi, and holoplankton groups such as siphonophores and cladocerans, also contributed to this major (axis 1) pattern of zooplankton variation, showing an opposite trend to that of the former taxa. The second component (axis 2) of the PCA accounted for a much lower percentage of the total variability (15.7%) and showed a similar interannual pattern of zooplankton variation at all salinity sites (Fig. 2a), which was mainly explained by the common patterns of increase in time of some Fig. 2 Year average scores on the first two axes of the PCAS of zooplankton taxa densities (log (x + 1)) at each salinity site in (a) the estuary of Bilbao and (b) the estuary of Urdaibai meroplankton (gastropod and bivalve larvae) and holoplankton (appendicularians) groups, and neritic copepods such as the PCPC assemblage and Oithona similis. Only freshwater copepods showed a slightly opposite trend to those taxa (Fig. 3).
In the estuary of Urdaibai, the first component (axis 1) of the PCA analysis accounted for a much lower variability of the zooplankton community (23.7% of the total variability) than in the estuary of Bilbao (Fig. 2). The plot of year scores along axis 1 (Fig. 2b) did not show clear trends of zooplankton change throughout the study period at any salinity and revealed that the main interannual fluctuations were similar at all salinities. The strongest change in the community structure of the entire estuary was that occurred in 2012 as compared to the rest of years. These fluctuations reflected a common pattern of variation of most neritic taxa, with the neritic copepods Euterpina acutifrons, Oithona nana, and the PCPC assemblage showing the highest contributions. Slightly opposite trends were observed mainly for the NIS A. tonsa. The second component (axis 2) accounted for 17.3% of the total variability and revealed an overall progressive change of zooplankton community from 1998-2002 to 2010-2015, which was similar at all salinities (Fig. 2b). This change was mainly produced by the density variations of the brackish NIS copepod A. tonsa, some meroplankton groups (larvae of gastropods, polychaetes, and cirripedes), and the brackish native copepods A. bifilosa and P. grani, which were opposite to those of some neritic holoplankters, i.e., siphonophores, doliolids, T. stylifera, O. plumifera, O. media, and the PCPC assemblage (Fig. 3).

Alpha Diversity
There were between-estuary and within-estuary (spatial) differences in the temporal variations of alpha diversity both for copepod species and zooplankton groups (Fig. 4). In the estuary of Bilbao, the values and seasonal variations (smallest diversity in late winter-early spring and highest in autumn for copepod species, and smallest in autumn-winter and highest in late spring-early summer for zooplankton groups) were rather similar across the entire estuary in the period 1998-2002. However, values and seasonal variations of alpha diversity showed higher differences landwards from B34, while they remained similar in the outer estuary (B35) in the next two periods. At the lowest salinity site (B30), alpha diversity for both copepod species and zooplankton groups decreased from 1998-2002 to 2003-2009, with the strongest decrease having occurred from 2002 to [2004][2005]. The seasonal pattern also showed interannual changes, and late summer-early autumn became the time of the year with lowest diversity values in the 2003-2009 period. In the last period (2010)(2011)(2012)(2013)(2014)(2015), however, at B30, the seasonal changes in diversity differed between copepod species and zooplankton groups. The diversity of copepod species increased and the seasonal pattern was similar to that of the previous period, while the diversity of zooplankton groups decreased, and the seasonal pattern showed higher variation between years. At intermediate salinities, the interannual and seasonal changes of alpha diversity for both copepod species and zooplankton groups at B33 were more similar to those at B30 while those of B34 were more similar to those of B35.
In the estuary of Urdaibai, alpha diversity of copepod species and zooplankton groups decreased, in general, from high salinity sites (U35 and U33) to the lowest salinity site (U26) from the beginning of the study period, and showed large interannual fluctuations, mainly at the innermost site (U26), but no clear trends over the study period. The seasonal patterns of alpha diversity of copepod species and zooplankton groups showed variations between years and periods, but such variations decreased with increasing salinity in the estuary. At the outermost site (U35) the seasonal pattern of diversity was rather similar over the study period, generally with lowest diversity in early spring for copepod species and autumn for zooplankton groups. The main Fig. 3 Taxa scores on the first two axes of the PCAs of zooplankton taxa densities (log (x + 1)) in the estuaries of Bilbao and Urdaibai. Abbreviation meanings can be found in Table 2. NIS copepods (red); A. bifilosa, C. aquaedulcis, and other brackish copepods (green); meroplankton groups (yellow); neritic groups and copepods (blue); freshwater copepods (black); and other categories (gray) differences between periods in the seasonal patterns of alpha diversity of copepod species and zooplankton groups were due to the decrease of diversity in spring from 1998-2002 to 2010-2015, which was more marked in copepods at intermediate salinities (U30 and U33).
The relationship between deseasonalized variations in copepod alpha diversity and deseasonalized density variations of copepod species differed between estuaries and salinities (Fig. 5). In the estuary of Bilbao, diversity showed a strong negative relationship to the density of the NIS A. tonsa at all salinities and also to the density of the NIS O. davisae in the innermost salinity zones (B30 and B33). Likewise, the density of P. marinus was found to influence copepod diversity positively at B30 and the density of C. aquaedulcis negatively at B33. In the outer salinity zones (B34 and B35), however, diversity was strongly and negatively related to the density of the neritic copepod A. clausi, but positively related to the density of the PCPC-calanus assemblage. The density of the congeneric species O. nana and O. similis species also showed positive relationships to copepod diversity at some salinity zones of this estuary. In contrast to the estuary of Bilbao, in the estuary of Urdaibai, copepod alpha diversity showed a strong negative relationship to the density of A. bifilosa at all salinity sites, except at U35, where it was inversely related to the densities of A. clausi and PCPC-calanus. In this estuary, copepod alpha diversity showed a negative relationship to the density of the NIS A. tonsa only at U26 and U30. Regarding zooplankton groups (Fig. 6), their alpha diversity was found to be negatively related to the density of copepods at all the salinity sites under study in both estuaries. In addition, in the estuary of Bilbao, alpha diversity of zooplankton groups showed positive relationships with the densities of cirripede larvae, gastropod larvae, and appendicularians at B35 and with those of cladocerans and appendicularians at B33 and B34.

Gamma and Beta Diversities
Gamma diversity values and patterns of variation for copepod species and zooplankton groups differed between estuaries (Fig. 7). Copepod gamma diversity was, in general, slightly higher in the estuary of Urdaibai, and its seasonal pattern was less variable along the study period in the estuary of Bilbao, where it showed lowest values in late winterearly spring and highest values in autumn. A slight trend of increase of copepod gamma diversity was observed in both estuaries over the study period, but the increase of diversity occurred mainly in late winter-early spring in the estuary of Bilbao and in autumn in the estuary of Urdaibai. In contrast, zooplankton gamma diversity was higher in the estuary of Bilbao than in the estuary of Urdaibai, and the seasonal pattern of zooplankton diversity was similar in the two estuaries, showing lowest values in autumn-winter and highest values in spring-summer. A clear decrease in zooplankton gamma diversity was observed from the first (1998)(1999)(2000)(2001)(2002) to the last (2010-2015) period in the estuary of Bilbao, but not in the estuary of Urdaibai.
Deseasonalized copepod gamma diversity was negatively related to the deseasonalized density of the neritic species A. clausi and the NIS A. tonsa in the estuary of Bilbao, but showed a strong negative relation to that of the brackish autochthonous species A. bifilosa in the estuary of Urdaibai. In contrast, it showed a strong positive relationship to the deseasonalized density of the neritic species O. nana and O. similis in both estuaries, of O. media in the estuary of Bilbao, and the PCPC-calanus assemblage in the estuary of Urdaibai (Fig. 8a). The deseasonalized gamma diversity of zooplankton groups was strongly and negatively related to the deseasonalized density of copepods in both estuaries, but it was also positively related to the density of cladocerans and appendicularians in the estuary of Bilbao (Fig. 8b).
Beta diversity (Fig. 7) increased both for copepod species and zooplankton groups from the 1998-2002 period to the following periods in both estuaries, but the increment was higher in the estuary of Bilbao. However, from 2003-2009 to 2010-2015, both beta diversities increased in the estuary of Bilbao, while they showed a slight decrease in the estuary of Urdaibai. The seasonal pattern of beta diversity for copepod species and zooplankton groups differed between years and periods in both estuaries. However, the most obvious variations of the seasonal pattern of beta diversity between periods was observed for Fig. 6 Dispersion plots of deseasonalized density (log (x + 1)) of key zooplankton groups vrs. the deseasonalized zooplankton groups alpha diversity (Shannon index) for each salinity site of the estuaries of Bilbao and Urdaibai. Straight lines show the tendency of the relationships. Cop, copepods; Cir, cirripede larvae; Gas, gastropod larvae; Cla, cladocerans; and App, appendicularians copepods in the estuary of Bilbao, where beta diversity values showed the most marked increase in summer.

Discussion
Overall, results revealed that the NIS copepods detected during the 18-year-study had a much higher impact on zooplankton community structure and diversity in the estuary of Bilbao than in the estuary of Urdaibai. This impact was highest in the uppermost zone and decreased with increasing salinity in both estuaries, being negligible at the highest salinity sites located at the lower estuary. Therefore, our results agree with other observations in which the lower salinity habitats of the upper reaches, characterized by a minimum of native species richness, are the most sensitive to biological invasions in estuaries (Nehring 2006). However, more importantly, they showed that the impact of arriving species may differ largely between estuaries due to a complex interaction of factors, including the history of human perturbations which may determine the suitability of a system to be colonized by new species.

Changes in Zooplankton Structure
In the estuary of Bilbao, the occurrences of the copepod NIS A. tonsa and O. davisae, and that of C. aquaedulcis, and to a lesser extent the occurrences of A. bifilosa and the  (Frisch et al. 2006), also affected the native zooplankton community structure. While in the lowest salinity habitat, the main mode of variation (PCA, axis 1) indicated a progressive increase of zooplankton in the presence of the brackish NIS and A. bifilosa and C. aquaedulcis, in intermediate salinity waters of 34, a slight return towards the initial community structure was observed, this denoting a landward increase in the influence of neritic taxa, to the detriment of the seaward advance of brackish species. This fact may be related to the changes in water quality during the study period, as we can infer from the increase of dissolved oxygen levels (Villate et al. 2013;Iriarte et al. 2016). Aravena et al. (2009) found that oxygen levels affected the proportion of A. tonsa to A. clausi densities at intermediate salinities of the estuary of Bilbao, because the density of A. clausi increased at higher oxygen levels, whereas A. tonsa showed a competitive advantage over A. clausi under low oxygen conditions due to its higher tolerance to hypoxic conditions (Marcus et al. 2004;Richmond et al. 2006). Some holoplankton and meroplankton groups, mainly appendicularians and larvae of gastropods and bivalves, showed a general interannual trend of increase along the entire estuary, as evidenced by the second main mode of variation of the zooplankton (PCA, axis 2), which supports the idea that the improvement of environmental conditions Fig. 8 Dispersion plots (a) of the deseasonalized density (log (x + 1)) of key species vrs. the deseasonalized copepod species gamma diversity (Shannon index) and (b) of the deseasonalized density (log (x + 1)) of the most abundant zooplankton groups vs. the deseason-alized zooplankton groups gamma diversity (Shannon index) of the estuaries of Bilbao and Urdaibai. Abbreviations of copepod species and zooplankton groups as in Figs. 5 and 6, respectively had a positive effect on zooplankton density in the estuary of Bilbao.
In the estuary of Urdaibai, unlike in the estuary of Bilbao, NIS copepods were not involved in the main changes of the zooplankton community structure, as shown by their low contribution to the PCA axis 1, and no clear differences in the trends of zooplankton change were identified between salinity habitats. The strong interannual variations observed were mainly driven by neritic copepods such as the PCPCcalanus assemblage, O. nana and E. acutifrons. This might be attributable to the stronger hydrodynamics and shallowness of the estuary of Urdaibai which give rise to a lower water column stability and higher sensitivity to climate factors , as it has been observed also in other small estuaries of southern Europe (Vieira et al. 2015). The strong changes in zooplankton observed in 2012 at all the salinity sites in the present study were already reported for the outer marine zone of both estuaries and were suggested to be linked to a hydro-climatic event (Fanjul et al. 2017). Our results suggest that such event could have affected the zooplankton structure along the entire estuary and that it had a smaller effect on the zooplankton of the estuary of Bilbao.
Brackish copepods, mainly A. tonsa and A. bifilosa, together with some meroplankton groups, were the taxa that contributed most to the second mode of zooplankton variation in the estuary of Urdaibai. Before the colonization by A. tonsa and O. davisae, the estuary of Urdaibai, unlike the estuary of Bilbao, had the typical brackish zooplankton community with a dominant species (A. bifilosa), and the dominance of A. bifilosa seems to have limited both the abundance (and relative abundance) of the congeneric NIS A. tonsa and the length of the season in which it is abundant Barroeta et al. 2020). As the NIS A. tonsa is functionally similar to the indigenous species A. bifilosa (David et al. 2007), no novel impact on zooplankton community that could otherwise cause ecosystem change happened, as found in other works too (see Doherty-Bone et al. 2019 and references therein). It is interesting to point out that meroplankton groups like gastropod and polychaeta larvae showed a similar level of contribution to that of the dominant brackish copepods to the changes in the zooplankton of the estuary of Urdaibai. In this system, gastropod larvae share salinity habitats with brackish copepods and may become the dominant taxa in the zooplankton assemblage of the inner estuary in summer (Villate et al. 1993;Villate 1997).

Spatial Changes in Both Estuaries
Overall, the alpha diversity of both copepod species and zooplankton groups decreased with decreasing salinity from the mouth to the head over the entire study period in the estuary of Urdaibai and since the establishment of a dominant brackish assemblage comprised by the NIS copepods in the estuary of Bilbao. This pattern of change in diversity with salinity is in agreement with the relationship between these two variables depicted in the Remane diagram (Remane 1934) for the salinity range analyzed in our study and has also being observed for copepods and zooplankton communities in other estuaries (e.g., Grindley 1981;Duggan et al. 2008;Whitfield et al. 2012). Therefore, the higher alpha diversity obtained in the inner estuary of Bilbao prior to the colonization by the brackish copepod NIS, with values similar to those measured in the outer estuary, may be considered as an indicator of impaired brackish estuarine communities, since the negligible presence of brackish copepods in the estuary of Bilbao during the first period of study (1998)(1999)(2000)(2001)(2002) allowed the occupation of the inner reaches of this estuary by a lower abundance but higher diversity assemblage of neritic zooplankton .

Interannual Changes in the Estuary of Bilbao
In the estuary of Bilbao, the main interannual changes in copepod alpha diversity occurred in the innermost zone and were mainly explained by A. tonsa and O. davisae density variations. The strong drop of diversity values from 2002-2003 to 2004-2005 coincided with the sudden increase of A. tonsa and to a lesser extent of O. davisae and the decrease of neritic species such as A. clausi (Barroeta et al. 2020). A similar case that relates a remarkable decrease in copepod diversity to a decrease in density of some species at the same time that other species increased in dominance was reported for a site in Arcachon Bay (Bay of Biscay) by Richirt et al. (2019). However, the progressive increase in diversity that occurred later, with highest values in the last period  in the estuary of Bilbao, was mainly due to the arrival and increase in abundance from 2008 to 2010 of C. aquaedulcis and P. marinus (Barroeta et al. 2020). This is an interesting finding, because it showed that the previous establishment of NIS that became dominant (A. tonsa and O. davisae) did not hinder the successful colonization of the inner estuary by other new arriving species, and suggests that this estuary is able to receive more NIS or range-expanding species which could contribute to increase even more the diversity in the future. The progressive colonization of the inner estuary by NIS and C. aquaedulcis had also an evident effect on the beta diversity, because prior to the occurrence of NIS in large numbers, the community was more similar at all salinities due to the dominance of neritic species along the entire estuary, but the recovery of an estuarine community dominated by brackish species increased species heterogeneity at the estuarine system level. Dispersion plots showed that A. tonsa contributed considerably to the changes in copepod gamma diversity in the estuary of Bilbao, but that neritic species were even more influential. Copepod gamma diversity, however, showed much smaller variations than alpha and beta diversities. In some cases, although a modification in species numbers or densities affected the alpha diversity of some sites, gamma diversity did not change because other sites maintained the species pool (Bonecker et al. 2013). In our case, the pool of neritic species in the estuary of Bilbao was much richer in species than the pool of colonizing brackish NIS copepods plus A. bifilosa and C. aquaedulcis, which did not affect the structure of the neritic assemblage in the lower estuary.

Seasonal Changes in the Estuary of Bilbao
The seasonal patterns of alpha diversity in zooplankton groups and copepod species at the inner estuary of Bilbao were also substantially modified by the progressive copepod colonization because NIS peaked in summer-autumn and C. aquaedulcis in spring-summer (Barroeta et al. 2020). Before the occurrence of NIS, the seasonal pattern of alpha diversity was similar at all salinities, being the typical seasonal pattern of the neritic community, with low values during the first part of the year (from January to June) and an increase during the second part (from June to December), also observed in other temperate coastal regions of the North Atlantic (Johnson et al. 2011). When A. tonsa and O. davisae became very abundant and dominant since 2003 (Barroeta et al. 2020), the diversity of copepods decreased during the second part of the year, in a similar way to observations in summer in the Uruguayan Solís Grande estuary (Gómez-Erache et al. 2000) and Doñana Park (Spain) artificial ponds (Frisch et al. 2006) due to A. tonsa increases. When A. bifilosa and C. aquaedulcis occurred in large numbers from 2010 onwards (Barroeta et al. 2020), the diversity of copepods increased in the first part of the year, reflecting a more estuarine type of seasonal pattern of diversity (Soetaert andVan Rijswijk 1993, Primo et al. 2009), and similar to that found in the inner estuary of Urdaibai. Seasonal differences in copepods' gamma diversity diminished over the study period mainly due to the progressive increase of diversity in the first half of the year from one period to the next. Copepods' diversity minima in the first part of the year in the estuary of Bilbao is attributable to the strong dominance of A. clausi in late winter-early spring (Fanjul et al. 2016) and the absence of abundant brackish copepods in the inner estuary; the increase in diversity in winter-spring in the entire estuary is attributable to the progressive occupation of the inner estuary by NIS, and finally, by C. aquaedulcis, which peaks in spring (Barroeta et al. 2020). The seasonal change in copepod beta diversity from lowest values in summer-autumn during the first period (1998)(1999)(2000)(2001)(2002) to highest values in the same season during the last period (2010)(2011)(2012)(2013)(2014)(2015) was related to the occupation of the inner estuary by the thermophilic NIS A. tonsa and O. davisae, which comprised a brackish community typical of the warm period, clearly segregated from the species-rich neritic assemblage within the estuary (Barroeta et al. 2020).

Interannual Changes in the Estuary of Urdaibai
In the estuary of Urdaibai, no clear trends in alpha diversity and no remarkable changes in the diversity of the inner estuary were observed after the arrival and settlement of A. tonsa, this denoting that the zooplankton community was resilient to the invasion-induced disruption in terms of changes in diversity. Regarding gamma diversity, the only copepod NIS that successfully established in this estuary, i.e., A. tonsa, clearly made a smaller contribution to the changes in diversity than the congeneric autochthonous brackish species A. bifilosa. However, an increase in beta diversity during the summerautumn period was observed in the estuary of Urdaibai after the occurrence of A. tonsa, whose presence in the estuary was limited to the warm season (Barroeta et al. 2020). Such increase in taxa heterogeneity was in agreement with the observed spatial segregation of the population of A. bifilosa, which was displaced seaward after the colonization of the inner estuary by A. tonsa . The seasonal and spatial segregation of these two congeneric species has also been reported in other European estuaries ( Baretta and Malschaert 1988;Soetaert and Van Rijswijk 1993;David et al. 2007).

Seasonal Changes in the Estuary of Urdaibai
In the estuary of Urdaibai, the observed decrease in alpha diversity in spring in waters of lower than 33 salinity was associated to the advance of the annual peak of A. bifilosa from summer to late spring after the establishment of the congeneric NIS which had its annual maxima in summer Barroeta et al. 2020). In the Westerschelde estuary, the occurrence of A. tonsa is also limited in time and space, occupying the upstream area in summer, this also causing A. bifilosa (a winter-spring species in this region) to be absent more upstream in summer, affecting the seasonal pattern of A. bifilosa and the copepod community diversity pattern (Soetaert and Van Rijswijk 1993). Overall, all our results revealed that A. bifilosa is the key species that has driven all the changes in zooplankton diversity in the estuary of Urdaibai, although its effect on the increase and seasonal variations of gamma and beta diversities over the study period were attributable to its seasonal and spatial displacements induced by the colonization of the inner estuary by A. tonsa.

The Role of Diversity and Community Structure Prior to the arrival of the New Copepod Species in the Differential Impact of These Species in the Two Estuaries
The diversity of copepod species and zooplankton groups was higher in the inner zone of the estuary of Bilbao than in the inner zone of the estuary of Urdaibai in the period previous to the arrival and establishment of the new copepod species in the inner zone of both estuaries. The colonization success and the resulting impact on the zooplankton structure and diversity of the arriving species were, however, much higher in the first estuary than in the second one. These findings apparently go against the Elton's hypothesis (Elton 1958), which suggests that resident diversity is an important determinant of invasion success because the high diversity increases the competitive environment of communities and makes those environments more difficult to be invaded (Naeem et al. 2000). In our case, however, the higher diversity of both copepod species and zooplankton groups in the inner estuary of Bilbao than in the inner estuary of Urdaibai could not be related to higher ecological resistance to invasions in the sense of the Elton's hypothesis, since the initial diversity measured in the inner estuary of Bilbao corresponded to a lower abundance but higher diversity assemblage of neritic copepod species and zooplankton groups that penetrated the inner estuary when the environmental improvement due to cleaning actions in the estuary allowed so . However, such assemblage of neritic species represents a post-disturbance community unable to compete with brackish invaders for the low salinity niches (Kneitel and Perrault 2006). If we focus on the brackish component of the zooplankton, the initial community of the inner zone of the estuary of Bilbao was less complex and diverse than that of the inner estuary of Urdaibai, the latter being dominated by the brackish species A. bifilosa (Barroeta et al. 2020), and therefore, it may be inferred that the lack of true estuarine competitors rather than the presence of a diverse neritic community determined the invasion success in the estuary. As the lack of typically brackish copepod species in the initial period in the estuary of Bilbao is attributable to the previous environmental deterioration of the inner estuary , our results are consistent with the refined version of the Elton's hypothesis which recognizes that habitat disturbances may facilitate invasions (Hobbs and Huenneke 1992;Reznick et al. 2020).