Shifts in Estuarine Macroinvertebrate Communities Associated With Water Quality and Climate Change

The present work aims to identify changes in the macroinvertebrate community of the Tagus estuary (Portugal) due to improvements in water quality and to climate change. Data was collected over a period of 16 years (1998–2014) from different sites located along the estuarine gradient. The AZTI Marine Biotic Index (AMBI) was used to assess the ecological quality status based on benthic invertebrate communities and identify possible variations associated with changes in water quality. The overall distribution of each species was examined to detect possible changes associated with climate, based on species’ affinity for more temperate or subtropical climates. Results demonstrate that there was an overall improvement of AMBI scores during the assessment period. The analysis of the geographical distribution of benthic species seems to indicate that there has been an increase of species which prefer subtropical climates in the shallower waters of the estuary, whereas in the deeper estuarine sections the propensity is for species that prefer temperate climates.


INTRODUCTION
By straddling the transitional point between rivers and the oceans, estuaries play an essential role in nature and provide a wide variety of essential ecosystem services (Savage et al., 2012;Layman et al., 2014;Gittman et al., 2016). Due to these ecosystem services and sheltered nature, estuaries have attracted human settlements since prehistoric times (López-Romero et al., 2021). Nowadays, these settlements have evolved into great metropolises where millions of people live and work. This rapid urban development in conjunction with the accompanying agricultural and industrial advances puts enormous pressure on estuaries (Lemley et al., 2017). However, due to their naturally high variability in environmental conditions which fluctuate daily, seasonally, and annually, it is very difficult to truly understand the origin of changes within these systems (Yin, 2002;Maes et al., 2004;Cloern and Jassby, 2010;Caffrey et al., 2014;Lanari et al., 2021). Benthic invertebrates are frequently used to assess these changes since they play an essential role in the productivity of the food webs and are also a crucial part of the process of nutrient cycling (van Oevelen et al., 2006). Since species included in this group have low mobility and a close relationship with sediments, they can incorporate biotic responses to a variety of contaminants, such as metals and organic matter, making them excellent indicators for assessing the environmental status of a water body (Nunes et al., 2008).
The AMBI-AZTI Marine Biotic Index (Borja et al., 2000) was developed to assess the ecological status of marine environments in accordance with the Water Framework Directive (WFD). This index has also been used to assess the environmental status and seabed integrity in the Marine Strategy Framework Directive (MSFD), and in general, to determine the status of soft-bottom marine benthic communities based on the sensitivity/tolerance of species to organic pollution. It classifies water bodies into five categories: Undisturbed, Slightly disturbed, Moderately disturbed, Heavily disturbed, and Extremely disturbed. It has been applied under different impact sources, demonstrating its usefulness in detecting specific localized impacts as well as "diffuse pollution." This index has been used broadly all over the world, as a tool to assess many different types of pressures (Borja et al., 2003(Borja et al., , 2019Carvalho et al., 2006;Belhaouari et al., 2019) in very different geographical locations such as the Baltic sea (Muxika et al., 2005), Persian gulf (Hoveizavy et al., 2012), and the Nandgoan coastal waters of India (Sivaraj et al., 2014). Furthermore, it has been adapted for use in United States coastal waters (US AMBI) (Pelletier et al., 2018) and the Arctic (Azovsky and Kokarev, 2019). This worldwide adoption is largely due to its suitability for the assessment of environments such as estuaries (Tweedley et al., 2014;Zhou et al., 2018), coastal lagoons (Carvalho et al., 2006;Reizopoulou et al., 2014), and fjords (Quiroga et al., 2013;Alve et al., 2016).
The soft bottom invertebrates that live in estuaries, are not only affected by human pressures, but also by natural changes such as seasonal and daily (e.g., tidal cycles) variations in the environmental conditions or even the occurrence of extreme stressful events, such as floods and droughts (Chainho et al., 2006). These extreme events have been more frequent because of climate change, a common problem for the planet (Jennerjahn and Mitchell, 2013;Robins et al., 2016). The change of average surface temperature, increased rainfall, acidification, and shifts in marine currents are some examples of the components which may bring on changes in estuarine communities and it is difficult to predict when these changes will occur and on what scale. However, the impact of increased temperature has been proven to be linked to low dissolved oxygen levels and a higher possibility of eutrophication. Increased rainfall has been linked to decreases in salinity which can affect species with a different osmotic regulation capability, while a decreased pH has negative effects on the metabolism of crustaceans and bivalves (Roessig et al., 2004;Pörtner et al., 2008;Bernardino et al., 2016;Goulding et al., 2017;Szalaj et al., 2017;Scanes et al., 2020). Scanes et al. (2020) concluded that effects due to climate change occur at a faster rate in estuaries due to lower depths, which enhances the influence of temperature increase. In a changing world it is of the utmost importance to be able to provide decision makers with the most accurate and up to date information on which to make policy decisions that will affect the lives of their fellow man (Wackernagel et al., 2004). The implementation of conservation and management decisions for ecosystems requires an in-depth knowledge of their biodiversity to inform essential policy (Dubois et al., 2016). Essentially, the longer the time series the clearer the picture, allowing shifts to be analyzed in a historical context by providing a comparative analysis of change of taxa, biomes, and geographic regions while also providing insights into the mechanisms involved (Dornelas et al., 2014;Cloern et al., 2016).
The Tagus estuary is the largest in Portugal, with an area of 325 km 2 and is situated in a biogeographic transition area whose border is located between the Carvoeiro cape in the North and the Tagus estuary in the South (Hayden et al., 1984). This means that this region is already associated with a mixture of flora and fauna from both Atlantic temperate regions and subtropical Mediterranean and Atlantic ones. This makes the Tagus estuary an ideal estuarine ecosystem to study the potential changes in biota resulting from climate change, especially because it exhibits many bays, channels and recesses providing a wide diversity of habitats that favor the colonization by species with different ecological characteristics. However, due to its location in the greater metropolitan area of Lisbon, which has approximately 2.5 million inhabitants, it is under considerable pressure from not only a variety of industrial operations and port facilities but also agricultural operations in the upper reaches . The inner reaches of the estuary present an almost delta-like morphology with low salinity (between 0 and 15), extensive mudflats and saltmarsh areas. The intermediary area is mostly subtidal and more saline (10-30), while the mouth of the estuary is under a high marine influence during most of the year (salinity 30-35) and is much deeper, reaching 40 m Gaudêncio and Cabral, 2007;Vaz et al., 2011). Vasconcelos et al. (2007) found that the Tagus estuary was the most anthropologically pressured in the country by aggregating a multi-metric index which quantitatively evaluates influences from key components: dams, population, industry, port activities, and resource exploitation. Cabral et al. (2012) used the Estuarine Fish Assessment Index (EFAI), to evaluate the ecological integrity of the Tagus estuary's fish communities and ascertains that it had a good Ecological Quality Ratio (EQR) for this element in 2006 and 2010. Chainho et al. (2008) used different benthic assessment tools to understand its robustness in evaluating the benthic ecological status as a response to different levels of human pressure and concluded that most indices and metrics were not sensitive enough to accurately distinguish between all ecological quality classes and to separate between natural and human pressures. Research on the effects of climate change to the biota of the Tagus estuary has been conducted for fish communities (Cabral et al., 2001;Vinagre et al., 2007) however similiar studies for the benthic macrofauna community are non-existent. On the contrary, the use of AMBI index is well documented and has been used to assess the integrity of the ecosystem based on benthic macroinvertebrates (Chainho et al., 2008Azeda et al., 2013).
Due to their geographical characteristics, we can find a much wider range of physicochemical parameters in estuaries than, for example, the open ocean. Therefore, in this environment of extremes we can expect to see faster changes due to climate change. The aim of this paper is to analyze spatial and temporal changes in the composition and structure of the benthic macroinvertebrate community in the Tagus estuary and to identify possible associations between those changes and improvements in water quality as well as influence of climate change on species distribution. Considering that over the last 20 years there have been a variety of interventions within the Tagus estuary to improve water quality, specifically the construction of several Wastewater Treatment Plants (WWTPs), the initial step of this study was to assess the ecological status of the studied sites over time using the AMBI index. The possible influence of climate change was based on Cabral et al. (2001) where they found that there had been an average increase of 1 • C in the Tagus estuary. This was assessed by investigating changes in the community composition based on their native distribution being from either North or South of the ecoregion in which the estuary is located.

