Introduction

In the second half of the twentieth century, humans changed the Earth's ecosystems more rapidly and extensively than in any comparable period of time in human history (Millennium Ecosystem Assessment 2005). In particular, continental waters (e.g., wetlands, streams, estuaries) represent the most threatened broad habitat types globally (Dixon et al. 2016; Janssen et al. 2016). As a consequence, it has been estimated that the rate of species loss in continental waters is up to 3.2% of the continental freshwater diversity (Burkhead 2012; and references herein), affecting the whole ecosystems functioning and even evolutionary processes (Naeem et al. 1994; Myers 1996; Lotze 2010). Far from saturation, the rate of species loss is expected to rise in the following decades (Johnson et al. 2017). Specifically, under the ongoing climate change scenario, it could be expected that global warming will push species to the limits of their distribution range, and species loss will be exacerbated, especially at lower latitudes (e.g., the Balkans and the Iberian Peninsula; Jarić et al. 2019).

Aquatic systems are subject to multiple, cumulative, long-term pressures, including pollution, habitat alteration and destruction, and invasive species introduction (Birk et al. 2020; Carvalho et al. 2019; Nõges et al. 2016; Thrush et al. 2020). The European Water Framework Directive (WFD; European Commission 2000) mandates the characterization of water body types to set the biological reference conditions based on normative indicator parameters (e.g. metrics based on measurements of richness, abundance and diversity), and the identification of pressures and impacts to develop classification systems to assess the ecological status of aquatic ecosystems. Aquatic ecosystems not achieving the environmental objectives of the “good” ecological status should have a program of measures implemented to restore them to this status. For instance, climate warning and water abstraction may reduce oxygen concentration in freshwaters and have a negative effect on species richness affecting, in consequence, community structure indicator parameters (Pardo and García 2016; Haubrock et al. 2020) that may consequently deteriorate the ecological status of the water bodies in question. In contrast, improving agricultural practices—e.g. the combined practice of deficit irrigation and fertilization reduction—improved water quality and had a significantly effect on indicator parameters, positively affecting macroinvertebrate communities (Stefanidis et al. 2016), and potentially improve the ecological status. In the case of species introduction, there are some general patterns in their impact. For instance, invasive species play a key role in altering the physical environment and can cause many other problems including genetic pollution, economic losses or ecologic alterations. In general, the capacity of invaders to transform habitats, also has an effect on the abundance of aquatic communities, particularly macrophytes, zooplankton and fish (Gallardo et al. 2016). Consequently, the biotic component in aquatic ecosystem responds to invasive species with taxonomic and functional changes in benthic macroinvertebrates (Fei et al. 2014). However, once established, invasive species are difficult to manage. Hence, in the context of the WFD, there is a need to determine if invasive species should be incorporated as a pressure or as a new indicator metric to be monitored within the macroinvertebrates biological element in the development of future ecosystem classification systems (Cardoso and Free 2008).

Among aquatic ecosystems, estuaries are not only one of the most productive, but also the one most invaded aquatic ecosystems in the world (Costanza et al. 1997; Preisler et al. 2009). Their susceptibility to new invasions has been related to (1) their position as entry gateways for invasive species previously to invade new freshwater environments, (2) their depauperated benthic biota, and (3) their vulnerability due to natural and human disturbances (Cohen and Carlton 1998; Basset et al. 2013). Along with the serious threat posed by invasive species (Williams and Grosholz 2008), the Asian clam, genus Corbicula Megerle von Mühlfeld, 1811, has deeply altered estuarine environments across the globe (Crespo et al. 2015). Ever since Corbicula was found in Europe in the 1980s, most likely, it was simultaneously introduced to major European rivers between 1970 and 1980 (Mouthon 1981; Ferreira-Rodríguez et al. 2020). The species may have been transported in ballast waters via the North America-Europe Atlantic invasion corridor (Pigneur et al. 2014). As a result of its abundance and broad geographic distribution, Corbicula has been particularly well studied in the invaded range [e.g., the works reported in references by Sousa et al. (2008)]. Despite the highest Corbicula densities in temperate regions are reached in estuarine habitats (Ferreira-Rodríguez and Pardo 2014), the ecological effects of Corbicula on the structure of benthic macroinvertebrate communities published to date are mostly focused on freshwaters (but see Franco et al. 2012; Sampaio and Rodil 2014). The effects of Corbicula on these environments are mostly signaled as negatives; e.g., weaker ability of sediments in retaining nutrients due to bioturbation (Chen et al. 2016), decreases in abundance of benthic flagellates and bacteria, and diatoms (Hakenkamp et al. 2001), and potential competition with native freshwater mussels (Ferreira-Rodríguez et al. 2018a, b). Moreover, the role of Corbicula in the provisioning of ecosystem services (or disservices) and functions has been also quantified (e.g., Sousa et al. 2009; Ferreira-Rodríguez et al. 2019b). For instance, Corbicula is able to drive large decreases in phytoplankton and seston in water bodies harboring large populations (Cohen et al. 1984; Leff et al. 1990; Beaver et al. 1991). Along with that, shell accumulation provides hard surfaces for species that prefer structured habitats (Werner and Rothhaupt 2007; Ilarri et al. 2014). In addition, filtration and excretion alter the water column composition by increasing the amounts of ammonium and phosphate, while excretion and mortality increase organic matter and calcium carbonate concentrations in sediments (Ferreira-Rodríguez et al. 2019b). Altogether may influence macroinvertebrate communities and, in consequence, ecological classification systems based on measurements of richness and diversity.

