Relationships Between Northern Hemisphere Teleconnection Patterns and Phytoplankton Productivity in the Neva Estuary (Northeastern Baltic Sea)

Teleconnection patterns can be an important tool for investigating the impact of climate change on biological communities. The aim of the study was, using 2003–2020 data on chlorophyll a concentrations (CHL) and plankton primary production (PP) in midsummer, to determine which of the teleconnection patterns have most pronounced effects on phytoplankton productivity in the estuary located on the border between western and eastern Europe. CHL correlated significantly with the winter values of the North Atlantic Oscillation (NAOw) and Scandinavia (SCANDw) indices, as well as with the values of the annual Polar/Eurasian (POLy) and annual Arctic Oscillation (AOy) indices. PP was significantly correlated with the values of POLy. East Atlantic/Western Russia pattern showed no significant correlation with both phytoplankton indicators. Stepwise multiple linear regressions were performed to determine the most influential indices affecting CHL and PP in the Neva Estuary. POLy, SCANDw, and NAOw appeared to be the main predictors in CHL multiple regression model, while the values of POLy and the July NAO and SCAND values were the main predictors in the PP model. According to our research, the productivity of phytoplankton in the Neva Estuary, located in the most northeastern part of the Baltic Sea, showed a significant relationship with the POL, which determines weather conditions in the northeastern regions of Eurasia. Possible mechanisms of the influence of these teleconnection patterns on phytoplankton productivity are discussed. Using the obtained multi-regression equations and the values of climatic indices, we calculated the values of CHL and PP for 1951–2002 and compared them with the results of field observations. The calculated and measured values of CHL and PP showed a significant increase in phytoplankton productivity in the Neva Estuary in the second half of the 2010s compared to earlier periods. In some years of the 1950s, 1980s, and late 1990s, CHL could also be above average and the low phytoplankton productivity should have been observed in the 1960s–1970s. This indicates a significant contribution of current climate change to fluctuation in phytoplankton productivity observed in recent decades, which should be taken into account when developing measures to protect aquatic ecosystems from eutrophication.