MATERIALS AND METHODS
This research was based on data collected from various projects carried out at MARE-Marine and Environmental Science Centre, Lisbon University Campus. Datasets used in this study were selected to address temporal changes in the benthic macrofauna community along different timescales and spatial gradients of the Tagus estuary and relate it to changes associated to (i) improvements on the ecological status of the system and (ii) climate change. Only data from the summer was used to discount seasonal variation and because this is the time of year when the benthic macroinvertebrate communities better reflect the impacts of anthropogenic pressures (Chainho et al., 2007). These samples were collected over 16 years (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) and covered the whole extent of the estuary with a total of 50 sampling stations. From these sampling locations, two types of data were obtained: time series that correspond to yearly samplings in Porto Buxo (PB); Portinho da Costa (PC); Parque das Nações Subtidal (PNS); and Parque das Nações Intertidal (PNI) but are spatially restricted; and spatial series which is restricted to two years (2003 and 2014) but whose sampling stations span the estuary consisting of Estuary Middle (EM), Estuary Bays (EB), and Estuary Lower (EL) (Figure 1).
Yearly samples were taken between 2002 and 2010 at PB and PC sites, using a Day grab with an area of 0.1 m 2 (1 replicate). The samples at PB were collected from 9 sampling stations radiating from the shore. At PC there were 12 sampling stations organized in circular transects approximately 50 m from the shoreline (Figure 1). This area is most like a marine environment, located in the deepest part of the estuary and consists of sandy sediments, with a coarser fraction present at PB and salinity in both sites is usually above 30 (Rodrigues et al., 2016). These samples were collected with the aim of monitoring the effects of the implementation of a municipal WWTP at PC which began operations in 2003, and the discontinuation of discharges of non-treated effluent at PB in the same year (Costa et al., 2018).
The samples from the north bank in the intermediate shallow area of the Tagus estuary were taken yearly, with the aid of a van Veen grab (0.05 m 2 ) (1 replicate) from the PNI (nine stations) and PNS (three stations) sites (Figure 1) between 1998 and 2013. Due to externalities, both sites were not sampled in the year 1999 and neither was PNS in 2002. The sediment is characterized by silts and the salinity varies daily and seasonally between 26 and 31 . Monitoring was carried out to assess the effectiveness of the rehabilitation (conducted until 1998) of the surrounding area at the old docks for the EXPO'98 exhibition carried out in the area currently known as Parque das Nações. Long term spatial changes across the estuarine gradient were assessed based on samples collected at 14 stations distributed randomly in the estuarine area in 2003 and 2014 using a van Veen grab (0.05 m 2 ) (1 replicate). These sampling stations were grouped into three sites according to their location in the estuary: EM (seven stations); EB (four stations); and EL (three stations) (Figure 1). The salinity for these sites can vary between 15 and 30 and as a general rule the sediments range from medium to fine sand in the upstream sites (EM) to silt and clay further downstream (EB and EL) (Chainho et al., 2008).
As mentioned previously, the salinity and temperature within the estuary vary considerably. For the sampling stations located within or closer to the mouth of the estuary (PB, PC, EL) and Estuary Bays (EB), the salinity varies between 32 and 35, while the temperature in PB, PC, and EL has a range around 17-19 • C, the Estuary Bays ranged between 17 and 23 • C. As for locations within the intermediate section of the estuary (PNS, PNI, EM), the salinity would vary between 25 and 32, and the temperature ranged between 18 and 21 • C. As a rule, the sediments range from fine sand (0.0250-0.0063 mm) to silt and clay (<0.063 mm) in the sampling locations located in the intermediate section of the estuary (PNS, PNI, and EM), to mostly fine sands and silt and clay further downstream (PB, PC, EB, EL). The percentage of organic matter across the estuary varies between 4 and 12%. The oxygen saturation in the water is between 81 and 113% for the intermediary section of the estuary (PNS, PNI, and EM), and 91 and 107% for the locations within or closer to the mouth of the estuary (PB, PC, EB, EL).
All sediment samples were fixed and preserved with 4% buffered formalin then sieved using a 0.5 mm mesh, preserved in 70% ethanol and stained with Rose Bengal in the laboratory. The specimens were identified to the lowest possible taxonomic level and then counted to determine the number of taxa and their respective abundance. Hayward and Ryland (1995) was used to identify most taxa, along with other scientific publications on specific taxonomic groups. Taxonomic classification followed WORMS (World Register of Marine Species) (available at www.marinespecies.org) online guidelines.
The AMBI index follows a model (Glémarec and Hily, 1981;Grall and Glémarec, 1997) which categorizes benthic invertebrates into five ecological groups (EGs), depending on their dominance along a gradient of organic enrichment and oxygen depletion. These groups include sensitive species (EG FIGURE 1 | Location of the sampling stations in the Tagus estuary: Porto do Buxo (PB); Portinho da Costa (PC); Parque das Nações Intertidal (PNI); Parque das Nações Subtidal (PNS); Estuary Middle (EM); Estuary Bays (EB) and Estuary Lower (EL). Ocean realms: I (variable westwarard current); II (weak and variable current); IIIw (westward currents); IIIe (strong equatorial currents). Coastal zones: A (subpolar); F (western temperate); G (western subtropical); according to Hayden et al. (1984). I), indifferent species (EG II), tolerant species (EG III), 2nd order opportunists (EG IV) and 1st order opportunists (EG V) (Supplementary Annexe 1). The AMBI value of each sampling station was calculated as proposed by Borja et al. (2000) and a yearly median value for each study site (PB; PC, PNS; PNI; EM; EB; EL) was obtained to understand if there were changes in the benthic ecological status associated with changes in anthropogenic pressures. The AMBI index was calculated using AMBI version 5.0 software, available from https://ambi.azti.es/.
To test the hypothesis that there is a change occurring in the benthic macroinvertebrate community of the Tagus estuary due to climate change, the global distribution of the species located at each study area was ascertained for every year that data was available. Spalding et al. (2007) definition of ecoregions was used, situating the Tagus estuary in the South European Atlantic Shelf Eco Region (SEASER). A classification of 0, 1, or 2 was attributed to each species, based on whether they were a native species to the ecoregions North of the SEASER (type 0), a native species to SEASER (type 1), or a native species to the ecoregions South of SEASER (type 2) (Supplementary Annexe 1). Following this the mean value of the replicates from each area and their deviation from 1 was used to assess whether there was a change in the benthic macroinvertebrate communities and whether this change was influenced by either southern or northern species. Furthermore, it was identified when some of these species started to occur in the estuary and whether their numbers seemed to increase or decrease over time.
A Principal Coordinate Analysis (PCoA) (Clarke and Gorley, 2006;Anderson et al., 2008) was conducted to assess the spatial and temporal patterns of benthic communities at each studied location. The software Primer 6 was used to conduct the ordinations on the triangular matrices of Bray-Curtis similarity coefficients for each sampling station and year. Datasets were previously square root transformed. Furthermore, to identify which species and factors were related to anthropogenic pressures and climate change, only vectors with a Pearson's correlation of at least r < 0.05 with the ordination axes were overlaid on the graphs.
To assess significant temporal differences, a PERMANOVA (Permutational Analysis of Variance) (Anderson et al., 2008) was performed, using a normalized dataset to build a similarity matrix based on Euclidean distances. This test ascertained whether there were significant temporal differences for the AMBI and geographic distribution of species present over the study period. Year and estuarine areas were used as fixed factors in the spatial and temporal gradient of the case study. If significant differences (r < 0.05) were found a paired t-test was performed. Furthermore, G-tests of independence (Sokal and Rohlf, 1995) were run using BIOMstat, version 3.0 (available from appliedbiostat.com) to determine whether distribution of the species (0 or 2) was independent of the sampling year.