Additionally, in estuarine environments, the gradient of chemical conditions (e.g., water salinity conditions) exerts a strong selective influence limiting local diversity values for some groups of organisms (Telesh and Khlebovich 2010). Furthermore, in a given estuary, some particular traits combination may be more resistant than others to pressures and impacts (e.g., invasive species introduction; Basset et al. 2013). Under this context, the objective of this study was to determine possible changes in macroinvertebrate composition attributable to Corbicula invasion in estuarine environments. For that, macroinvertebrates were sampled along the salinity gradient in northwestern Iberian invaded and uninvaded estuaries (39 sampling sites) where Corbicula can reach vast extensions of thousands of individuals per square meter (Ferreira-Rodríguez and Pardo 2014). We propose two possibilities to analyze the results, the incorporation of Corbicula as a pressure to be managed or the incorporation of Corbicula as a biological element to be monitored in ecological classification systems based on biological elements (e.g., macroinvertebrate communities). First, given that biotic indices based on macroinvertebrates are regularly used for ecological status assessment, if diversity indices reflect changes attributable to Corbicula colonization, we argue that Corbicula should be incorporated as a pressure and, in consequence, the focus of efforts to restore for achieving a good ecological status according to the WFD. Second, as eradication is virtually impossible in systems where Corbicula reach high densities, if multivariate community measures are affected by Corbicula inclusion or exclusion from the analysis, we argue that Corbicula should be incorporated as a biological element to be monitored.

Materials and methods

Site description

This study was performed in northwestern Iberian estuaries in spring 2014. The climate domain of the region is wet oceanic (Köppen classification; AEMET-IM 2011), with tendency to dryness during summer; the mean annual precipitation is around 1400 mm, and the average temperature is 14.3 °C. The region is characterized by an igneous geology and siliceous rivers that conform an extensive hydrographic network (Sabater et al. 2021). The topography determines the existence of short watercourses with river dominance over the estuaries in the estuaries draining to the North into the Cantabrian Sea, and larger watercourses with marine dominance over the estuaries in the estuaries draining to the West into the Atlantic Ocean.