INTRODUCTION
Eutrophication is recognized as a serious environmental and economic problem in coastal areas around the world and in the Baltic Sea in particular (Heisler et al., 2008;Holt et al., 2016;Damar et al., 2020).The problem of the mass development of phytoplankton in recent decades has become especially acute due to climate change, which can aggravate the negative consequences for the water ecosystems from human activities (e.g., Behrenfeld et al., 2006;Doney et al., 2012;Golubkov and Golubkov, 2020;Golubkov, 2021).However, the mechanisms of the eutrophication process as a result of climate change are not well understood (Le et al., 2019).The impact of climate variability on marine and estuarine ecosystem has been extensively discussed in recent decades (e.g., Stenseth et al., 2002;Doney, 2006;Doney et al., 2012;Bogatov and Fedorovskiy, 2016;Le et al., 2019;Golubkov and Golubkov, 2020).Climate fluctuations are exogenous and hidden driving forces that are causing profound large-scale changes in marine ecosystems, affecting the state of their environment and biotic interactions on interannual and longer time scales (Andersen et al., 2011;Doney et al., 2012;Kashkooli et al., 2017;Golubkov, 2021).For instance, interannual fluctuations in air temperature and atmospheric precipitation lead to a change in the runoff of nutrients into the Baltic Sea from the catchment area, as well as to alterations in the composition and productivity of phytoplankton (Kotta et al., 2009;Teutschbein et al., 2017;Golubkov and Golubkov, 2020;Golubkov et al., 2021).
Due to the holistic nature of the climate system, increasing attention has been given to large-scale patterns of climate variability with marked ecological impacts on interannual and longer time scales (Stenseth et al., 2002).It is known that atmospheric circulation and weather conditions in different parts of the world are largely determined by teleconnection patterns, and their fluctuations are reasonably well predicted using climate indices (e.g., Krichak and Alpert, 2005;Baldwin et al., 2007;Popova, 2007;Ceglar et al., 2017;Yao and Luo, 2018).Due to this, a relationship between various biological processes and climatic indices has been found in many aquatic ecosystems (Straile and Adrian, 2000;Blenckner and Hillerbrand, 2002;Weyhenmeyer, 2008;Sharov et al., 2014;Gubelit, 2015;Kashkooli et al., 2017;Greaves et al., 2020;Maximov et al., 2021).Climatic indices are used to predict the productivity of open ocean waters (Le et al., 2019;Greaves et al., 2020) and inland seas (Kashkooli et al., 2017).Teleconnection patterns can also be used for paleoreconstructions (Diz et al., 2018), since there are examples of reconstructing data of their values back to 1675 (Luterbacher et al., 1999) or even to 1500 (Luterbacher et al., 2004).Such studies allow us to look into the past and determine how unique the modern period is, or the same conditions were observed earlier and the modern situation is just a manifestation of a certain cyclical process (Naef-Daenzer et al., 2012;Diz et al., 2018).In addition to this, knowledge of the relationship between the indicators of phytoplankton productivity and atmospheric circulation makes it possible to restore historical data in those waters where there are no field observations.In some cases, these patterns have helped to determine the causes of regime shifts that occurred almost simultaneously in different marine ecosystems.For example, Alheit et al. (2005), determined that synchronous regime shifts in the North and Baltic Seas that occurred in the late 1980s were caused by a change in the phase of the North Atlantic Oscillation teleconnection pattern.
The Neva Estuary is located at the top of the Gulf of Finland and is the most northeastern part of the Baltic Sea (Golubkov and Alimov, 2010).The primary productivity and biomasses of autotrophic organisms in the estuary are high, mainly due to eutrophic effects of the large nutrient inflow from the Neva River, which is the major contributor of freshwater to the Baltic Sea (Golubkov, 2009;Golubkov and Alimov, 2010;Golubkov et al., 2017).Phytoplankton communities of the estuary include freshwater and marine species; eurytopic species also make up a significant proportion (Golubkov et al., 2021).A significant increase in plankton primary production and proportion of freshwater species of algae were observed in the 2010s (Golubkov et al., 2017(Golubkov et al., , 2021)).Apparently, this was due not only to the anthropogenic nutrient load, but also to changes in weather conditions in recent years because of the global warming, which manifests regionally in warm winters and cool rainy summer seasons (Golubkov and Golubkov, 2020).Most regional climate models predict future increases in winter and summer air temperatures and precipitation in the northern Baltic regions (Meier et al., 2012;Teutschbein et al., 2017).Changes in weather conditions affect water temperatures and salinity, nutrient concentrations and plankton primary production (Friedland et al., 2012;Holt et al., 2016;Golubkov and Golubkov, 2020).The changes are especially important for transit waters, including the Neva Estuary, where the water residence time ranges from several days to several weeks (Golubkov and Golubkov, 2020).These ecosystems appear to respond more quickly to changing weather conditions than ecosystems of offshore seawaters (Golubkov and Golubkov, 2020).
The set of teleconnection patterns that affect weather conditions in Europe in general, and the Baltic Sea region in particular, is diverse.The North Atlantic Oscillation (NAO) in the Baltic Sea region affects changes in air temperature, cloudiness and precipitation (Trigo et al., 2002), and water temperature (Girjatowicz and Małgorzata, 2019).This, in turn, affects river flow velocity (Trigo et al., 2004;Uvo et al., 2021), the water level in lakes in the catchment area of the Baltic Sea (Filatov et al., 2016;Plewa et al., 2019) and also on the duration of the growing season and the total plankton productivity in aquatic ecosystems (Belgrano et al., 1999;Straile et al., 2003;Sharov et al., 2014).Changes in the Arctic Oscillation (AO) pattern affect the time of snow cover melting (Schaefer et al., 2004), air temperature (Wang et al., 2005), and amount of precipitation (Zveryaev and Arkhipkin, 2019;Bednorz and Tomczyk, 2021) in northern Europe, as well as the time of the release of the Baltic Sea from the ice cover (Jevrejeva et al., 2003).For Europe, the relationships between the Scandinavia (SCAND) teleconnection pattern and the amount of precipitation (Casanueva et al., 2014), the depth of snow cover (Popova, 2007;Bednorz and Wibig, 2008), and the air temperature in winter (Kerr, 1997;Belgrano et al., 1999;Zveryaev and Arkhipkin, 2019) have been confirmed.East Atlantic/Western Russia (EAWR) teleconnection pattern in the Baltic Sea region affects the air temperature (Ionita, 2014), water temperature in lakes (Ptak et al., 2018), the amount of atmospheric precipitation (Krichak and Alpert, 2005) and the degree of soil moisture (Zveryaev and Arkhipkin, 2019).Polar/Eurasian (POL) teleconnection pattern is closely related to the interaction between the activity of polar vortices in the Arctic moving along 60 • N and mid-latitude circulation over the Asian continent (Barnston and Livezey, 1987).Balling and Goodrich (2011) suggested that the POL is closely related to NAO, and both patterns relate to changes in polar vortex intensity.The POL values correlate with the thickness of the snow cover, the degree of soil moisture, and air temperature in northwestern Russia (Popova, 2007;Zveryaev and Arkhipkin, 2019).Against the background of global warming, the warming trend in the polar region is most pronounced (Bekryaev et al., 2010), and since the Neva Estuary is located on the border of the temperate and subpolar zones (Meteoblue, 2021), changes in POL can be significant for weather conditions and the state of aquatic ecosystems in this region.
The aim of this study was to determine the impact of a number of teleconnection patterns (NAO, AO, EAWR, SCAND, POL) on phytoplankton productivity in the northeastern estuary of the Baltic Sea, located on the border between Western and Eastern Europe.As indicators of the productivity of phytoplankton, we used the concentration of chlorophyll a and the primary production of plankton, which are the popular and important indicators of eutrophication of marine coastal waters (Andersen et al., 2006;Smith, 2007).We tested the hypothesis that due to the eastern location of the Neva Estuary in the Baltic Sea, the set of teleconnection patterns that have a significant effect on phytoplankton productivity will differ from those patterns that are most significant for the regions of Northwestern Europe (e.g., NAO, SCAND).We also tried to determine which teleconnection patterns are most applicable for forecasting and paleoreconstruction of phytoplankton productivity indicators in the estuary.