RESULTS
A total of 67,188 specimens were identified. Most of the specimens (45,154) were collected in PB (21,557) and PC (23,597) ( Table 1). The dominant taxonomic group at all locations, except for PNI, was Polychaeta in both abundance and richness. At PNI the most abundant taxonomic group was Clitellata, while Polychaeta was the group with the highest richness. Bivalves were the 2nd most dominant taxonomic group in PB and PC. At LE and EB the class Malacostraca was the 2nd most representative taxa ( Table 1). The dominant taxa at PB and PC were Mediomastus fragilis (Rasmussen, 1973) (26 and 23.1% respectively), Pomatoceros sp. (11.4 and 20% respectively), and Oligochaeta (9.5 and 4,8% respectively). At PNS the dominant taxa were Oligochaeta (25%), Streblospio shrubsolii (Buchanan, 1890) (23.9%), and Tharyx sp. The AMBI scores for PNS ( Figure 2C) and PNI ( Figure 2D) do not reveal a change in disturbance levels during the study period and fluctuated between 2.78 and 3.00 in PNI, and 1.47 and 3.01 in PNS. Both reflect SD communities. The AMBI values for EM, EB, and EL ( Figure 2E) demonstrated clear improvements from MD to SD communities between 2003 and 2014. EL went from 3.69 to 3.07, EB from 3.74 to 2.19 and EM from 4.00 to 2.19. The whole estuary data for 2003 and 2014 indicated that there was an improvement in the ecological status over that period with AMBI values decreasing in different estuarine habitats (EM, EB, and EL). This corroborates with the general trend in all the sites during this study demonstrated by the fact that, when compared, the AMBI index is always lower in the last year of sampling compared to the first, with the only exception being PC, where the values remain similar.