For this study, 60 sampling sites covering 12 estuaries across the region were pre-selected by examining satellite imagery on Google Earth (Google, Mountain View, CA, USA). The first sampling site in each estuary was located at the upper limit of the Juncus maritimus Lamarck zonation (i.e., limit between the euhaline and mesohaline salinity conditions; Mercado et al. 2012). From this site upstream, five sampling sites were established in each water body 1 km apart from each other (Fig. 1). Two additional upstream sampling sites were included in the River Miño to cover the salinity gradient. Some sampling sites were discarded for accessibility reasons (sites considered unsafe or with physical barriers). Finally, a total of 39 sampling sites were included in the study (Supplementary material 1). With the aim to broadly classify these sites, a water category and a water body type were assigned following the WFD implementation in the region. The water bodies’ environmental features and the criteria for their delimitation were extracted from public available online information from river basin management plans (http://mapas.xunta.gal/visores/dhgc/). First, sites within each estuary were classified as either belonging to freshwater (i.e., oligohaline sections) or transitional water (i.e., mesohaline sections) categories. Second, sites in freshwater sections were classified as either Type 21 (n = 10; Cantabric-Atlantic siliceous rivers) or Type 28 (n = 7; Main Cantabric-Atlantic siliceous river axis). Sites in transitional water sections were classified as either Type 8 (n = 15; intertidal Atlantic estuaries with dominance of the river over the estuary), Type 9 (n = 6; intertidal Atlantic estuaries with marine dominance) or Type 10 (n = 1; subtidal Atlantic estuaries). Finally, Corbicula absence (n = 29) or presence (n = 10) was also used to classify each sampling site.

Fig. 1
figure 1

Study area in northwestern Iberian Peninsula. Only main rivers are depicted. Sampled rivers (n = 12) are indicated in boxes

Full size image

Species composition

The aim of this study was to acquire a homogeneously sampled database, to be able to compare samples from rivers with those from transitional waters. We used a methodology already developed and used successfully by the authors in shallow transitional waters of the Balearic Islands, (Lucena-Moya and Pardo 2012) similar to the one used in rivers (Pardo et al. 2014), as many estuaries are shallow environments dominated by macrophytes we thought it was an appropriate sampling. At each sampling site, macroinvertebrates were collected adhering to the WFD standards for Spanish invertebrate sampling (Pardo et al. 2014). Macroinvertebrates were sampled with a kick-net (250 μm mesh size) in littoral zones at low tide, when estuaries were at base flow, each sample containing 20 subsamples collected in a 100 m transect. The subsampled units were selected following a multi-habitat sampling approach in which major habitats (representing > 5% of the total sampled area) were sampled according to their proportional distribution in the transect. Each sub-sample was defined by an area of 0.5 m × 0.25 m, and thus, the 20 subsamples covered a total of 2.5 m2.

All samples were immediately preserved in ethanol 70% or fixed in 4% formalin in the field. Macroinvertebrates were identified in the laboratory to the lowest feasible taxonomic level [usually species and genus with the exception of Acari (Class), Oligochaeta (Family), Diptera (Family, Subfamily or Tribe) and Nematoda (Phyllum)] under a stereo-microscope using dichotomic keys (Tachet et al. 2002). All specimens collected were finally labelled and fixed in 70% alcohol and deposited in the Reference Collection of the Limnology Laboratory, University of Vigo.

Taxa-trait composition

Macroinvertebrate families were classified according to: (1) the functional feeding style (Schmera et al. 2015) and (2) the habit or mode of existence (Heino 2008). As estuaries represent transitional areas between freshwater and marine ecosystems, we arbitrary choose freshwater terminology to define functional feeding styles and habits using specialized literature (e.g., Merritt and Cummins 1996; Tachet et al. 2002) and expert opinion for minor groups of macroinvertebrates. Because we had multiple sources of information, some macroinvertebrates could have been allocated to several groups, but we selected the most common traits signaled in the specialized literature for each taxon. We were able to use this approach to homogenize the variability of traits that often are referred in the literature for a single taxon. First, functional feeding styles are based on morpho-behavioral mechanisms of food acquisition coupled with the ingestion of a wide range of food items and include shredders, collectors-gatherers (collect FPOM from the stream bottom), collector-filterers (collect FPOM from the water column), scrapers, piercers and predators. Secondly, habits characterize taxa locomotion, attachment, concealment, and include burrowers, climbers, clingers, divers, sprawlers and swimmers. Combination of both grouping features (functional feeding style and habit trait) is of high relevance with regard to the functional roles of macroinvertebrates aquatic ecosystems (Heino 2005). Hence, macroinvertebrate families were finally classified into 1 of 23 categories, termed subsequently as functional groups, resulting from the combination of the 2 selected biological traits grouping features (Supplementary material 2).

Physical and chemical measurements

In parallel to the assessment of taxa composition at each sampling site, several in-stream variables that might influence macroinvertebrate assemblages were measured. Dissolved oxygen (mg L−1; DO), oxygen saturation (%; O2), hydrogen potential (unitless; pH), electric conductivity (µs cm−1; EC), temperature (°C; T°), and salinity were measured in situ with portable meters. Oxygen was measured with an oxygen meter YSI 556 MPS, pH with a Thermo Orion 290 + , and EC and salinity with an Orion Model 115 at 25 °C. To analyze water composition, water samples were collected in triplicates and immediately frozen in a portable freezer into polyethylene bottles. Samples were externally analyzed for nutrients and calcium in the Center of Scientific and Technological Research Support (CACTI, University of Vigo). Concentration (μg L−1) of oxidized nitrogen (nitrite N–NO2− and nitrate N–NO3−), ammonium (N–NH4+), and phosphate (P–PO43−) were determined following the specific ISO standard methods for water samples by means of a continuous-flow analyzer (Auto-Analyzer 3, Bran + Luebbe, Germany). Inductively coupled plasma optical emission spectrometry (ICP-OES) was used for the analysis of environmental Ca2+ concentration (mg L−1).

Statistical analysis

Distance-based linear model (DistLM) (Legendre and Anderson 1999; McArdle and Anderson 2001), using first all variables and secondly the best AICc selection criterion (i.e., the model with the lowest AICc value being considered the best), were used to identify the predictor variables which best explain the variation in macroinvertebrate composition and functional groups in the studied estuaries. The predictor variables used in the DistLM routine included N–NO3−, P–PO43−, N–NH4+, Ca2+, pH, EC, and O2. N–NO2− was suppressed from the analyses because of its high correlation with N–NO3− (Spearman coefficient: r = 0.87). Distance-based redundancy analysis (dbRDA) (Anderson et al. 2008) was used to provide a visual representation of the existing relationships between macroinvertebrate assemblages and functional groups and the significant predictor environmental variables.

We performed SIMPER analysis (SIMilarity PERcentage) to detect representative taxa and functional groups (cut-off of 70%) that characterize each water category and water body type, and to estimate the individual contribution and the importance of each taxon to the global similarity between functional groups. The method uses the Bray–Curtis measure of similarity, comparing in turn, each sample in one site with each sample in a second site. The Bray–Curtis similarity is a common index to examine patterns of community composition in ecological status classification systems (e.g., Delgado et al. 2010). Bray–Curtis index operates at the family or functional group level, and therefore the mean similarity between the first and second site can be obtained for each pair of samples. Differences in macroinvertebrate assemblages among the freshwater and transitional water sections were tested with one-way Analysis of Similarities (ANOSIM). Analysis were performed with and without including Corbicula (i.e., Corbicula as a biological element to be monitored). Data from transitional water Type 10 were excluded from the analysis due to the lack of replication. All the aforesaid analyses were performed in PRIMER7 v.7 software (Clarke and Gorley 2015) with the PERMANOVA + add-on package (Anderson et al. 2008).

Finally, alpha-diversity measures (Margalef's index and Shannon–Wiener index) were calculated to assess changes in macroinvertebrate composition related to Corbicula (i.e., Corbicula as a pressure). A one-way Analysis of Variance (ANOVA) test (at a significance level of p = 0.05) was carried out in IBM SPSS (version 19.0; Armonk, NY, US: IBM Corp) to prove whether differences in these indices between invaded and uninvaded estuaries were statistically significant. Levene's test for equality of variances and Kolmogorov–Smirnov’s test to verify the normal distribution of the data had been performed prior the analysis.

Results

Macroinvertebrate composition

A total of 698,229 individuals of aquatic macroinvertebrates belonging to 103 families were collected (Supplementary material 2). In general, Hydrobiidae, Gammmaridae and Corophidae were the most representative taxa (i.e., contribution percentage for average similarity > 10%), and shredders-swimmers and filterers-sprawlers were the representative functional groups. Specifically, in transitional water sections, Corophidae, Hydrobiidae, Gammaridae and Anthuridae were the most represented, with filterers-sprawlers, scrapers-climbers, shredders-swimmers and predators-sprawlers as representative functional groups. In contrast, in freshwater sections, Hydrobiidae, Chironomini and Elmidae were most represented, and scrapers-climbers, gatherers-burrowers and scrapers-clingers were the representative functional groups (Supplementary material 3, Supplementary material 4).

Variation in macroinvertebrate composition across the estuarine environmental gradient

The environmental matrix for the estuaries included the composition and abundance of 93 families, corresponding to 23 functional groups, and 9 environmental variables which ranges are shown in Supplementary material 5. A closer examination of the environmental variables (Ca2+ and EC values) and of their faunal assemblages, together with a first initial dbRDA, indicated that the water category assignation for some sites was erroneous (we used the official classification of water categories in Galicia not adjusted for invertebrate composition), for that four sites we reclassified: site LE5 (Lérez estuary) showed EC values higher than expected for freshwater sections (1787 µs cm−1), thus was reclassified as transitional water section; SO4 from the Sor river (EC of 64.5 µs cm−1), and UL1 and UL2 (EC of 62.5 and 62.4 µs cm−1, respectively) had EC within the range of oligohaline waters in the region, and thus were considered freshwater sections instead of transitional water sections.

Differences in macroinvertebrate assemblages and functional groups were fitted to environmental variables resulting in similar (37% and 39%, respectively) percentages of variance explained. The DistLM model for taxonomic composition including all variables had an R2 = 0.37, but best AICc models reached lower values of explained variance (R2 = 0.27) and anyhow included the same most important variables (EC, P–PO43−, O2 and T°), that were related with the first two dbRDA axes performed with all variables. The first two dbRDA axes with all variables explained 62.9% of total variation of fitted model (Table 1, Fig. 2). The dbRDA1 axis (42.2%) represented the main gradient in EC existing within the estuaries, ranging from the lowest EC values corresponding to freshwater sections on the negative part of the axis, towards most downstream sites in the estuaries with higher EC corresponding to transitional water sections on the positive part (Fig. 2). The dbRDA2 axis (20.8%) identified sites for their higher T° and N–NH4+ contents corresponding mostly to Main Cantabric-Atlantic siliceous river axis towards the negative part from sites from Cantabric-Atlantic siliceous rivers in its positive part. In contrast, the dbRDA2 axis did not differentiate intertidal Atlantic estuaries with dominance of the river over the estuary from intertidal Atlantic estuaries with marine dominance.

Table 1 Proportion of variance in (a) macroinvertebrate composition and (b) functional groups explained by environmental variables (Nitrate nitrogen N–NO3−, Phosphate phosphorous P–PO43−, Ammonium nitrogen N–NH4+, Calcium Ca2+, hydrogen potential pH, temperature T°, electric conductivity EC, and oxygen saturation O2) in forward DistLM sequential tests
Full size table
Fig. 2
figure 2

Distance-based redundancy analysis (dbRDA) ordination plot of invertebrate assemblages composition by sampling sites. Predictor variables selected by the linear model are shown as vectors (Wtemp = temperature, EC = electric conductivity, N–NH4 = ammonium, P–PO4 = phosphorous, N–NO3 = nitrate, Ox_Sat = oxygen saturation, pH = water pH, Ca2+  = calcium). Explained variation in the fitted model and total explained variation are indicated for each axis. Dark circles = freshwater types, empty circles = transitional water types

Full size image

The analyses of the functional groups in the estuaries provided slightly higher values of variance explained as for invertebrate composition. The best AICc model for functional groups (R2 = 0.29) included same most relevant variables as for the model with all variables (EC, P–PO43−, O2 and N–NH4+). The first two dbRDA axes of the analyses performed with all variables explained 71.1% of total variation of fitted model (Table 1, Fig. 3). The dbRDA1 axis (43.2%) represented the most important gradient in EC and P–PO43− that gradually increase towards the right and separated sites in freshwater sections on the left of the ordination from sites in transitional water sections towards the right part (Fig. 3). The dbRDA2 axis (27.9%) identified sites for their higher N–NO3− and T° corresponding mostly to Main Cantabric-Atlantic siliceous river axis towards the negative part from sites from Cantabric-Atlantic siliceous rivers in its positive part. In contrast, the dbRDA2 axis did not differentiate intertidal Atlantic estuaries with dominance of the river over the estuary from intertidal Atlantic estuaries with marine dominance.

Fig. 3
figure 3

Distance-based redundancy analysis (dbRDA) ordination plot of functional group composition by sampling sites. Predictor variables selected by the linear model are shown as vectors (Wtemp = temperature, EC = electric conductivity, N–NH4+  = ammonium, P–PO43− = phosphate, N–NO3− = nitrate, Ox_Sat = oxygen saturation, pH = water pH, Ca++  = calcium). Explained variation in the fitted model and total explained variation are indicated for each axis. Dark circles = freshwater types, empty circles = transitional water types

Full size image

Multivariate community measures (ANOSIM)

The ANOSIM showed differences in macroinvertebrate assemblage composition between freshwater and transitional water sections, with a high dissimilarity between their compositions (Table 2). Similarly, within the freshwater sections, types 21 and 28 showed different macroinvertebrate assemblages. However, invertebrate composition did not differ between types 8 and 9 in transitional water sections. When Corbicula was excluded from the analysis, macroinvertebrate assemblages continued to differ between freshwater and transitional water sections and among typologies in freshwater sections.

Table 2 Summary of SIMPER and ANOSIM analysis results (with Corbicula/without Corbicula) for macroinvertebrate and functional group composition
Full size table

ANOSIM also showed different functional group composition between freshwater and transitional water sections (Table 2). Moreover, functional groups did not differ in composition between the two types 8 and 9 in transitional water sections and types 21 and 28 in freshwater sections. However, when Corbicula was excluded from the analysis, ANOSIM indicated the existence of significant differences in the functional groups representation types in freshwater sections, thus the absence of Corbicula from the community increased the differences in composition between those types. A closer look at the percentage representation of the main five functional feeding groups within type 28 (the one with Corbicula) identified a change in the representation of collector-gatherers for collector-filterers in invaded estuaries (Fig. 4).

Fig. 4
figure 4

Contribution measured as percentage representation of the main five functional feeding groups within invaded and uninvaded main Cantabric-Atlantic siliceous river axis (bars represent standard error)

Full size image

Corbicula invasion yielded to slightly richer and more diverse communities as measured by Margalef index and Shannon–Wiener index (Table 3). However, not significant differences were found between uninvaded and invaded estuaries in alpha-diversity measures (ANOVA: p > 0.05).

Table 3 Average alpha-diversities (± SE) values obtained for taxonomic composition in uninvaded and invaded estuaries by water category and water body type
Full size table

Discussion

Different estuaries and different estuarine sections may harbor specific faunas representative of the biotic and abiotic conditions (WFD; European Commission 2000). Specifically, presumably resulting from changes attributable to Corbicula invasion (i.e., increasing organic matter and calcium concentration; Ferreira-Rodríguez et al. 2019b), there were significant changes in functional feeding group composition in the Cantabric-Atlantic siliceous river main axis type (Type 28). In this regard, habitat alteration is one of the drivers responsible of biotic homogenization attributable to invasive species (Rahel 2000). In uninvaded freshwater sections, collector-gatherers contribution was consistently higher. However, in invaded sections, collectors-filterers significantly increased their contribution at the expenses of collector-gatherers. In general, collector‐gatherers feeds upon amorphous detritus obtained from biofilms (Benstead and Pringle 2004). Moreover, Corbicula presence may impede collector‐gatherers access to soft bottoms and, in consequence, their decline. However, there are contradictory results in the literature. In one side, Selego et al. (2012) also reported collector-gatherers replacement by collector-filterers after river restoration efforts. In contrast, Pereira et al. (2017) reported an apparent association between both functional feeding groups. Hence, conclusions in this regard would be speculative and further manipulative work may be needed to clarify functional feeding groups association.

From the physico-chemical point of view, an estuary represents a boundary of progressive change between freshwater and marine environments (Attrill and Rundle 2002). Present results have shown that, within the same estuary, different macroinvertebrate assemblages are found between water body categories; i.e., freshwater and transitional water sections. In the estuarine salinity gradient, definable macroinvertebrate assemblages (i.e., interacting species populations occurring together in space; sensu Stroud et al. 2015) can be identified at the regional scale. For instance, in terms of abundance, Gammaridae and Hydrobiidae are among the representative families in freshwater as well as in transitional water sections. In contrast, Anthuridae and Corophiidae abundances differentiate freshwater from transitional water sections. Both taxa are sprawlers and wait motionless for food to come within reach (White et al. 1980). However, in the particular case of the filter-feeding family Corophiidae, the great availability of suspended food resources in transitional water sections may allow this family to be found among the dominant families in benthic assemblages.

Taking into account the macroinvertebrate assemblages in northwestern Iberian estuaries, Corbicula invasion slightly increased species richness and diversity. However, these assemblages did not differ significantly. In contrast, at smaller spatial scale, positive correlations were found between Corbicula and macroinvertebrate diversity indices by Ilarri et al. (2012). Such relationships have been attributed to increased availability of detritus (food resource) and shell reefs creation (habitat resource). In agreement with this previous work, our results have shown an increment of the contribution of gatherers (Oligochaeta; which feed upon the benthic fine particulate organic matter) and burrowers (Nematoda; which live within the substrate) in invaded estuaries. However, despite increasing species richness and diversity may be attributable to Corbicula invasion in unstructured estuarine habitats, those effects should not be extrapolate to other systems where native shell reefs can support richer macroinvertebrate assemblages than invasive species (Ilarri et al. 2015). Moreover, it should be noted that our study focuses in one particular habitat—i.e., littoral zone—restricted to the upper limit of the J. maritimus zonation. Hence, a different model could explain the macroinvertebrate assemblage composition if further studies are performed covering the whole cross-section and the entire salinity gradient. In one side, main channels with high sediment mobility may be characterized by fauna tolerant to sediment movement. In these mobile sediments, fauna which require stable sediments, such as burrowing bivalves (e.g., Corbicula; Ferreira-Rodríguez and Pardo 2014), are not likely to flourish due to the potential for smothering and feeding interference. In the other side, a different assemblage composition of salinity tolerant species may characterize the lower estuarine sections (i.e., salinity ranges from 15 to 35; Maze et al. 1993). In addition, it should be taken into account that the methodology used is this study was designed for freshwater systems, not for estuarine systems. This means that most probably the benthic assemblages in the estuaries are not fairly represented. Hence, a specific methodology for sampling the macroinvertebrate community in northwestern Iberian estuaries should be included in the development of future ecosystem classification systems.

Concerning the invasive potential of Corbicula, in lower estuarine sections, the species penetrates in increasingly more stressful conditions to the point where condition become too harsh to survive (e.g., salinity ~ 15 in the River Miño estuary; Ferreira-Rodríguez and Pardo 2016). Changes in the salinity gradient may influence Corbicula distribution because populations may be constrained or overcome the natural barriers that control their successful establishment. At the edge of its distribution, the species is able to survive but not to sustain a self-maintained population (Crespo et al. 2017). However, it should be noted that Corbicula has first invaded intertidal Atlantic estuaries with marine dominance (i.e., Type 9 estuaries). Hence, it is difficult to predict the eventual consequences over the macroinvertebrate assemblages in northwestern Iberian estuaries with dominance of the river over the estuary (i.e., Type 8 estuaries). Based on the invasion history of Corbicula in the region (see Ferreira-Rodríguez et al. 2019a), estuary-river vulnerability seems presumably a function of their hydrodynamic and chemical characteristics; e.g., higher calcium concentration in Atlantic estuarine environments with marine dominance (Ferreira-Rodríguez et al. 2017). Under the current climate change scenario, as summer precipitations decrease (predicted decrease of 9.9%; Schneider et al. 2013), it is likely to expect reduced river flow regimes and greater marine influence in estuaries. In this regard, in Atlantic estuaries with marine dominance, reduced precipitations may constrain the species towards upper parts in the estuary. However, in estuaries dominated by the freshwater influence, reduced precipitations could increase their vulnerability to Corbicula invasion by increasing the marine influence.

Conclusions

The purpose of the current study was to assess changes in macroinvertebrate composition related to Corbicula invasion in northwestern Iberian estuaries. This study has shown that considering Corbicula as a pressure at the regional scale did not provide different results in diversity indices calculation. In contrast, Corbicula inclusion or exclusion from multivariate community measures significantly affected the functional groups representation in the Cantabric-Atlantic siliceous river axis type. Hence, as classification systems are lacking for northwestern Iberian estuaries, these findings suggest that Corbicula should be incorporated as an indicator parameter (i.e. metrics based on measurements of richness, abundance and diversity) or a new metric should be developed penalizing the presence or abundance of invasive species.