Study Area
The Neva Estuary, which is located at the top of the Gulf of Finland (Figure 1), receives water from the Neva River, the most full-flowing river of the Baltic region, whose flow averages 2,492 m 3 s −1 (78.6 km 3 year −1 ).The Neva Estuary is located on the border of the subpolar and temperate zones (Meteoblue, 2021), in the northwestern part of Russia near the northeastern border of the European Union (Figure 1).Type of climate according to Köppen-Geiger climate classification (Kottek et al., 2006) is Dfc-Snow climate, fully humid, cool summer.A number of features common to other major Baltic estuaries generally characterizes this estuary.As most of them, the Neva Estuary is brackish-water, non-tidal and shallow.It is the recipient of discharges of treated and untreated wastewaters from St. Petersburg City, which is the largest megalopolis in the Baltic region with a population of more than 5 million citizens (Golubkov et al., 2019).The Flood Protective Facility (Dam) has separated the upper part from the middle part of the estuary since the end of 1980s.It consists of 11 dams separated by broad water passages and ship gates in its southern and northern parts (Ryabchuk et al., 2017).There is no temperature stratification in this upper part of the estuary.Low water transparency, which does not exceed 1.8 m of Secchi depth in summer time, constrains the distribution of bottom vegetation.The slightly brackish-water part of the Neva Estuary is located between the Kotlin Island and a longitude of ca.29 • 10 E. There is a temperature stratification in summer.The general environmental characteristics of the Neva Estuary are presented in Table 1.A more detailed description of the estuary is given in Golubkov and Golubkov (2020).
Phytoplankton in the Neva Estuary is represented by 174 species and forms from eight taxonomic classes (Golubkov et al., 2021).The most diverse and abundant groups of phytoplankton are cyanobacteria, green algae, and diatoms.The green algae and diatoms are dominant groups in the upper part of the estuary east of Kotlin Island.Cyanobacteria dominate in the middle part of the estuary west of Kotlin Island.Phytoplankton assemblage of the Neva Estuary and its relation to environmental factors were described in Golubkov et al. (2021).

Data and Methods
The study analyzed data on the concentration of chlorophyll a and primary production of plankton, which were determined in the course of a long-term study of the ecosystem of the Neva Estuary conducted by the Zoological Institute of the Russian Academy of Sciences.Samples were collected at 17 stations in the Neva Estuary in mid-summer in the period from late July to early August 2003-2020 (Figure 1).Chlorophyll a concentration was determined spectrophotometrically with preliminary extraction with 90% acetone (Grasshoff et al., 1999) and using a C6-multisensor platform with Cyclop-7 submersible sensor (TurnersDesigns, United States).The primary production of plankton (PP) in the water column were measured by the oxygen method of light and dark bottles (Hall et al., 2007;Vernet and Smith, 2007).Detailed description of the method, experimental design, data variability, and relation to environmental factors is given previously (Golubkov et al., 2017;Golubkov and Golubkov, 2020).

Northern Hemisphere Teleconnection Patterns
We used the values of the Arctic Oscillation (AO), North Atlantic Oscillation (NAO), East Atlantic/Western Russia (EAWR), Scandinavia (SCAND), and Polar/Eurasian (POL) indices to determine the relationships between the teleconnection patterns and indicators of phytoplankton productivity.
The AO is a large-scale regime of climate variability, also called the Northern Hemisphere annular regime.AO is a climatic pattern characterized by counterclockwise winds circulating around the Arctic at 55 • north latitude (Thompson andWallace, 1998, 2000).The NAO is based on the pressure difference at sea level between the subtropical (Azores) maximum and the subpolar minimum.The NAO index is often considered a regional manifestation of the AO index   (Thompson and Wallace, 1998).NAO phases are associated with basin-scale changes in the intensity and location of the North Atlantic jet stream and storm trajectory, as well as with large-scale modulations of the normal patterns of zonal and meridional heat and moisture transport (NOAA, 2021).EAWR pattern operates in Eurasia throughout the year (Lim, 2015).The main surface temperature anomalies associated with the positive phase of the EAWR index reflect temperatures below average in large areas of western Russia and north-eastern Africa (Barnston and Livezey, 1987).The POL is closely related to the interaction between polar vortex activity in the Arctic and mid-latitude circulation over the Eurasian continent (Barnston and Livezey, 1987).The SCAND consists of a primary circulation center over Scandinavia with weaker centers of the opposite sign over Western Europe and eastern Russia/Western Mongolia.The positive phase of this model is associated with positive anomalies in the height of the isobaric surface, sometimes reflecting the main blocking anticyclones, over Scandinavia and western Russia (NOAA, 2021).
Monthly teleconnection patterns for 1950-2020 were obtained from the Climate Prediction Center (CPC) database of the National Oceanic and Atmospheric Administration (NOAA) (NOAA, 2021).Their values were averaged over the year (January-December), over the winter months (December-February), and the July values were also used.In this way, a set of fifteen index values was obtained, which were used to determine the relationship between teleconnection patterns and the concentration of chlorophyll a and plankton primary production in the Neva Estuary.

