Climatic and Eutrophication Effects on the North Aegean Sea Productivity and Anchovy (Engraulis Encrasicolus) Stock

Climate change is expected to have a strong effect in the Mediterranean Sea [1]. Even though climate projections may be characterised by significant uncertainties, current climate model simulations have all indicated a significant warming in the Mediterranean Sea [2-4]. An increase of sea surface temperature is expected due to global warming, while decreased precipitation and river runoff might potentially result in an increase of salinity with an opposite effect on stratification, particularly in coastal river influenced areas. An increase of stratification may reduce the productivity through reduced enrichment of euphotic zone with deep-water nutrients. On the other hand, increased stratification may influence the North Aegean circulation that is dominated by the Black Sea Water (BSW) pathways and Levantine water inflow Figure 1, as well as the thermohaline circulation that is related to the N. Aegean nutrient balance. In particular, increasing stratification results in a weakened thermohaline circulation characterised by decreased Levantine water inflow and southward export of deep N. Aegean waters, caused by the reduced dense water formation. The decreased export of deep nutrient rich N. Aegean waters result in the slight enrichment of the nutrient pool and primary production in open sea areas [5].


Introduction
Climate change is expected to have a strong effect in the Mediterranean Sea [1]. Even though climate projections may be characterised by significant uncertainties, current climate model simulations have all indicated a significant warming in the Mediterranean Sea [2][3][4]. An increase of sea surface temperature is expected due to global warming, while decreased precipitation and river runoff might potentially result in an increase of salinity with an opposite effect on stratification, particularly in coastal river influenced areas. An increase of stratification may reduce the productivity through reduced enrichment of euphotic zone with deep-water nutrients. On the other hand, increased stratification may influence the North Aegean circulation that is dominated by the Black Sea Water (BSW) pathways and Levantine water inflow Figure 1, as well as the thermohaline circulation that is related to the N. Aegean nutrient balance.
In particular, increasing stratification results in a weakened thermohaline circulation characterised by decreased Levantine water inflow and southward export of deep N. Aegean waters, caused by the reduced dense water formation. The decreased export of deep nutrient rich N. Aegean waters result in the slight enrichment of the nutrient pool and primary production in open sea areas [5].
Increasing temperature might also affect the metabolic rates (growth, respiration) of both low and high trophic level (ectothermic) organisms. Small pelagic fish respond rapidly to changing ocean environmental conditions, due to their short life span and plankton-based diet. Therefore, climate induced changes in different abiotic environmental factors, such as temperature and biotic factors such as zooplankton may have a significant impact on the N. Aegean anchovy stock (Engraulis Oceanography & Fisheries Open access Journal encrasicolus). The main anchovy habitats in the N. Aegean are the more productive coastal areas that receive nutrient inputs from river and BSW discharge [6]. Consequently, the N. Aegean anchovy stock is strongly related to the productivity in riverinfluenced areas, which are mainly controlled by river nutrient loads.
The application of a suite of models is a useful tool to predict changes in the abundance of the N. Aegean anchovy stock, if climatic change scenarios will be applied. The POM-ERSEM Lower Trophic Level (LTL) model has been coupled to an anchovy full life cycle Individual Based Model [7], describing explicitly all important processes, such as spawning, growth, natural and fishing mortality and fish movement in order to test climate change simulations using the IPSL-CM4 and SINTEX-G climate forcing [8,9], implementing the baseline "Business As Usual" (BAU) future river loads scenario and computing longterm climatic simulations .