Individuals (Types 0 and 2) G-Independence Test Results for Dependency Among Sampling Years
The G-test of independence was run in order to ascertain if the abundance of type 0 and type 2 individuals was independent of sample year for PB, PC, PNS (Table 2), EB, and EM.
For EM and EB independence was found between the two sample years (2003-2014) for individuals type 0 (p < 0.05) and for individuals type 2 (p < 0.05).

Spatial and Temporal Patterns
The first two axes of the PCoA for PB ( Figure 4A) explained 30% of the data variation. It appears that the years 2005, 2006, 2009, and 2010 are separated from the remaining years along the first axis, with these years being positively associated with PCO1. The earlier years (2002 and 2004) were associated to the negative PCO2 axis and AMBI and D were negatively correlated with PCO2, while the earlier years (2002 and 2004) are separated from the others across PCO2. All represented species and AMBI index are positively correlated with PCO1 while geographic distribution (D) is negatively correlated. Similarly, most species are positively correlated with PCO2 excluding Oligochaeta and M. fragilis. On the other hand, AMBI and D are negatively correlated with PCO2. The species vectors are all limited to species which had a Pearsons correlation of r > 0.5. All species represented by the vectors are type 1.
The first two axes in the PCoA for PC ( Figure 4B) account for 27.9% of the data variation. There appears to be a clear separation from the first two sampling years to the last two years along the first axis, with 2002 and 2004 being negatively associated to PCO1 and 2009 and 2010 being positively associated with it. The species vectors are all limited to species which had a Pearson correlation of r > 0.5 with the ordination axes. For PB, the species vectors were all limited to species which had a Pearsons correlation of at least r > 0.5. All the species vectors represented in Figure 4 are positively correlated with PCO1 and associated with the later years of the study (2009 and 2010) except for Nephtys cirrosa (Ehlers, 1868) which is negatively correlated with them. All species and AMBI are negatively correlated with PCO2 and D has a negligible association to both axes.
The first two axes in the PCoA for PNS ( Figure 5A) account for 36.5% of the data variability. Years from 1998 to 2009 were gradually separated along the first axis. The later years (2011, 2012, and 2013) were associated to the negative PCO2 and AMBI and D were negatively correlated with this axis. The species vectors are limited to species which had a Pearson correlation of r > 0.4 with the ordination axes. Most of the species represented in Figure 5A and D are positively correlated with PCO1 with the exception of Corophium sp., S. shrubsolii and AMBI which are negatively correlated. All species as well as AMBI and D are negatively correlated with PCO2.
The first two axes in the PCoA for PNI ( Figure 5B) explain 48.5% of the data variation (31.5 for PCO1 and 17% for PCO2). According to the PCoA, earlier years (negatively associated with PCO1) seemed to be separated from the later years (positively associated with PCO1). The years 2000 and 2008 were negatively associated with PCO2. The species vectors are all limited to species which had a Pearsons correlation of r > 0.3 with the ordination axes. All the species in Figure 5B are positively correlated with PCO1 and PCO2. AMBI is negatively correlated with PCO2, while D is positively correlated, and both are negatively correlated with PCO1. The first two axes in the PCoA for Estuary-EM, EB, and EL (Figure 6) account for 29.6% of the data variation. A separation between 2003, which is negatively associated with PCO1, and 2014, which is positively associated, seemed to occur. The species vectors are all limited to species which had a Pearsons correlation of r > 0.4 with the ordination axes. From the species shown in the Figure 6, S. shrubsoli, Cyatura carinata (Krøyer, 1847), C. setosa, and Boccardiella ligerica (Ferronnière, 1898) and AMBI are negatively correlated with PCO1, while D and Nephtys sp., Parvicardium pinnulatum (Conrad, 1831), N. cirrosa, Heteromastus filiformis (Claparède,

DISCUSSION
The results of the long-term assessment of the benthic communities in the Tagus estuary seem to indicate two main findings: 1. There was an overall improvement in the ecological status over the study period, with AMBI values decreasing in different estuarine habitats (EM, EB, EL, and PB) from MD to SD during the study period. 2. There appears to be a propensity for species which prefer northern climates in the deeper waters at the mouth of the estuary (PB and PC), more southerly distributed species in the mid depths (EM; EL; and PNS) and an absence of either in the shallower more extreme environments (PNI and EB) Studies on this sort of timescale (16 years) regarding benthic macroinvertebrate communities in estuaries are limited (Dauvin, 2000;Frid et al., 2009;Wildsmith et al., 2011) and this can be due to a variety of factors, not least the time required to process samples and the cost involved. Although some other studies using AMBI and MAMBI show that seasonal variability does not influence trends in the data (Borja and Tunberg, 2011), others have demonstrated that there is in fact an effect in estuarine systems where extreme events, such as floods and droughts, are frequent (Reiss and Kröncke, 2005;Chainho et al., 2010). Therefore, it is important to choose seasons in which benthic communities reflect the impacts of anthropogenic pressure on the system rather than natural variability (e.g., seasonal) (Chainho et al., 2007). Another factor which makes the temporal change difficult to measure is the variability of benthic communities along the estuarine gradient (spatial variability) (Chainho et al., 2008;Teixeira et al., 2009). The current study attempts to identify temporal changes in the benthic macroinvertebrate community along different timescales and spatial gradients of the Tagus estuary using data collected for the environmental monitoring of the estuary. This approach identified spatial variations of benthic communities throughout the study area.
At the mouth of the estuary (PB and PC), even though there was a general improvement of the AMBI at PB after the elimination of the non-treated effluent in 2003, no significant changes were noted in the make up of the benthic community during the study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). This feature probably results from the fact that these sites are in the main estuarine channel, with a strong tidal influence, where strong water flows are responsible for the dispersion of pollutants. In fact, even before the elimination of the untreated effluent in PB and construction of the WWTP in PC, the PB area was only moderately polluted and after that it became slightly polluted. There were no significant changes for the AMBI score at the Parque das Nações (PNI and PNS) for the period 1998-2013, mostly because the major improvements in the area had occurred before this time series began (Chainho et al., 2008). The AMBI results for Parque das Nações seemed to indicate that there is a dominance of EG III species which equates to organisms which are more tolerant toward organic enrichment (Grall and Glémarec, 1997). With regards to the spatial series (EM; EB; EL) the results clearly show an improvement in the ecological quality between the two sampling years. Other studies which indicate there has been an improvement of the ecological quality in the Tagus estuary by using other ecological assessment tools  further confirm these results.
The use of benthic indices to measure the state of an ecosystem is fraught with challenges which require the adaptation of existing tools (Chainho et al., 2008) and the use of indicators which integrate a great variety of metrics (Dauvin, 2000;Quintino et al., 2006;Borja et al., 2009;Pinto et al., 2009). The response of these species to variable and unpredictable conditions in transition systems makes it even more difficult to separate the natural and anthropogenic origins of stress on the biota (Elliott and Quintino, 2007). Furthermore, in the current study, the robustness of the AMBI results could be questioned in the Parque das Nações area due to the reduced number of collected specimens, mainly in the intertidal zone, which was predictable due to the poor biological richness that these areas have when compared to subtidal environments. The benthic assemblage of estuaries can change quite quickly due to a variety of environmental factors as well as the potential colonization of alien species. This can severely affect the distribution and abundance of species in an ecosystem (Holt, 1990). For example, between 1930 and 1950, there was a reduction in the abundance of cirripedes on the British coast, because of an increase in water temperature of 0.5 • C (Southward and Crisp, 1954). According to Austin et al. (2001), the number of invasive species in estuaries has increased in the last 20 years due to climate change. These small variations in temperature seem to be enough to affect the abundance of species, especially in systems like estuaries which suffer these effects at a faster rate than the open sea as well as when the organisms are living at the limit of their geographical distribution (Scanes et al., 2020). In the case of the Tagus estuary, the water temperature increased by 1 • C in a 15 year period between 1980 and the mid 90s (Cabral et al., 2001) and this trend is predicted to continue well into the 21st century. In addition, this estuary is located on the border of two biogeographic zones (Hayden et al., 1984), and therefore it would be expected that a change in the macrobenthos community due to these increased temperatures would be reflected in a change of abundances of local species, by increasing the number of species from southern latitudes and reducing the abundance of northern species.
The differences found in species origin and subsequent distribution in the estuary seem to indicate that their geographical distribution is dependent on the depth of different areas within the system. It appears that within this systems shallower areas such as EB and PNI, where the influence of the air temperature is more of a factor, it is too cold in the winter and too warm in the summer for either group of southern or northerly distributed species to successfully colonize these areas and thus allow native fauna to proliferate. In contrast the deeper areas within the estuary such as PNS and EM seem to be more susceptible to an increase in abundance of southerly species, as a gradual water temperature increase and absence of the extremes found in shallower areas (such as PNI) allow some of these southerly distributed species to thrive. Corophium orientale (Schellenberg, 1928) is a good indicator of this kind of pattern. This species was sampled in Parque das Nações during several years between 2005 and 2013, and only in some of the sampling locations that took place in 2003 and 2014 across the whole estuary. C. orientale lives in brackish or estuarine waters, and is found mostly in fine muddy sediments, and occasionally in sandy bottoms (Diviacco and Bianchi, 1987). It is commonly found in the Mediterranean, Black, and Azov Seas (Bousfield and Hoover, 1997), but in recent years, this species has expanded its distribution. It was registered in the Guadalquivir Estuary, in the South of Spain, between 1991and 1994(Cuesta et al., 1996, as well as in southern Portugal (Rodrigues and Dauvin, 1987). More recently was registered in 2003 in the Mira estuary (Chainho et al., 2008). The current distribution of C. orientale may be related to the increase of water temperature within its past range due to climate change, creating the conditions for a natural change of its distribution. This species is very sensitive to variations in temperature, with its mortality increasing at higher temperatures (Bigongiari et al., 2004). This species also seems to adapt well to extreme events associated to climate change, such as floods and droughts, as shown by Gamito et al. (2010), who modeled the evolution of C. orientale populations in the Mira estuary, showing that, these amphipods recover and resume their previous situation shortly after an extreme climatic event, contributing to the resilience of the system.
At the mouth of the estuary where the depth can exceed 40 m the propensity is for the occurrence of more northerly distributed species, and this has not changed since the beginning of the study period. This is in line with the rest of our findings as change in temperature will occur at a much slower rate at these depths, therefore further monitoring over a longer time period would be necessary to see any significant change in species composition in the future. Although there were no significant differences between species with predominantly southerly or northerly distributions over the time series, the latter was always dominant over the former. However, a group of three species that belong to the genus Medicorophium are worth mentioning: Medicorophium aculeatum (Chevreux, 1908), Medicorophium annulatum (Chevreux, 1908), and Medicorophium minimum (Schiecke, 1978). The reason for this is that species belonging to this genus are native to Mediterranean waters, and only in recent years were reported to be found outside of the Mediterranean and Black Seas (Myers and Riera, 2013). M. aculeatum was registered at PB in 2011, 2015, 2016, and 2017 as well as at the Parque das Nações area in 2011 and 2012 (unpublished data). However, it is worth noting that the number of individuals found in PB each year was always one. Endemic to the Mediterranean Sea, this species was registered in 2000-2004 in the Odiel-Tinto estuary (Sánchez-Moyano and García-Asencio, 2010b) in Spain, and near the coastal area of the Guadiana estuary in 2000 on the south-eastern Portuguese border with Spain (Sánchez-Moyano and García-Asencio, 2010a). M. annulatum was registered in the Tagus estuary only at PB and PC sites. This species was found at PB in , and at PC in 2004. According to Myers (1982), this species is usually found in deeper waters, which could explain why the species was only found at the mouth of the estuary. Although endemic to the Mediterranean Sea, details on its distribution are scarce due to the low number of records. The species was also registered at the Sado estuary in 1999 (Carvalho et al., 2001). M. minimum was registered in the Tagus estuary at PB in , and at PC in 2010, 2015. Outside of the Mediterranean Sea, this species has been registered in 2011 in the Canary Islands (Tuya et al., 2014), and in the Algarve (Sampaio et al., 2016).
The change in the distribution of a species is not only conditioned on temperature as there are a plethora of variables that influence a species distribution such as, pH, water quality, dispersion mechanisms, life cycle, biotic interactions, and physical barriers to name just a few (Clark and Frid, 2001;Hiscock et al., 2004;Thuiller, 2004). This is further reinforced due to the demonstration by the G-tests of independence that even though there are significant differences in the number of type 2 individuals in some points of the estuary the same cannot be said for the number of species. This appears to have been the case of the polychaete Cossura coasta (Kitamori, 1960) registered in Parque das Nações in , and in PB in 2009. This species was also found in the samples collected in 2003 and 2014. However, except for the sampling that took place in 2003 where C. coasta was found in considerable numbers in some stations, the number of individuals found across all other locations (Parque das Nações and PB) or years (2014) was always one. C. coasta lives in brackish or estuarine waters and can be found in depths up to 200 m (Jayaraj et al., 2007). This polychaete is native to Indo-Pacific waters, where it can be found from Japan to India. The species was first registered in the Mediterranean in 1983 (Bogdanos and Fredj, 1983), although some claim the identification is questionable (Read, 2000;Zenetos et al., 2005). However, the species was also registered by Chainho et al. (2008) in the Mira and Tagus estuaries, Portugal in 2003.
Changes in climate occur on long time scales and are sometimes hard to truly perceive. Therefore, with a time series of 16 years it is difficult to make any concrete conclusions as the results in the distribution of the species were overall not significant. In Barry et al. (1995) the time series, ranging between 1932 and 1993, revealed significant differences in macrobenthos communities with a shift to more southerly species in Monterey Bay in California with a temperature increase of 0.75 • C. Even with a relatively small time series the same sort of shift could be occurring in the Tagus estuary as the system is also experiencing an increase in temperature (Cabral et al., 2001) and a change in the makeup of benthic biota due to the recent occurrence of non-native species (Moura et al., 2017).
These results show that there has been a change within the benthic macroinvertebrate community in the Tagus estuary. This appears to be related to an increase in water quality and with a shift in species geographic range with an increase in the presence of southernly distributed species in the mid depths of the estuary, the shallower areas being too extreme for the proliferation of either northern or southern distributed species, and no change in the predominantly northerly distributed species in the deeper waters at the mouth of the system. However, a longer time series is necessary to confirm these results, and future studies (from a distribution perspective), should include the winter months as well as summer to ascertain the impact of seasonal variability.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
TG revised and wrote the manuscript. PS conceptualized the study, wrote the original draft, and reviewed manuscript. GS and JM collected samples, performed taxonomical identification, contributed to interpertation of results, and reviewed manuscript. FC performed taxonomical identification. PC and JC conceptualized the study, performed the statistical analyses, interperted the results, and reviewed the manuscript. IM, CF, and NL provided resources from local authorities and contributed for the interpretation of the results. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The authors are thankful to Francisco Ferreira, Miguel Letra, and the MOR and Pérola Branca crews for their support during the sampling surveys. We would also like to thank Sara Cabral, Erica Sá, Maria João Tavares, and Carla Azeda for their assistance during field and laboratory work.