Statistical Analyses and Data Modeling
Statistical analyses were performed using R software (version 3.6.0)(R Development Core Team, 2021).For statistical analysis, data on chlorophyll a concentration and primary plankton production were averaged over 17 stations for each year.Data from 2006 to 2007 were excluded from statistical analysis because of the strong anthropogenic effect on the estuary in those years.The high concentration of phosphorus and low Secchi depth in 2006-2007 was caused by intensive dredging activity related to the construction of a new passenger port and the creation of new lands in the Neva Delta (Golubkov et al., 2008;Golubkov and Alimov, 2010).For primary production, the data for 2005 were also excluded, because in that year its values were not determined from all 17 stations.For consistency, we did not use these data in our statistical analysis.Two matrices were created.One matrix consisting of 16 rows contained chlorophyll a concentrations and climatic indices in different years.The second matrix, consisting of 15 rows, contained the values of the primary production and indices.Both matrices used for the analysis are presented in Supplementary Material.
The pairwise Pearson correlation coefficients between phytoplankton productivity indicators and climatic indices were produced using Microsoft Excel.Stepwise multiple linear regressions were performed to determine the most influential indexes affecting CHL and PP in the Neva Estuary.We have built "model.null"separately for CHL and PP and "model.full"with all indexes separately for CHL and PP.As a result, two models (null and full) was done for CHL and two for PP.Then the "step" function (direction = "both") was used to add and remove indices from model and to find the model with the lowest Akaike Information Criterion (AIC).AIC was used as an estimator of out-of-sample prediction error and thereby relative quality of statistical models for a given set of data.The lower the value, the better for the AIC (Mangiafico, 2021).When the model with lowest AIC was found it called final.model.This procedure was done separate for CHL and PP, and two final models for CHL and PP were chosen and combinations of the most influential teleconnection patterns was found separate for CHL and PP.In the next step, we calculated the R 2 for these final models.After that, we calculated the adjusted R 2 (Adj R 2 ) to avoid the error of overestimating the quality of the models, arising from the addition of each additional variable to the models.Adj R 2 is a version of R 2 , which is calculated to take into account the number of variables in the multiple regressions model.Adj R 2 is always lower than the R-squared.
We then ran an analysis of variance for the individual factors included in the final models by function "Anova" in the car R package (Fox and Weisberg, 2021) to determine which one is most important for prediction of CHL and PP in each final models.The residuals plots were checked for two final models.It showed that models were unbiased, and homoscedastic.Stepwise multiple linear regressions analysis is detailed in Mangiafico (2021).The full script of statistical analyses can be found in Supplementary Material.

Analysis of Pearson's Pairwise Correlations Between Phytoplankton Productivity Indicators and the Indices of Teleconnection Patterns
Pearson's pairwise correlation analysis between the indicators of phytoplankton productivity in the estuary and the values of climatic indices shows the presence of statistically significant positive and negative correlation coefficients (Table 2).The concentration of chlorophyll a (CHL) positively correlated with the annual average values of the Arctic oscillation (AOy) and negatively with the values of the annual Polar/Eurasia pattern (POLy) (Table 2).The annual average values of North Atlantic oscillation (NAOy), East Atlantic/Western Russia (EAWRy), Scandinavia (SCANDy) patterns did not significantly correlate with CHL.The values of the second phytoplankton productivity indicator, the primary production of plankton (PP), did not positively correlate with any of the annual values of the indices, but negatively correlated with POLy, and the relationship with the latter was the strongest and the p-value was even < 0.01 (Table 2).
Summer CHL also synchronously positively correlated with the winter value of the Arctic (AOw) and North Atlantic (NAOw) oscillations, with the strongest and most reliable correlation with NAOw (Table 2).An even stronger but negative relationship, with a p-value of < 0.01, was between CHL and the winter mean of the Scandinavia teleconnection pattern (SCANDw).Midsummer chlorophyll concentrations in the Neva Estuary were higher in years when AOw and NAOw values were higher and SCANDw values lower.Plankton primary production in midsummer did not statistically significantly correlated with any of the winter averaged indices (Table 2).Similarly, CHL and PP did not significantly correlate with any of the July indices values (Table 2).
Overall, the largest number of statistically significant correlations CHL showed with the winter averaged various indices.Therefore, AOw, NAOw, and SCANDw indices are best used to predict or reconstruct midsummer CHL values in this estuary.On the contrary, the plankton primary production correlated closely and reliably with POLy from of all the indices used.
Since chlorophyll a concentration and plankton primary production correlated differently with different climatic parameters, we analyzed the relationship between these parameters.As shown by regression analysis, these two parameters were related, the Pearson's correlation coefficient between them was statistically significant, but not very high (Figure 2).

Stepwise Multiple Regression Analyses Between Phytoplankton Productivity Indicators and Teleconnection Patterns
Analysis of Pearson's pairwise correlations showed that the concentration of chlorophyll a correlates differently with various climatic indices.On the other hand, all these teleconnection patterns are part of a single atmospheric circulation system.Therefore, it is also important to look for the relationships between the indicators of phytoplankton productivity and the set of teleconnection patterns.To achieve this goal, the method of stepwise multiple regression analysis was used.According to the graphs of the residuals of the predicted values of CHL and PP (Figure 3), both models are unbiased and homoscedastic.For the CHL model, the values of three of the fifteen climate indices we used, NAOw, SCANDw and POLy, were significant predictors (Table 3).The coefficient of determination was high (R 2 = 0.55).However, in order to avoid the error of overestimating it, which arises as a result of adding each additional predictor to the model, we calculated the adjusted R 2 (Adj R 2 ), which also showed a high value of 0.44.This means that the resulting regression model is in good agreement with the observed CHL values.This is also evidenced by the high level of significance (p-value) of the model, 0.0196, i.e., the chances of getting a statistically significant result with this model were high.Two predictors, POLy and SCANDw, had a negative relationship with CHL values, and NAOw had a positive one.The resulting multi-regression equation describing the change in CHL values in the Neva Estuary depending on significant climatic indices is shown in Table 3.
To clarify the significance of the predictors and to find which of them most strongly affect the concentration of chlorophyll a in the estuary, we performed Student's t-test and Fisher's F-test.These tests showed, the addition or removal of each of these three predictors did not affect the reliability of the model (Table 3).It was important that the values of these three predictors change synchronously together.
For the primary production model, three of the fifteen predictors were also significant (Table 4).The resulting PP model had a higher linear determination coefficient (R 2 ) and adjusted R 2 than the CHL model, 0.74 and 0.67, respectively.In addition, the p-value of the PP model was an order of magnitude higher (p = 0.0014) than for the CHL model (p = 0.0196).As in the case of CHL, PP was negatively associated with the average annual values of the Polar/Eurasian pattern.Values of this pattern were the most significant predictor in the model, as shown by the t and F-tests (Table 4).As for the CHL model, the values the NAO and SCAND patterns were also additional predictors of the PP  model, not for the winter months, but for July (NAOj, SCANDj), i.e., for the last month preceding the sampling period.At the same time, as t and F-tests showed, a change in each of the predictors separately affected the reliability of the model, in contrast to the model obtained for CHL.The resulting regression equation is given in Table 4.