Materials and Methods
The 3D ecosystem model consists of three, highly portable, on-line coupled sub-models: the 3D Princeton Ocean Model (POM) [10], which describes the hydrodynamics of the area, the ecosystem model [11,5] based on the European Regional Seas Ecosystem Model (ERSEM) [12], describing the biogeochemical cycles and the full-life cycle small pelagic fish model developed by Politikos et al. [6] that simulates the dynamics of the North Aegean Sea anchovy stock.
POM is a primitive equation, time dependent, σ-coordinate, free surface, split-mode time step model. It calculates the following equations for the velocity Ui= (U, V, W), temperature T and salinity S.
It contains an embedded second-moment turbulence closure sub model [13], which gives the vertical eddy diffusivity parameters KM and KH. The analogous horizontal parameters FU, FV, FT and FS are calculated through the Smagorinsky [14] formulation. The density (T, S, P) is calculated from the UNESCO equation of state adapted by Mellor [15].
ERSEM is designed to simulate the biogeochemical dynamics of a realistic marine ecosystem [12]. Since its initial development the model has received regular updates and this work makes use the version presented by Petihakis et al. [16]. The model uses a functional type approach to describe the dynamics of the low trophic levels of the ecosystem and the primary producers are split into four phytoplankton functional types. These include three categories based on relative size-picoplankton, nanophytoplankton and dinoflagellates -in addition to diatoms, which are included separately due to their unique incorporation of silicate. Each functional type uses four corresponding variables that represent the content of its constituents: carbon, nitrogen, phosphate, silicate and oxygen. Furthermore, the functional types are characterized by stoichiometric ratios of carbon-to-nutrients, which may vary dynamically between individual ranges defined in the model parameterization by Blackford et al. [17]. The phytoplankton groups differ because of the different values of the parameters characterizing the specific lysis and specific rest respiration.
The zooplankton functional types consist of mesozooplankton, microzooplankton and heterotrophic Nano flagellates. These also separated by relative size because this determines the type of phytoplankton they graze on and the organic matter they excrete. Each zooplankton group is represented by variable corresponding to its biomass (carbon content), as well as its nitrate and phosphate content. Also, the revised version of bacterial sub-model has been used [11]. Pelagic bacteria are assumed to be free-living heterotrophs utilizing particulate and dissolved organic material, produced by the excretion, lysis and mortality of primary and secondary producers as food. Bacterial growth is controlled by the availability of Dissolved Organic Carbon (DOC), by the availability of dissolved organic and inorganic nutrients, which allow them to assimilate DOC and by protozoan grazing [18].
The following equation is solved for the concentration of C for each functional group of the pelagic system: U, V, W represent the velocity field, A H the horizontal viscosity coefficient and K H the vertical eddy mixing coefficient, provided by the POM. ΣBF stands for the total biochemical flux, calculated by ERSEM, for each pelagic group. The equation is approximated by a finite-difference scheme and is solved in two-time steps an explicit conservative scheme for the advection and an implicit one for the vertical diffusion [15,19,20]. The benthic-pelagic coupling is described by a simple first order benthic returns module, which includes the settling of organic detritus into the benthos and diffusional nutrient fluxes into and out of the sediment.
The full-life cycle small pelagic fish Individual Based Model (IBM) developed by Politikos et al. [6] was applied to describe the dynamics of the North Aegean Sea anchovy stock. The fish model was coupled with a lower trophic level (LTL) hydrodynamic/

Oceanography & Fisheries Open access Journal
biogeochemical model [5], providing the zooplankton fields (mesozooplankton, microzooplankton) that are used as prey by the fish, as well as the temperature and ocean currents. In order to investigate the effect of the future climatic conditions on the N. Aegean anchovy stock, long-term climatic simulations  were performed with the full-life cycle anchovy IBM coupled to the hydrodynamic/biogeochemical POM-ERSEM LTL model [16], under the Intergovernmental Panel on Climate Change (IPCC) A1B scenario (balance across all sources -A1B) for CO 2 emissions. The climatic atmospheric forcing was obtained from the Institute Pierre-Simon Laplace (IPSL) IPSL-CM4 global climate model output [8].
Given the significant uncertainties of the climate projections, the climate simulations were repeated with an additional IPCC climatic forcing [9], showing a weaker warming signal (+1.2 o C) as compared to the IPSL-CM4 (+2.7 o C). The simulations mentioned above were performed adopting present river loads in order to isolate the climate effect. An additional simulation, using the IPSL-CM4 climate forcing was performed, implementing the baseline "Business as Usual" (BAU) future river loads scenario that assumes a decreased phosphate load following present trends. In all cases, the open boundary conditions for the hydrodynamics were obtained from Mediterranean basin scale model simulations, using the IPSL-CM4 and SINTEX-G climate forcing, over the same periods. The open boundary conditions for biogeochemical variables were obtained from a hindcast simulation multi-year average of the Mediterranean coupled POM-ERSEM model.

Results and Discussion
In the future climate (IPSL-CM4 global climate model output), the annual mean net heat flux is increased (average=+18W/ m 2 ), particularly in the Eastern Mediterranean, resulting in a significant increase of sea surface temperature (+2.15 ºC, on average) and the decrease of vertical mixing Figure 1. Sea surface salinity (not shown) is slightly increased (+0.09 psu on average) due to the overall increase of evaporation-precipitation (E-P) flux.

