Assessing the effects of temperature and salinity oscillations on a key mesopredator fish from European coastal systems

A population dynamics model was developed to assess the short and long-term effects of temperature and salinity variations in the common goby Pomatoschistus microps in a Portuguese estuary (Minho estuary, NW Portugal). The population was divided into juveniles, females and males, which constituted the model's state variables. Linear regressions between the observed and the predicted density of juveniles, females and the total population were significant. Parameter's sensitivity and uncertainty analysis were estimated. The model was able to satisfactory describe the P. microps population dynamics, and thus was used to simulate the effects of climatic changes on the fish population. Simulations indicated that the common goby population is sensitive to both temperature and salinity changes. Overall, scenarios of more than 3 °C increase caused significant population decreases. Similarly, increased salinities led to a population shrinkage, whereas scenarios of salinity decrease generated an opposite variation on the population. According to the IPCC predictions for climatic tendencies, the population of the common goby will tend to decrease in the near future, experiencing marked oscillations (decrease or increase) during climatic extremes, namely droughts and floods, respectively. These results may be a useful for future planning and management of estuarine systems given that the common goby is an important species of estuarine food webs in many temperate ecosystems.

temperature and salinity regimes will shift in coastal and transitional waters (IPCC, 2014).
These changes are already occurring and influencing all trophic levels, from phytoplankton to carnivorous fishes (Beaugrand, 2009). However, up to date most studies on the effects of climate change on fishes have focused on commercially exploited species (e.g. Hare et al., 2010;MacKenzie et al., 2012). In contrast, studies that focus on small-sized fishes with no commercial interest are still scant. These species have an essential role in the structure and dynamics of food webs because they provide a link between species on lower and higher trophic levels (Doornbos, 1984;Moreira et al., 1992;Cabral, 2000). Their intermediate trophic position within food webs creates a "wasp-waist" flow control, that can be amplified in systems harboring these species in very high densities (Coll and Libralato, 2012;Cury et al., 2000). In a global climate scenario it is expected that temperature and salinity will change in the near future (IPCC, 2014), and for this reason it is important to understand how changes in these two parameters will affect small-sized fish populations that occupy intermediate positions within trophic food webs.
The common goby (Pomatoschistus microps) is a very abundant and widespreaded smallsized fish occurring in temperate estuaries. This species is highly tolerant to environmental constraints, being able to tolerate wide ranges of temperature and salinity, and thrive in harsh environmental conditions (Fonds, 1973;Moreira et al., 1992;Rigal et al., 2008). Despite this several biological and ecological traits of P. microps, such as growth, reproduction, migration and mortality are highly dependent on the conditions of these two parameters (Jones and Miller, Towards this end, a modelling approach was implemented to test the response of the small-sized fish species P. microps under different climatic scenarios. To preclude the possible effects of climate change in the common goby population a system where the species is highly dense and productive was chosen (Souza et al., 2014). This study aims at understanding how a small-sized fish species that occupy intermediate positions within trophic food webs will be able to cope with changes on the climatic conditions.

Study area
This study was conducted in Minho estuary (NW Iberian Peninsula -41º53'N 8º50'0W), which ranges up to 40 km (considering the upstream limit of spring tides), covering a total area of 23 km 2 . This estuary is a shallow system (Moreno et al., 2005), with a mean depth of 2.6 m and maximum width of 2 km (Sousa et al., 2005;Freitas et al., 2009); and is characterized as a mesotidal and partially mixed system, although it tends towards a salt wedge estuary during periods of high river flow (Sousa et al., 2005).
The estuarine fauna is dominated by the European green crab (Carcinus maenas) and the common goby on the epibenthic compartment (Dias et al., 2010;Dolbeth et al., 2010;Souza et al., 2014;Mota et al., 2014), while two non-indigenous species (NIS), the Asian clam (Corbicula fluminea) and the red swamp crayfish (Procambarus clarkii), are the dominant macroinvertebrates in the study area (Sousa et al., 2008a;. In fact, Minho estuary have been invaded by several aquatic NIS in the last decades, which impacted the system in various ways (Sousa et al., 2008b;Mota et al., 2014;Novais et al., 2015;, Ilarri et al., 2015a and females was estimated based on the percentage of contribution of each group within the fifty randomly selected fishes from each sample.
The water temperature and salinity at the bottom were measured each month in the three sampling stations using a multiparameter probe YSI 6820.
The daylight duration data were obtained from NOAA website (http://www.esrl.noaa.gov/ gmd/grad/solcalc/sunrise.html). Monthly daylight duration (minutes) was used as a proxy for the variation in day length at the study site.

Conceptualization and formulation of the model based on P. microps biological and ecological traits
Stage structured models have been proven to be advantageous because they account for different kinetics and parameters that regulate the dynamics and physiology of different life stages of a given species (e.g. juveniles, males and females) (Batchelder and Miller, 1989;Labat, 1991).
In the present model, the estuarine P. microps population was divided into three groups: juveniles, females and males, which are the state variables of the model. The flows between state variables are individuals per unit of time, while the units of the state variables are individuals.100m -2 . The processes that regulate the number of individuals in each group over time are: growth, (the number of individuals transferred from one group to the next), death (the number of individuals subtracted to each group by mortality), migration (the number of individuals subtracted to females and males by the overwintering migration (Jones and Miller, 1966)), and recruitment(the input of juveniles to the population).
The model forcing functions are daylight duration, salinity and water temperature, which affect recruitment, growth, mortality and migration. The model was written in STELLA (Structural Thinking, Experimental Learning Laboratory with Animation) 5.0 software, an object-oriented graphical programming language designed specifically for modeling dynamic systems (Jørgensen and Bendoricchio, 2001), which translates the graphical representation of the model into ordinary differential equations (ODE). The model used a time step (i.e. temporal 156 resolution) of one month to 12 months, chosen to allow a direct comparison with the data obtained on the field (Souza et al., 2014). A simplified conceptual diagram of the model is shown in Fig. 2. The parameters and equations that regulate the number of individuals in each population group are presented in Tables 1 and 2 Calibration refers to the systematic adjustment of model parameter estimates so that the model outputs reflect more accurately the observed dynamic behavior of the system. This procedure is applied when the available information for the parameters is likely to deviate from the normal behavior of the dynamic model. Calibration is a modeling tool often applied when the data for the parameter is adapted from a different system, the population displays heterogeneity and/or is subject to change through time (Beaudouin et al., 2008).

Sensitivity and uncertainty analysis
The sensitivity analysis of the model was estimated for variations of ± 10% on each parameter at a time (i.e. all the other parameters were kept unchanged according to the one-stepat-a-time (OAT) approach). This method explores the parameter space and provides a robust sensitivity measure in the presence of nonlinearity and interactions among the parameters (Wainwright et al., 2014), being widely used in ordinary differential equations models (ODE) due to its simplicity and efficiency.
To estimate the sensitivity of parameters, the following expression (Jørgensen, 1994) was used: Where Y' is the sensitivity of the model outputs to variations on parameters (Xi). The output of any model can be affect by different sources of uncertainty, including input data, choice of parameters, or calibration method, and thus it is important to have an explicit measure quantification of how much the uncertainty affects models outputs (Confalonieri et al., 2016;Roux et al., 2014). To understand how much the model outputs could have been affected by the uncertainty in the measurements of the state variables we performed an uncertainty analysis (UA), on the four model parameters scoring highest in the sensitivity analysis (Table 3).
These parameters are those most likely to affect the results of the model, and thus we assessed accuracy of the model by estimating the relative root square mean error (RRMSE) (Confalonieri et al., 2016). The lower the value of RRMSE, the lower the influence of uncertainty on the outputs of the model and higher is its accuracy. The UA allows to quantify the propagation of uncertainty in the model output that could be caused by natural variation and potential errors associated with the measurement of the state variables (i.e. the density of individuals in each group of the P. microps population) used in the model calibration. To this end, we used a Latin hypercube sampling (LHS) strategy (McKay et al., 1979) to generate a series of virtual observations for each state variable, and assuming a Poisson distribution. The LHS expands the concept of a Latin square for any number of dimensions. The distribution of each variable is divided in equally probable "n" number of intervals (strata) . For each variable a sample is randomly drawn at each interval (McKay et al., 1979), and the values of each variable are then randomly paired to each other. This is a type of stratified random sampling procedure that can be understood as a compromise between a random and stratified sampling techniques that The virtual data series was generated with function randomLHS from the package lhs (Carnell, 2018), the RRMSE calculated with function rrmse from the package Fgmutils (Fraga- Filho et al., 2016), and the RRMSE minimization with function optim from the package stats, for R software.

Studied species
The common goby Pomatoschistus microps is a widely distributed estuarine species spanning ca. 44º in latitudinal range, occurring from Norway to Mauritania, including the Canary Islands, western Mediterranean and Baltic Sea (Froese and Pauly, 2016). This species is often reported as one of the most abundant fish in northern Atlantic estuaries (Martinho et al., (Nyman, 1953;Jones and Reynolds, 1999;Pampoulie, 2001). Males fertilize the eggs, fan and guard them until hatching (Svensson et al., 1998;Jones and Reynolds, 1999;Pampoulie, 2001).
During nest guarding behavior, males often prey on their own brood (Magnhagen, 1992) removing ca. 30% of the egg mass of a clutch (Forsgren et al., 1996).
Common gobies are known to have a high individual fecundity (Bouchereau and Guelorget, 1998), with each mature female being able to generate from 460 to 3400 eggs (Miller, 1986;Bouchereau et al., 1989;Bouchereau and Guelorget, 1998), but the mortality rate during the early stages of fish development is also very high (Leis, 2007). In fact, the survival rate of marine and diadromous fish larvae varies between 6.7 x 10 -5 and 0.1% (Dahlberg, 1979). No information was found in the literature regarding the mortality rate of P. microps larvae, and due to this, the value of larvae mortality used in the model was obtained through calibration and based on values of other marine and estuarine fishes.
The percentage of mature females on the population during spawning season was estimated as the ratio between the number of females in advanced stages of gonadal development and the total number of females.
The lag between spawning and recruitment was established in one-time-step, given that the species has a very short larval phase (2 to 10 days), and the recruitment likely occurs shortly after (Riley, 2003).

Mortality
One of the most important shortcomings in the knowledge of estuarine fishes is the lack of estimates on the source of mortality for any life history stage (Houde, 2008). Even where mortality estimates have been made for estuarine species, the influence of confounding factors (i.e. gear avoidance, inaccessible habitats, etc.) makes it difficult to determine mortality rates (Able and Fahay, 2010). As far as we know, there is no published paper addressing the mortality rate of the common goby in nature, therefore mortality rates used in the model were obtained differently, we assume that their mortality rate is also different (Magnhagen, 1992;Svensson et al., 1998).
Notwithstanding, mortality rate may vary throughout the year in temperate estuarine fishes (Able and Fahay, 2010). In fact, during winter, small and relatively immobile fish experience an increase in their mortality rates, due to net energy deficits caused by low temperatures and food scarcity (Sogard, 1997;Hurst et al., 2000;Hales and Able, 2001;Hurst, 2007). This may lead to an increase in the mortality of estuarine fish of about 33% during winter (Able and Fahay, 2010). The seasonal variation of P. microps mortality was taken into account in the model by assuming an increment of 30% in the mortality rate of all population groups when water temperature was decreases below 10ºC.
The number of P. microps individuals subtracted to each population group was defined by: ( Where Mortalityi = mortality of the population group i; MortRatei = mortality rate of the population group i; Di = density of the population group i.
The parameters values and the equations of each population groups are presented in Tables 1 and 2, respectively.

Migration
The typical life cycle of the common goby lasts for one year, with adults migrating to warmer waters during winter (Jones and Miller, 1966;Muus, 1967). Given that the common goby presents a dynamic and plastic behavior in several life traits Pampoulie et al., 2000;Heubel et al., 2008), it is expected that the temperature level which triggers seasonal migration in estuarine populations should also be different across the geographical range of the species (Jones and Miller, 1966). For instance, Jones and Miller (1966) reported that migration is triggered when temperature is lower than 7ºC, while Claridge  (1985) mentioned that at 5ºC migration is triggered. For other estuarine overwintering migrating species in nearby systems, it is argued that 10ºC is responsible for triggering seasonal migration (Gomes, 1991). Given the scarcity of information regarding the temperature level that triggers overwintering migration of common gobies in Southern European estuaries, we have considered reasonable to assume that temperatures lower than 10ºC induce P. microps migration in Minho estuary.
Moreover, migration can also be triggered by other environmental cues, such as precipitation, drought, water discharge and photoperiod . In a recent study, McNamara et al. (2011) suggested that photoperiod is probably the most prominent and universal variable, indicating that the time of the year can also be relevant to several organisms.
Photoperiod is a reliable indicator of the time of the year, and thus, can be a useful predictor of the phenology of resources . In this context, the photoperiod was also taken into account in the migration equation. The number of migrating P. microps individuals in each population group was defined by: Where Migrai = migration of the population group i; MigraRatei = migration rate of the population group i; Di = density of the population group i.

Effect of temperature and salinity on growth
Since P. microps is not able to control its body temperature to a significant degree, the typical response is that its metabolic rate varies directly with ambient temperature (von Oertzen, 1983). The common goby presents a relatively wide tolerance range for temperature variation, and is able to cope with temperatures ranging from -1ºC to 24ºC (Fonds, 1973;Moreira et al., 1992). Freitas et al. (2010) assumed that the optimal temperature for P. microps growth is 20ºC, but since this value could not be experimentally validated, the value used in the model was Salinity is one of the most important environmental factors affecting the growth and survival of aquatic organisms, influencing both physiological and ecological processes (Poizat et al., 2004;Nordlie, 2006), and many studies have demonstrated the influence of external salinity on growth capacities of fishes (Boeuf and Payan, 2001). The metabolic rate of P.
microps varies directly with salinity (Rigal et al., 2008), and the species has a relatively wide tolerance range for salinity variation, withstanding salinities ranging from 0 to 51 psu (Rigal et al., 2008), though better physiological performances occur at low salinities (Pampoulie et al ., 2000;Rigal et al., 2008). To cope with this, the model uses an optimum curve to describe the effect of salinity on P. microps growth. The optimum salinity value for the species (SOpt) was obtained in the literature and then by calibratied (see Table 1) The effect of temperature and salinity on P. microps growth was described as an optimum-type curve (Martins et al., 2008), where: (4) Where i = temperature/salinity; iopt = optimum temperature/salinity for growth; imin = minimum temperature/salinity at which growth ceases; imax = maximum temperature/salinity at which growth ceases. Long run simulations (240 months) were performed to test the stability of the model.

IPCC predictions
The IPCC (Intergovernmental Panel on Climate Change) Fifth Assessment Report (AR5) predicted that, surface air warming in the 21 st century will range from 1.1 to 6.4ºC (IPCC, 2014). Also, the annual temperature over Europe will warm at a rate of 0.1 and 0.4 Cº per decade, and warming will be greater in southern Europe and northeast Europe (IPCC, 2014).
The IPCC projections show that the annual precipitation will decrease across southern Europe (maximum 1% per decade), resulting in drier summers and wetter winters (IPCC, 2014).
This is likely to cause changes on the salinity levels of estuarine systems, since droughts and floods events will be more frequently triggered in these systems, as recently reported (Cardoso et al., 2008;Dolbeth et al., 2010;Santos et al., 2010;Ilarri et al., 2011).
In this context, several scenarios of temperature and salinity variations in Minho estuary were simulated. Four levels of water temperature increasing (+1, +2, +3 and +4ºC) and four different levels of salinity change (-5 psu, +5 psu, +10 psu, and oscillatory (-5 psu from November to April and +5 psu from May to October)) were simulated. Additionally, the combined effects of temperature and salinity variations were also simulated.
Finally, we performed projection simulations (for 20 years) to assess the extended effects of expected temperature and salinity variations in Southwestern Europe under climatic change scenarios on the common goby population in Minho estuary. Two different rates of temperature increase were simulated: slow (+0.01ºC per year) and rapid (+0.04ºC per year) combined with different scenarios of salinity (normal, -5 psu, +5 psu, and oscillatory). The initial conditions of the simulations followed the conditions measured in the field accompanied by the modification related to the scenarios of temperature and salinity tested in each simulation. The density of P. microps juveniles predicted by the model followed the same pattern as the observed variation, with a marked peak of abundance in December (Fig. 3). The density predicted for females, males and total population also followed similar patterns than those of the observed data, with density continuously increasing after spring and reaching a peak in December or January (Fig. 3).

Model stability, sensitivity and uncertainty analysis
The model showed long-term stability, which supports the internal logic of the model (Jørgensen, 1994). The sensitivity analysis identified the parameters related to reproduction (egg loss, fecundity, larval dispersal and mortality, mature females) to be the most sensitive (Table 3).
The uncertainty analysis carried out through the recalibration of the most sensitive parameters with the data series generated by the LHS delivered a very narrow frequency distribution of the RRMSE values, ranging from 17.10% and 18.20%. Most of these were lower than that obtained with the empirical data (17.82%, Fig. 4). However, the RRMSE value of the empirical model was still low thus well within the range obtained with the virtual data series (Fig. 4), indicating that uncertainty had low influence on the empirical model outputs.

Climatic change simulations
Once the correlation between the model outputs and real data was shown to be satisfactory (Table 4), the model was considered suitable to simulate the effects of the forthcoming climatic changes on the common goby population during a year cycle (12 months).
water temperature scenarios, with the juveniles recruiting earlier in the year in all scenarios except in +1ºC. In +4ºC scenario, the density peak of juveniles will change from December to June, while females and males peaking in July instead of January (Fig. 5).

Salinity variations
Simulations accounting for salinity variations suggest that P. microps population would be benefited by a decrease in salinity (19% increase of the total density in one year), while an increase (+5 and +10 psu) or an oscillatory pattern in salinity would lead to a decrease in P.

Combined effects of temperature and salinity variations
Overall, the combined effects of temperature and salinity increase would lead to a decrease in P. microps population in all scenarios, except the +3ºC combined with a decrease of 5 psu in salinity (19% increase). When oscillatory pattern in salinity is combined with temperature increase fish population would decline only in +1ºC and +4ºC (21% and 33%, respectively), while an increase of 8% and 19% would be observed for +2ºC and +3ºC scenarios, respectively. On the other hand, a temperature increase combined with a salinity decrease would cause a noticeable increase in population levels of P. microps for all scenarios (ranging from 23% to 61%), but it is in +4ºC that the population would decline by nearly 30% in a year cycle. (Fig. 7).

Projection simulations
According to results, temperature increase for longer periods of time would have significant consequences for P. microps population in Minho estuary, with a continuous decrease in population density throughout time in all scenarios (Fig. 8). Similarly an oscillatory pattern of salinity or a salinity increase would lead to a marked decrease in P. microps population, while a salinity decrease would have the opposite effect (Fig. 9). The combined effects of temperature increase and salinity variation will cause an even faster decrease of P. microps density levels in all scenarios accounting for an oscillatory salinity pattern or salinity increases. On contrary, with salinity decreases, the common goby population will initially decrease, recovering after 20 years on the slow IPCC scenario. Conversely, on the rapid IPCC scenario, the population would immediately increase, reaching density values 5 times higher than when compared to the present situation (Fig. 10).

Discussion
The model was capable to satisfactory simulate the variation of P. microps density and dynamics at the Minho estuary. Projection simulations indicated that P. microps population will be highly sensitive to changes in both temperature and salinity. According to predictions, rises in water temperature will cause long-term detrimental effects on P. microps population, with harsher scenarios affecting P. microps more severely.
Furthermore, predictions also suggest that the spawning season might change due to increasing water temperature. In milder scenarios, changes in spawning season might be associated to an extension of the recruitment season, with common gobies spawning earlier in the year. However, in harsher scenarios, the spawning season will be greatly altered, with juveniles starting to recruit in winter but with a marked shortage in the duration of the recruitment season. According to experimental evidence, the duration of spawning seasons has a major effect on P. microps populations (Bouchereau and Guelorget, 1998), and it may be one of the reasons behind the high density of the species in Minho estuary, once in this system, the reproduction season appears to be longer than in other estuaries (Souza et al., 2014). Freitas et al. (2010) assumed that the optimal temperature for P. microps growth is 20ºC.
Nevertheless, previous empirical observations (Dolbeth et al., 2010) and the results from the present model indicated that the species is more abundant and productive at lower temperatures.
However, unless specific experimental studies are conducted to determine the optimal temperature for the growth of P. microps, all other values are assumed and may need to be reviewed in future studies. .

17
The common goby population also responds negatively to salinity increases, indicating that droughts may cause a shrinkage in P. microps populations, which are in line with the results reported by Dolbeth et al. (2010), who observed a decrease in P. microps secondary production after drought events in the same studied site. On the other hand, model outputs suggested that P.
microps population would be largely benefited by flood events, due to the decrease in salinity within the estuary. This agrees with Pampoulie et al. (2000), who described an increased reproductive investment by P. microps after a high freshwater inflow in a coastal lagoon in France. Also, the common goby seemed to be further benefited by the reduction of competitors such as, the sand goby P. minutus, within the lagoon (Pampoulie et al., 2000). Similarly, in a long-term study of P. microps population dynamics Nyitrai et al. (2013) showed that the species peaks in years with higher precipitation, which further suggests that the species is benefited in scenarios of salinity decrease. Notwithstanding, the model showed that the effect of a reduced salinity in winter is voided when accompanied by an increased salinity during summer, suggesting that P. microps populations would decrease in the next years, if the IPCC predictions of wetter winters and drier summers are accurate.
However, it is important to consider that the subsequent effects of climatic extremes may have opposite trends and negative feedback processes (IPCC, 2014). For instance, a massive die-off of bivalves after droughts  and floods (Sousa et al., 2012), may lead to a significant increase on the quantity of empty shells in the river bottom that might be used for P. microps reproduction in the next breeding season, which may led to an increase in the population density after one or more generations due to the persistence of these shells in the system for years (Ilarri et al., 2015b). Actually, the reproduction of common gobies seems to be limited by the presence of nest substrates (Nyman, 1953;Magnhagen, 1998) and their abundance and availability can directly influence the number of breeding males (Breitburg, 1987;Lindström, 1988).
The model was able to predict more accurately the dynamics of juveniles and females, whilst the predicted male dynamics differed more from real data, which may be related to the nest guarding behavior of males, that makes them difficult to be caught within estuaries (Miller, (Bouchereau et al., 1993;Fouda et al., 1993;Koutrakis and Tskliras, 2009). In fact, most of the dissimilarities between the observed and the predicted variation of male density occurred during the breeding season, which supports the idea that male guarding behavior may have influenced the results and lead to such dissimilarities, that are partially because this behavior is not accounted by the present model.
The model was most sensitive to variations in the reproduction parameters. This was somehow expected given that the common gobies present high plasticity on their reproductive traits Pampoulie et al., 2000;Heubel et al., 2008), and suggests that the species can rapidly respond to environmental constrains and rapidly adapt to new environmental conditions. The uncertainty analysis showed that the model output is somewhat sensitive to uncertainty in the measurements of the data used to perform the model calibration (RRMSE = 17.82%). However the range of RRMSE obtained by recalibrating the model with the virtual data series was very narrow, with the RRMSE of the empirical model was well within that range, and thus can be considered to be accurate (Confalonieri et al., 2016).
Despite of the IPCC predictions referring to temperature increase in the air, it should be expected that the water temperature will also increase due to global changes in climate (Bates et al., 2008). Nevertheless, it is unlikely that water temperature will increase at the same rate of the atmospheric temperature, given the differences in the thermal properties between the two fluids; and hence, temperature increase in water probably would be smaller than in air. There are uncertainties in projected changes in hydrological systems since it often depends on a number of variables such as precipitation, evapotranspiration, soil moisture and runoff (Bates et al ., 2008).
In this context, we opted to use the IPCC projections for air temperature increase despite of knowing that the temperature increase in water would be smaller. Nonetheless, it is unlikely that the water temperature would increase as much as the most extreme IPCC scenarios, therefore, the projections on the P. microps population dynamics at +3ºC and +4ºC should be seen with caution and understood as predictions for extreme climatic scenarios. The model predicted that for every tested scenario of temperature increase, the P. microps population would experience a gradual decrease in projection simulations. Also, the most likely scenario of salinity change (oscillatory pattern) in extreme climatic events would lead to a sharp decrease in P. microps density. In this context, it is probable that during the next decades at Minho estuary, common gobies may experience population shrinkage. Given the trophic position and abundance of the species, this could cascade through the estuarine biological community, especially in a system where the species is remarkably abundant such as in Minho estuary (Souza et al., 2014). As a mesopredator, the common goby connects low and high levels of the food-web of fishes (Doornbos, 1984;Moreira et al., 1992;Cabral, 2000); therefore, changes in P. microps population would affect both higher and lower trophic levels, with its trophic role being even more relevant in systems where it achieves higher densities (Pockberger et al., 2014). Nevertheless, the real ecological impact of the P. microps population reduction is hard to predict, since the sympatric species P. minutus may play a similar ecological role (Salgado et al., 2004) to provide functional redundancy (Ives, 1995) and creating an "insurance effect" in the system (Yachi and Loreau, 1999;Loreau et al., 2003). Actually, both species are morphologically and ecologically similar, differing mostly on salinity preferences, with P. minutus preferring to inhabit saltier waters compared to P. microps (Leitão et al., 2006;Dolbeth et al., 2007). Also, both species can often compete for food and space (Złoch and Sapota, 2010) and hence, it is reasonable to assume that P. minutus may perhaps be benefited by a decrease in P. microps population, and potentially fulfill the ecological gaps left by the common goby.
Notwithstanding, given the uncertainty about the ecological effects that a decrease in P. microps density might trigger, it would be interesting to perform further studies on the interactions between P. microps and P. minutus particularly at different conditions of temperature, salinity and density.
Given that the P. microps geographical range of occurrence is wide, and our study was conducted in a system located nearer to the southern edge of the species distribution (Froese and Pauly, 2016), the populations inhabiting systems at higher latitudes and thus subjected to colder further south may suffer more serious consequences. Still, given the plasticity of P. microps, each population may respond differently to environmental changes, and hence, each system should be treated as a unique case of study, despite the trend presented in this study, which predicts a decrease in P. microps density caused by warming waters.
The use of ecological models has been increasing in the last decades, with significant developments in the software tools available and also in their accuracy. Nonetheless, modeling approach still have limitations, which also include the IPCC projections themselves (Hollowed et al., 2013;Cheung et al., 2016). Population dynamics models are widely used but they require a good data set containing homogeneously distributed data. Additionally, the calibration of parameters in population dynamics models are especially difficult (Chatzinikolaou, 2012). The model we developed was tested against a dataset of 12 data entries, which is not a long timeseries for this type of model, but is reasonable enough considering the life cycle of the species, logistic constrains related to the sampling and the time-frame of the project. In addition, the model showed not to be affected by the uncertainty of the state variables and therefore was accurate in its outputs (i.e. low values of RRMSE = relatively high accuracy = model with relatively low susceptibility to uncertainty). Furthermore, it is also important to state that several parameters inputted into the model were obtained from different species and/or localities due to the lack of information in the literature about the common goby and the Minho estuary. These probably influenced the outputs of model, and for this reason, the outcome of our model needs to be seen with caution. Despite of these issues, the robustness of model and its design allowed us to drawn good and cautious interpretations regarding the direction and the magnitude of the shifts in the population dynamics of P. microps.
Our study did not account for limiting factors in the environmental carrying capacity to sustain a population increase of P. microps, therefore the model outputs ought to be seen with the consequences for the ecosystem. For instance, a decrease on the common goby population might have detrimental consequences for the fishery yield, given that the high abundance of P.
microps certainly provides resources for carnivorous fishes targeted by fishermen, but further studies are needed in order to better comprehend the inter-specific responses towards the decline of P. microps population and its consequences for fishery.

Conclusion
The model for P. microps population dynamics seems to be effective in simulating the performance of the common goby in the Minho estuary when submitted to changes in temperature and salinity conditions. The obtained simulations are relevant in the context of the global climate (IPCC, 2014) since they demonstrated that the populations of P. microps in scenarios of temperature and salinity increase responded with a population decrease. While in scenarios of a decrease in salinity, the population will experience a substantial increase in terms of density.
The obtained results presents a projection approach on how a core species will cope with climatic change in the near future. This type of approach represents a useful tool for future planning and management of estuarine systems, once the results predict how P. microps, an important component of estuarine biological communities, will vary with global effects of climate change.

Acknowledgements
The authors would like to thanks Eduardo Martins and Prof. Carlos Antunes for their help during the field campaign in Minho estuary. We would also like to thank Fabiana Freitas, Felipe         psu from November to April, and +5 psu from May to October. Slow scenario = +0.01ºC.y -1 and Rapid scenario = +0.04ºC.y -1 .