Reconstruction of Chlorophyll a Concentration and Plankton Primary Production Based on Teleconnection Patterns
Using the obtained multi-regression equations and the values of the climatic indices from 1951 to 2002, we calculated the CHL and PP values for this period and plotted them on the diagrams together with the data of modern observations (Figure 4).The figures also show the calculated mean values of trophic indices for each of the obtained CHL and PP data series.The charts show CHL and PP values as positive or negative deviations relative to the calculated long-term average values for each indicator.
The diagrams show that the values of phytoplankton productivity indicators in second part of the 2010s are significantly higher than the average long-term ones (Figure 4).The calculation using the multi-regression equation did not show such high CHL concentrations before 2000s as in second part of the 2010s, although, in some years of the 1950s, 1980s, and late 1990s, CHL could also be above average.According to our observations, CHL was noticeably below the average from 2003 to 2013, and based on the multi-regression equation, this situation could prevail throughout the second half of the twentieth century.In accordance with the teleconnection patterns, a long period with CHL values below the long-term average value could be from the late 1950s to the early 1980s, and in the 1960s, the concentration of chlorophyll a could be almost two times lower than its minimum values in the 2000s (Figure 4A).The values of plankton primary production, like CHL, were higher than the long-term average in 2016-2020, although the same high value was in 2004 (Figure 4B).The calculations showed that similar high PP values could be observed, for example, in 1985 and 1998.Based on the obtained regression, the 1980s and most of the 1950s as a whole could be characterized by PP values higher than the long-term average (Figure 4B).As in the case of CHL, PP was below the long-term average for much of the 2000s, with a minimum in 2007.According to the values of the climatic indices, the low PP could have been in 1963, 1968, 1976, and 1997, but not as low as in 2007.A fairly long period with PP values below the long-term average could be during the entire 1970s (Figure 4B).