Oceanography & Fisheries Open access Journal
Salinity changes are also related to the BSW pathway that appears to acquire a more western direction south of Limnos Island in the 2080-2100 simulation, rather than to the North Eastern Aegean Figure 2. This change in the BSW pathway is related to a reduced Levantine water inflow that results in a southward displacement of the thermohaline front Figure 2.
This effect is part of a weakened thermohaline circulation, characterized by decreased Levantine water inflow and decreased southward export of deep N. Aegean waters due to the reduced dense water formation caused by the increased stratification particularly noticed in Skiros basin, which is a primary dense water formation site [21].

Oceanography & Fisheries Open access Journal
The mean integrated Chlorophyll a (Chl-a), mesozooplankton biomass and net primary production, along with mean Sea Surface Temperature (SST), Sea Surface Salinity (SSS) and mixed layer depth (MLD) over the present and future climate simulations are shown in Figures 3 and 4 respectively. Looking at the differences between the two simulations Figure 5 one can notice the significant increase of the SST (+2.5 C on average). The impact of the temperature increase on stratification (MLD) is mostly noticed in open sea areas, whereas the effect of the salinity increase (0.05 psu on average) due to decreasing river runoff and precipitation is predominant in coastal river influenced areas. Salinity changes are also related to the BSW pathway that appears to acquire a more western direction south of Limnos Island in the 2080-2100 simulation, rather than to the North Eastern Aegean Figure 1. This change in the BSW pathway is related to a reduced Levantine water inflow that results in a southward displacement of the thermohaline front Figure 1. This effect is part of a weakened thermohaline circulation, characterised by decreased Levantine water inflow and decreased southward export of deep N. Aegean waters due to the reduced dense water formation caused by the increased stratification. The decreased export of deep nutrient rich N. Aegean waters (due to thermohaline circulation changes) [5], results in the slight enrichment of the nutrient pool (not shown), particularly in open sea areas, triggering a slight increase of net primary production Figure 5 despite the increasing stratification. Changes in net primary production (netPP) and plankton biomass around Limnos Island are related to changes in the BSW pathway, mentioned above. The reduction of plankton biomass and primary production in coastal areas, particularly in Thermaikos gulf, is related to the increased temperature, affecting the plankton metabolic rates (increased growth, respiration).

Oceanography & Fisheries Open access Journal
As shown in Figure 6, the N. Aegean anchovy stock is decreased by 15 % in the future climate simulation when the SINTEX-G climate forcing is used and considering a temperature increase of +1.2 o C in the future. When the IPSL-CM4 climate forcing is used, with a temperature increase of +2.7 o C, a severe decrease (-80%) of the N. Aegean anchovy stock is observed (projections not shown). The anchovy biomass decrease is related to the zooplankton decrease in coastal areas and more importantly to the future temperature increase that affects the fish metabolic rates [ Figure 5]. According to the bioenergetics model formulation [6,7], the temperature dependence of plankton consumption by anchovy is dome-shaped [22], i.e. it is maximized at an optimum temperature, whereas the respiration increases continuously with temperature Figure 7. The latter results in the decrease of net fish somatic growth at the increased water temperatures of the future. Consequently, egg production (that depends on fish weight) decreases, whereas the early life stages experience higher starvation mortalities at the warmer temperatures due to the increased energy needed to meet maintenance costs.  Although changes in the phytoplankton and zooplankton biomass exhibit a similar behavior to the Net Primary Production, they are affected more by the simulated temperature increase which has a direct influence on growth and respiration rates. The final result is a lower biomass, where a higher temperature increase results in a decreased phytoplankton biomass comparable to that due to decreased river loads. It is interesting to note that although the two mechanisms are very different, one acting on the physiology of the organism and the other on the bottom up control, both produce almost identical effects in terms of phytoplankton biomass. Moreover, in warmer conditions the average size of the organisms in a community would decrease as a consequence of the Temperature Size Rule (TSR) and because smaller organisms have lower absolute energy requirements the number of phytoplankton cells that can be sustained will be higher.