DISCUSSION
Teleconnection patterns can be an important tool for investigating the impact of climate change on biological communities because of the holistic nature of the climate system (Stenseth et al., 2002).Moreover, these patterns allow to reconstruct their possible ecological characteristics in the past (Naef-Daenzer et al., 2012;Diz et al., 2018).The analysis of Pearson's pairwise correlations between midsummer chlorophyll a concentration and plankton primary production in the Neva Estuary and the values of climatic indices showed the presence of statistically significant relationships between them (Table 2).However, these two indicators of phytoplankton productivity showed different relationships with various climatic indices.Although these two parameters are linearly related (Figure 2), the correlation coefficient between them is lower than one.The reason is that the primary production is a more physiologically mediated indicator (Andersen et al., 2006;Smith, 2007), depending not only on the concentration of pigments in algae, but also on the illumination, which changes rapidly depending on cloudy or sunny weather.Windy or calm weather can also affect the ability of algae to photosynthesize through the depth of water mixing.With weak wind mixing of the water, the algae sink to a greater depth, where the irradiation is lower (Reynolds, 1999).In other words, the concentration of chlorophyll is a more conservative indicator, which responds more slowly to changes in environmental conditions, compared to primary production.In addition, phytoplankton biomass and the concentration of chlorophyll in water depends on the sedimentation rate of phytoplankton, which differs with different phytoplankton composition (Spilling et al., 2018) and its grazing by zooplankton (Reynolds, 2006).Various groups of phytoplankton in the Neva Estuary in recent decades have shown different trends associated with fluctuations in weather conditions (Golubkov et al., 2021).As a result, the response of different characteristics of phytoplankton productivity to changes in weather conditions may differ.Therefore, when determining trends in eutrophication, it is better to use different indicators of this process including primary production (Andersen et al., 2006;Smith, 2007).
Studies of the relationship between the productivity of marine phytoplankton and the NAO pattern are the most numerous in comparison with other teleconnection patterns that we have used in our analysis (e.g., Belgrano et al., 1999;Villate et al., 2008;Boyce et al., 2010;Gladan et al., 2010).In our study, the winter NAO and SCAND indices showed statistically significant Pearson's correlations with phytoplankton productivity (Table 2) and acted as additional predictors in the multiple regression model for chlorophyll a concentration (Table 3).In 2010s, the NAOw index was mostly in a positive phase, while the SCANDw index showed no definite trend (Figure 5).Moreover, the SCAND in the winter months have fluctuated over the past 70 years and stable positive or negative phases were not very pronounced (Figure 5B).In contrast, since the 1970s, the positive phase of the winter NAO index prevailed, and its positive values increased even more in 2012-2020 (Figure 5A).
Multiple regression models obtained on the basis of 2003-2020 data allowed us to roughly reconstruct the changes in CHL and PP values in the Neva Estuary over the past 70 years, provided that their changes were controlled by the same set of factors as in the modern period.Such a reconstruction, probably, does not give a completely objective assessment of changes in these indicators in the past, since it does not take into account the influence of other environmental factors, for example, anthropogenic ones.However, it allows one to assess the state of the system before the onset of rapid climatic changes and highlight possible past periods with high or low CHL and PP values.The analysis shows that the second half of 2010s are characterized by the values of CHL and PP higher than the long-term average (Figure 4).This is especially true for CHL, as the calculation using multiple regression did not show the possibility of such high CHL values in the past.This result is consistent with the available direct measurements of chlorophyll a and plankton primary production in the Neva Estuary in the 1980s and 1990s.For example, the concentration of chlorophyll a in the middle part of the estuary was 10.57 ± 1.66 mg m −3 , and the plankton primary production was 0.72 ± 0.29 gC m −2 day −1 in August 1996 (Telesh et al., 1999), which is close to their values calculated using our multi-regression equations (Figure 4).The average concentration of chlorophyll in this part of the estuary for 1983-1995 was 13.5 ± 1.04 mg m −3 (Silina, 1997), which is also close to the values according to the multi-regression equation (Figure 3).
The mechanisms of the influence of the NAO pattern on CHL and PP are very diverse.With a positive NAO phase, stronger winds are also observed in the Baltic, which leads to a stronger mixing of the water column and to the supply of additional nutrients to the surface waters (Neumann and Schernewski, 2008).In addition, with a positive NAO phase, there is an increase in the surface runoff flow from the Baltic Sea and the bottom counterflow from the North Sea (Jędrasik and Kowalewski, 2019).This can provide an influx of nutrients from the rich bottom layers to the surface waters.
A positive NAO phase and a negative SCAND phase means that warm winters with little snow cover prevail in Scandinavia and other parts of Northern Europe (Kerr, 1997;Belgrano et al., 1999;Popova, 2007).Warm winters in this region lead to an increase in the runoff of nutrients from the catchment area as a result of strong leaching of soils due to their less freezing and winter precipitation in the form of rain rather than snow (Kotta et al., 2009;Teutschbein et al., 2017;Golubkov and Golubkov, 2020).In turn, an increased runoff of nutrients can lead to an augmentation in the primary production of plankton and the concentration of chlorophyll in water.For example, warm winters led to an increase in primary production in Swedish fjords located in the Kattegat straits connecting the Baltic and North Seas; for these fjords, positive correlation coefficients were obtained between the primary production of plankton and the value of the winter NAO index (Lindahl et al., 1998;Belgrano et al., 1999).An increase in the productivity of summer phytoplankton in years with warm winters was also observed in the Neva Estuary (Golubkov and Golubkov, 2020).
In contrast to the northern regions, in the coastal part of the Adriatic Sea in southern Europe, the productivity of phytoplankton was positively correlated with NAO only in winter, but not in warmer periods, when it was negatively correlated with the values of this teleconnection pattern (Gladan et al., 2010).Unlike the Neva Estuary, which freezes in winter, photosynthesis of plankton in the Adriatic goes on all year round and is especially intense in winter, when water masses are mixed as result of strong winds during a positive NAO, as well as due to a decrease in temperature and an increase in the density of surface waters.On the contrary, dry and hot weather prevails in the Mediterranean in summer during the positive NAO phase (Hurrell et al., 2003).This increases stratification and reduces mixing of the water column and the supply of nutrients from deep layers to the surface.In spring, NAO also negatively correlated with plankton primary production, because, due to the prevalence of rainy and cloudy weather, the amount of light falling on the sea surface decreased (Gladan et al., 2010).In addition, an increase in precipitation in the European Mediterranean does not lead to an increase in the runoff of nutrients from the catchment area, because soils and river waters are rich in carbonates (Viličić et al., 2009), and their high concentrations lead to the deposition of phosphorus from solutions (Rozan et al., 2002).A negative relationship between NAO and phytoplankton productivity was also shown for the open waters of the North Atlantic, where algal productivity decreased in the positive NAO phase (Boyce et al., 2010).When NAO is positive, the region is dominated by westerly winds, which increase the temperature of the ocean surface waters and the authors hypothesize that this leads to an increase in the biomass of zooplankton, whose press reduces the biomass of phytoplankton (Boyce et al., 2010).
There are also two opposite opinions about the reasons for the current decline in phytoplankton productivity in the equatorial and subequatorial zones of the oceans.Some researchers suggest that this occurs due to the stronger stratification of oceanic waters as a result of an increase in the temperature of the upper water layer (Boyce et al., 2010;Moore et al., 2018).While others believe that modern conditions lead to a change in the strength and direction of winds affecting ocean currents and the delivery of nutrients for phytoplankton to certain zones of the oceans (Lozier et al., 2011;Whitney, 2015).These examples show that, depending on the geographic location and characteristics of the water system, as well as its catchment area, the same global patterns of atmospheric circulation can affect phytoplankton in different ways.Therefore, more research from different regions and types of waterbodies is needed to understand how modern climate change leading to temperature changes (Boyce et al., 2010;Moore et al., 2018), and the intensity and directions of atmospheric circulations (Lozier et al., 2011;Whitney, 2015) affect the productivity of phytoplankton.
The July values of the NAO and SCAND indices, which turned out to be predictors for the PP model (Table 3), are usually not used to determine the relationship between phytoplankton productivity and teleconnection patterns.Previously, it was shown that the negative SCAND phase in summer leads to lower temperatures and more precipitation in the northwest of the Russian Federation (Zveryaev and Arkhipkin, 2019).On the contrary, the positive SCAND phase in summer causes an anticyclone over Scadinavia and the appearance of an upwelling along the Polish Baltic coast due to the prevailing northeasterly winds (Bednorz et al., 2013).The SCAND pattern also strongly affects the level of the Baltic Sea, because the cyclonic center over Scandinavia, during the negative SCAND phase, enhances the western circulation, which triggers the inflow of water from the North Sea through the Danish straits into the Baltic Sea and can cause surges from the west in the Gulf of Finland (Bednorz and Tomczyk, 2021).Since the water from the Neva River flows from east to west, with a strong westerly wind, there is a surge from the open part of the Gulf of Finland and a backwater of constantly incoming river waters.This lead to the flood and increases the residence time of water in the Neva Estuary, which leads to the accumulation of phytoplankton biomass (Golubkov, 2009).Thus, weather conditions with a positive NAOj phase and a negative SCANDj phase, on the one hand, lead to an increase in the amount of precipitation and the inflow of nutrients, and, on the other hand, to the prevalence of westerly winds, which cause an increase in the water residence time in the estuary.Apparently, owing to these two mechanisms, the productivity of midsummer phytoplankton positively related with July NAO values and negatively with SCAND.This is especially true for plankton primary production, the multi-regression model of which uses NAOj and SCANDj as additional predictors (Table 4).
Pearson's pairwise correlation analysis and stepwise multiregression analysis showed the decisive role of the Polar/ European teleconnection pattern in the phytoplankton productivity in the Neva Estuary in summer (Tables 2-4).Student's t-test and Fisher's F-test showed that POLy was the main predictor in multi-regression models for both indicators of phytoplankton productivity (Tables 3, 4).With negative POLy values, the concentration of chlorophyll a and planktonic primary production in the Neva Estuary increased (Tables 3, 4).Since the mid-1990s, POLy values have been mostly in the negative phase, while positive values predominated in the second half of the twentieth century (Figure 6).Therefore, it is natural that in the modern period there is an increase in the productivity of estuarine phytoplankton (Golubkov, 2009;Golubkov et al., 2017), especially in the concentration of chlorophyll a (Figure 4A) and biomasses of most phytoplankton groups (Golubkov et al., 2021).
We were unable to find any studies of the relationship between the Polar/Eurasian pattern and phytoplankton productivity indicators.However, it has been shown that the POL pattern exhibits a statistically significant negative correlation with the degree of soil moisture and a positive correlation with air temperature in June and July in northwestern European part of Russia (Zveryaev and Arkhipkin, 2019).It is also known that POL is one of the patterns that control the "Vangengeim-Girs" circulation, which forms over the Atlantic-Eurasian region and affects the westerly winds in summer (Gao et al., 2017).Strong westerly winds limit the extension of the North Polar vortex, which usually leads to its weakening and the transition of the POL pattern to a negative phase.During this phase of the POL pattern, westerly winds, as in the case of the positive phase of NAO, may lead to an extension of the water residence time in the Neva Estuary, and to an increase in the plankton primary production.
AO teleconnection pattern did not appear to be predictors of phytoplankton productivity in multi-regression analysis (Tables 3, 4).However, Pearson's pairwise correlation analysis showed that this pattern for different periods also show statistically significant but not very high correlations with midsummer CHL (Table 2).For example, the chlorophyll a concentration in the Neva Estuary in summer correlated positively with the winter AO values.The influence of the AO teleconnection pattern on the productivity of coastal waters has not been studied very often.However, it is known that when the AO values are positive in winter, the number of cyclones arriving in the Baltic region from the North Atlantic and bringing a large amount of precipitation increases (Bednorz and Tomczyk, 2021).This, in turn, may cause an increase in phytoplankton production in the Neva Estuary (Golubkov and Golubkov, 2020).On the contrary, the negative rather than positive AO phase led to the massive development of algae in the coastal zone of the Atlantic Ocean near the Iberian Peninsula located in the southwestern part of Europe (Báez et al., 2014).There is no contradiction in these results, since it is known that the positive AO phase leads to stormy and rainy weather in northern Europe, and the negative AO phase leads to similar weather conditions in southern Europe (Thompson and Wallace, 1998).
There are even fewer studies on the effect of the EAWR teleconnection pattern on phytoplankton productivity.For example, an attempt was made to link the shift regime in the ecosystem of the Caspian Sea, which occurred in the late 1990s and mid-2000s with changes in the EAWR index, but statistically reliable results could not be obtained (Kashkooli et al., 2017).In our study, we also could not have obtained statistically significant correlations between phytoplankton productivity indicators in the Neva Estuary and EAWR teleconnection pattern (Table 2).
Our study showed that the current increase in phytoplankton productivity in the Neva Estuary is higher than it could have been in the past, provided that the regression relationships found were relevant in the past (Figure 4).This is especially true for the concentration of chlorophyll a.In the case of primary production, although the same high values could be observed earlier, but unlike in the modern period, no longer than 1 or 2 years (Figure 4B).Earlier studies have shown that ecosystem shifts occur periodically in different parts of the Baltic Sea.For example, Lindegren et al. (2012), showed a regime shift for the Kattegat straits that occurred at the turn of the 1980s and 1990s.In the mid-1980s, high primary production was observed in this part of the Baltic Sea.Then, significant changes occurred at the turn of the 1980s-1990s, and from 1992 to 2008, there was a period of low phytoplankton productivity.According to our regression models, in the mid-1980s, high plankton productivity should also have been observed in the Neva Estuary, and the 1990s and 2000s corresponded to the period of its low productivity (Figure 4).Accordingly, although these areas of the Baltic Sea are located far from each other, it is possible that the same factors could have been the trigger mechanisms for a decrease in phytoplankton productivity at that time, and such factors could be associated with atmospheric circulation.Such large-scale changes in marine ecosystems caused by atmospheric circulation have been described for the Baltic Sea (Möllmann et al., 2009) and for other marine areas (Möllmann and Diekmann, 2012;Blenckner and Niiranen, 2013).
Our analysis and the obtained multi-regression models can further make it possible to assess to what extent the changes in the ecosystem of the Neva Estuary are consistent with largescale changes in other aquatic ecosystems that have occurred in the past.At the same time, we agree with colleagues who urge caution in such reconstructions and refine the resulting models as new evidence and data sets emerge (Murphy et al., 2020).In addition, it should be taken into account that different climatic patterns can have the most significant impact on ecosystems in certain regions.For example, according to our studies, phytoplankton productivity in the Neva Estuary located in the most northeastern part of the Baltic Sea, showed the more significant relationships with the Polar/Eurasian teleconnection pattern, which determines weather conditions in more northeastern regions of Eurasia, as compared to NAO teleconnection pattern, which mainly determines weather conditions in western Europe.Some other climatic indices also showed significant correlations with concentration of chlorophyll a and plankton primary production in the Neva Estuary.Therefore, it is important to continue research in this direction, finding out which teleconnection patterns are most informative for a particular region.
It is also worth noting that in this work we did not consider the role of the anthropogenic factor, which has a strong effect on the ecosystem of the Neva Estuary (Golubkov et al., 2019).However, unlike weather fluctuations, this factor was rather high and constant from the beginning of twentieth century and changed relatively slowly over time (Golubkov et al., 2019), making up a certain background against which changes caused by climate fluctuations took place.The study also shows that, despite all human activities transforming the biosphere, global planetary processes still have a significant impact on the environment, even in ecosystems like the Neva Estuary, which are heavily affected by anthropogenic influence.This impact should be taken into account when developing measures to protect aquatic ecosystems from eutrophication, since climate fluctuations can enhance or weaken the consequences of anthropogenic influence.

FIGURE 1 |
FIGURE 1 | The upper and middle parts of the Neva Estuary with indication of sampling stations (1-17).Black lines: isobaths of 5, 10, and 20 m.Areas with dots indicate dense reeds.C1, C2-gates for vessels; D1-D6-waters gates in the St. Petersburg Flood Protection Facility.Red rectangles-the location of the Neva Estuary.Two-letter country codes are given according to ISO 3166-1 alpha-2 (International Organization for Standardization, 2020).

FIGURE 2 |
FIGURE 2 | Relationship between midsummer chlorophyll a concentrations and plankton primary production in the Neva Estuary.

FIGURE 3 |
FIGURE 3 | Plots of residuals vs. predicted values for chlorophyll a (A) and plankton primary production (B) models.

FIGURE 4 |
FIGURE 4 | Concentration of chlorophyll a (A) and plankton primary production (B) in the Neva Estuary in 1951-2020.Data for 1951-2002 calculated according to the obtained multi-regression dependences, data for 2003-2020 are obtained from field observations of the authors.

FIGURE 5 |
FIGURE 5 | Winter averaged values of the North Atlantic Oscillation (A) and Scandinavia (B) indices from 1951 to 2020 (calculated according to NOAA data, 2021).

FIGURE 6 |
FIGURE 6 | Averaged annual values of the Polar/Eurasian Index from 1951 to 2020 (calculated according to NOAA data, 2021).

TABLE 1 |
Averaged environmental variables for the upper and middle parts of the Neva Estuary during the study period.
ns, not significant.Teleconnections patterns abbreviations are given in Materials and Methods, see section "Northern Hemisphere Teleconnection Patterns."y, annual averaged, w, for winter months averaged, j, July values.

TABLE 3 |
Results of stepwise multiple regression analysis between the concentration of chlorophyll a (CHL, variable Y, predictand) and values of climatic indices (variables X, predictors).
POLy is the Polar/Eurasian Index for the year; NAOw and SCANDw are the North Atlantic Oscillation and Scandinavian indices for the winter months.

TABLE 4 |
Results of stepwise multiple regression analysis between plankton primary production (PP) (variable Y, predictand) and values of climatic indices (variables X, predictors).POLy, Polar/Eurasian Index for the year; NAOj, North Atlantic Oscillation Index for July; SCANDj, Scandinavia Index for July.