Oceanography & Fisheries Open access Journal
To gain a better understanding of the various factors affecting the anchovy stock in the future climate, the first three years of the simulations, using the IPSL-CM4 forcing, are shown in Figure 8. In order to examine the effect of eutrophication, related to river nutrient loads, another future climate simulation was performed, adopting the baseline "Business As Usual" (BAU) future river loads scenario that assumes a decreased phosphate load following present trends. It may be seen Figure  8 that the three simulations (present, future, future+river BAU) show very similar adult anchovy biomass, initially, after the first year's recruitment (March 1980). For the remaining of the first year, the adult anchovy mean weight is decreased in the future climate simulations, resulting in reduced egg production, which is a function of the fish trophic status/weight. Consequently, the future climate simulations show reduced abundance in larvae and juvenile, which results in a decrease of the next year's recruitment. Figure 8: Total N. Aegean anchovy (Engraulis encrasicolus) adult biomass and mean weight, juvenile/larvae abundance, egg production and larvae/juvenile starvation mortality, in the present (1980-2000, blue lines), the future (2080-2100, red lines) and the future climate simulation with reduced river loads (2080-2100 River BAU, green lines), using the IPSL-CM4 climate forcing. An additional reason for the decrease of the larvae/juveniles is related to the increase of starvation mortality during these stages Figures 8,9 dues to the increased energy needed to meet maintenance costs in the warmer temperatures. After the simulation of the first year, one may see that the mean adult weight is higher in the future climate simulations, which is associated to the relative increase of zooplankton Figure 10, due to the lower predation by the significantly reduced anchovy biomass. Therefore, it appears that the effect of reduced fish predation on the zooplankton is stronger than the effect from climate change and/or the reduced river nutrient loads on the zooplankton biomass.
The decrease in the anchovy stock is stronger in the future climate simulation with reduced river nutrient loads, as in this case there is also a stronger decrease of plankton productivity in river influenced areas, which is particularly reflected in the reduction of the adult biomass/weight, since the adult population is more confined in the coastal river influenced areas, while larvae are more spread in offshore waters. Therefore, the additional decrease in the anchovy stock in this simulation is mostly related to the reduced egg production due to the decrease in the adult weight. As shown in Figure 9, the two future climate simulations show similar starvation mortalities for larvae/ juveniles. Adaptation strategies to climate change in marine systems shows that there is ample evidence in the literature that shows that marine species are adapting to climate change through shifting distributions and timing of biological events, while evidence for adaptation through evolutionary processes is limited [23]. The limits to plasticity in thermal tolerance with future warming of the seas are virtually unknown [24].
The decrease of the anchovy stock in the future climate Figure 6 may be attributed to the reduced zooplankton biomass in coastal river influenced areas Figure 5 that is stronger in the case of reduced river nutrient loads Figure 10 and to the effect of warmer temperatures on the anchovy metabolic rates.

Oceanography & Fisheries Open access Journal
The reduced anchovy growth affects a) The egg production that is linked to fish weight b) The starvation mortality of larvae/juveniles that affects recruitment.
In order to investigate the importance of the two processes, the fractional change of egg production is compared Figure 11, followed by the fractional change in the total number of larvae becoming juveniles and the total number of juveniles that are recruited to the adult biomass [17]. The difference between the egg production change and the number of larvae becoming juveniles can be attributed to larvae starvation mortality, while the difference between the latter and the next year's recruitment may be attributed to juvenile starvation mortality [25]. It may be seen that the stronger relative change is found for the egg production, followed by the effect of juvenile starvation mortality that is increased in the future climate simulations. In the simulation with reduced river loads the effect of reduced egg production appears even stronger [26,27].

Oceanography & Fisheries Open access Journal Conclusion
In this work, a set of linked models has been applied to predict changes in the abundance of the North Aegean anchovy stock using two climatic change scenarios. In the future climate simulation, a stronger decrease of the anchovy stock is found when river nutrient loads are reduced resulting in a decreased plankton productivity, which affects the adult anchovy biomass/ weight. This additional decrease in the anchovy stock, is mostly related to the reduced egg production due to the decrease in the adult weight. Both future climate simulations show similar starvation mortalities for larvae/juveniles. The increased temperature (+1.2 o C till+2.7 o C) in the future climate, results in a decrease in the anchovy stock from 15% to 80%, which is related to the zooplankton decrease and particularly to the effect of the temperature to the fish metabolic rate (increased respiration). In any case, it is evident that a temperature increase in the future climate change will significantly affect the anchovy stock with apparent negative effects in the socioeconomics of the fishery's dependent coastal areas and communities. Developing tools that can describe the most important physical and biochemical processes that combined together determine the dynamics of the ecosystem is therefore essential to predict the effects of climate change in fish populations. Given the complexity of these processes and their interactions, mathematical models can be regarded as unique tools to deliver integrated approaches and better understand the resilience of marine aquatic organisms to climate warming.
The existence of numerical models may efficiently describe the ecosystem and the anchovy dynamics of the N. Aegean Sea, establishes the numerical basis for the development of a forecasting system capable of supporting coastal zone management issues. Such a system will use numerical models in conjunction with observational data and data assimilation techniques.