Comparison of biological and ecological long-term trends related to northern hemisphere climate in different marine ecosystems

Data from five sites of the International Long Term Ecological Research (ILTER) network in the NorthEastern Pacific, Western Arctic Ocean, Northern Baltic Sea, South-Eastern North Sea and in the Western Mediterranean Sea were analyzed by dynamic factor analysis (DFA) to trace common multi-year trends in abundance and composition of phytoplankton, benthic fauna and temperate reef fish. Multiannual trends were related to climate and environmental variables to study interactions. Two common trends in biological responses were detected, with temperature and climate indices as explanatory variables in four of the five LTER sites considered. Only one trend was observed at the fifth site, the Northern Baltic Sea, where Nature Conservation 34: 311–341 (2019) doi: 10.3897/natureconservation.34.30209 http://natureconservation.pensoft.net Copyright Ingrid Kröncke et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Launched to accelerate biodiversity conservation A peer-reviewed open-access journal


Introduction
Long-term ecological research aiming at understanding sources of natural variability in species composition, dominant species and functional diversity of marine communities is essential for disentangling the effects of natural environmental drivers and climate variation from the direct effects of local human activities (Giron-Nava et al. 2017;Mirtl et al. 2018).Many marine ecosystems around the world have undergone dramatic changes in species composition or a reduction in species diversity (Perry et al. 2005;Worm et al. 2006;Steffen et al. 2007;Suikkanen et al. 2007Suikkanen et al. , 2013;;Rockström et al. 2009;Kröncke et al. 2011) due to human activities and climate-driven ocean warming in the northern Pacific and Atlantic Oceans during recent decades (Deser et al. 2000;IPCC 2014).
An analysis of changes across different biological components and geographic regions in relation to climate indices has the potential to provide a novel insight into the role of climate in influencing trends in marine biodiversity (Dornelas et al. 2014;Bowler et al. 2017a, b).Changes and shifts in biodiversity of single or different biological components can be part of biological regime shifts (BRS).BRS are defined mainly on the basis of changes in the ecosystem as a whole (Scheffer et al. 2001;Collie et al. 2004;Conversi et al. 2015) and are typically characterized by infrequent and abrupt changes in ecosystem structure and function, occurring at multiple trophic levels and on large geographic scales (Collie et al. 2004;Cury and Shannon 2004;deYoung et al. 2008;Bakun 2005;Lees et al. 2006).Marine biological regime shifts have been explained mainly as a result of multiple causes, e.g.climate regime shifts (CRS), overexploitation of resources or a combination of climate and volcanic eruptions (Scheffer and Carpenter 2003;Fisher et al. 2015;Gardmark et al. 2015;Levin and Möllmann 2015;Reid et al. 2016).Of particular note was the 1988 BRS, which was detected in marine ecosystems globally (e.g., Reid et al. 2016;Conversi et al. 2010), while the 2000 CRS/BRS has been subject to only a few long-term studies in specific marine systems (Kröncke et al. 2013a;Beaugrand et al. 2014;Dippner and Kröncke 2015;Perretti et al. 2017).All these studies have shown that signals of climate variability can be detected by an increase in warm-temperate species, a decrease in cold-temperate species at various trophic levels of the marine ecosystem, and subsequently in changes in species composition and species richness.The trends in biodiversity described above were strongly correlated with increasing sea surface temperatures (SST) and regional climate variability (Dornelas et al. 2014).These variables have been proven useful for characterizing basin-scale variation in the world's oceans and offer great potential for investigating inter-annual and inter-decadal climate oscillations, biological and ecological regime shifts, and climate teleconnections.Various climatic indices have been developed for these purposes and their modes typically involve different mechanisms over different regions.For example, the North Atlantic Oscillation (NAO) (Hurrell 1995) and the North Pacific Index (NPI) (Trenberth and Hurrell 1994) are the leading climate modes of sea level pressure variability over the northern Atlantic and Pacific Oceans, respectively.These indices along with the Pacific Decadal Oscillation (PDO) (Mantua 1999) and the El Niño/Southern Oscillation (ENSO) (Philander 1990) are the four dominant modes of climate variability in the Northern Hemisphere (Tsonis et al. 2007).NAO has been widely used to characterize the role of climate variability in the responses of terrestrial (Mysterud et al. 2003), freshwater (Straile et al. 2003) and marine ecosystems (Drinkwater et al. 2003;Dippner 2006).Hameed and Piontkovski (2004) developed an indicator alternative approach to the NAO that uses the pressures at the atmospheric centers of action (Azores High and Icelandic Low, IL) as independent variables and showed that the standardized time series of the pressure, zonal and meridional position of the IL is a good predictor of various ecosystem attributes such as the north wall position of the Gulf Stream and the zooplankton abundance in the Gulf of Maine.In the southern North Sea, the NAO was a good predictor of the structure of macrofauna communities in the following spring (Kröncke et al. 1998), which allowed the development of forecast equations (Dippner and Kröncke 2003).However, the NAO lost its persistency around the year 2000, when a climate regime shift (CRS) occurred (Dippner et al. 2014).Thus, the North Sea Environmental index (NSE), a multivariate predictor constructed from a mixture of global and regional climatic time series, was developed to better predict inter-annual and inter-decadal climate variability in the North Sea (Dippner and Kröncke 2015).Similar climate predictors have been developed for other regions, such as the Arctic Oscillation (AO) (Thompson and Wallace 1998), the Pacific North American pattern (PNA) (Barnston and Livezey 1987), the Indian Ocean Dipole (IOD) (Saji et al. 1999, Webster et al. 1999), and the Western Mediterranean Oscillation (WeMO) (Martin-Vide and Lopez-Bustins 2006).
Combining different climate modes with vertical propagating Rossby waves to and from the stratosphere (Ineson and Scaife 2009) results in a coupled "super-loop" NAO → PDO → ENSO → PNA → stratosphere → NAO, which contains the major signals of low-frequency variability in Northern Hemisphere (Wang et al. 2009).The analyses of climate modes are used for the identification of climate phenomena such as teleconnections and regime shifts.For example, high spectral phase coherence between two climate modes indicates the existence of a teleconnection.Such teleconnections have been identified between the NAO and ENSO (Huang et al. 1998), ENSO and IOD (Saji and Yamagata 2003;Tozuka et al. 2007), and ENSO and AO (Jevrejeva et al. 2003).
CRS are often considered as changes in the trend of global mean air temperature.Yasunaka and Hanawa (2002) used empirical orthogonal functions' analysis to evaluate the SST and sea level pressure (SLP) fields in the northern hemisphere and identified different CRS in the 20 th century.Swanson and Tsonis (2009) proposed the following mechanism for CRS.When different modes of climate variability are synchronized (defined as the root-mean-square correlation coefficient between all pair of climate modes) and the coupling strength (defined as the phase of an individual climate mode relative to the phases of other modes) between the modes increase simultaneously, the climate system becomes unstable and moves into a new state (Tsonis et al. 2007;Swanson and Tsonis 2009;Wang et al. 2009).Such synchronization of CRS occurred four times during the 20 th century (1910-1920, 1938-1945, 1956-1960, and 1976-1981) and an increase in coupling strength occurred three times (1910-1920, 1938-1945, and 1976-1981).The consequence was a change in the trend of global mean air temperature (Swanson and Tsonis 2009).The last identified synchronization and increase in coupling strength causing a climate shift was identified in 2001/2002 (Swanson and Tsonis 2009).These identified periods were related to the turning points of the Atlantic Multidecadal Oscillation (AMO; Sutton and Hodson 2005).The AMO has a periodicity of ~60-80 years (Schlesinger and Ramankutty 1994).Thus, a relatively cold period at the beginning of the Twentieth Century was followed by a warm period in the 1940s and 1950s, another cold period in the 1970s and 1980s, and a warm period in the 1990s (Ting et al. 2013).
Comparative ecosystem analysis has been effectively used to improve our understanding of the processes controlling the biodiversity, productivity, and resilience of marine ecosystems (Murawski et al. 2009), and develops general principles of biological shifts (Conversi et al. 2015;Beaugrand et al. 2015).In this paper, we analyzed longterm trends in biological and climate data from five sites belonging to the European (LTER-Europe) and International (ILTER) Long Term Ecological Research networks (www.ilter.network),located in different marine basins along a latitudinal gradient in the northern hemisphere (i.e., the Northeastern Pacific, Western Arctic Ocean, Northern Baltic Sea, Southern North Sea and Western Mediterranean Sea; Fig. 1).The main objectives of the study were to investigate: (1) whether common trends in abundance or biomass across different biological groups or ecological parameters of pelagic and benthic components of marine communities occurred in the last two decades, and (2) whether biological shifts corresponded to shifts in ocean or air temperature and regional climate variability in the different marine regions.We also sought (3) to determine, whether any such changes could be explained by species-specific temperature affinities during cold and warm phases.

Material and methods
The International Long-term Ecological Research Network (ILTER), founded in 1993, comprises 44 active member LTER networks representing 700 LTER Sites and ~80 LTSER Platforms across all continents, active in the fields of ecosystem, criti-cal zone and socio-ecological research (Mirtl et al. 2018).ILTER maintains global environmental observation and ecological research networks and infrastructures for excellent science.The role of ILTER is to investigate ecosystem structure, function, and services in response to a wide range of environmental strengths using long-term, place-based research.

Long-term ecological research sites
The Santa Barbara Coastal LTER (https://deims.org/dbd399ed-9c26-4621-b479-7ab505c8aa35) in southern California (USA, Fig. 1) was established in 2000 to investigate the roles of a wide range of marine and terrestrial processes in structuring kelp forest ecosystems in the face of changing climate and human use.The study site extends 120 km and is bounded by the continental coastline and an archipelago of four islands located 30 -40 km offshore.Its east-west orientation contributes to complex circulation patterns that are influenced by the intersection of the cold southward-flowing California Current and the warm, northward-flowing southern California Counter Current (Hendershott and Winant 1996;Harms and Winant 1998).This circulation and resulting strong temperature gradients create a major biogeographic boundary at the western border of the Channel for a wide variety of marine taxa (Valentine 1966;Briggs 1974), and make the Santa Barbara Channel an ideal location to examine shifts in species composition driven by climate variability.Details in sampling and methods are given by Reed et al. (2016).
The arctic LTER Observatory HAUSGARTEN (https://deims.org/f6d9ed12-6bc1-47fb-8e81-ef24e9579596) in the Fram Strait between NE Greenland and the Svalbard archipelago (Fig. 1) was established in 1999 to detect and track the impact of large-scale environmental changes on the marine ecosystem in the transition zone between the northern North Atlantic and the central Arctic Ocean (Soltwedel et al. 2005(Soltwedel et al. , 2016)).Today, HAUSGARTEN observatory constitutes a network of 21 permanent sampling sites, the majority of which are located along a bathymetric transect between ~250 and ~5,500 m water depth at about 79N from the Kongsfjorden (Svalbard) in the east, along the Vestnesa Ridge towards the Molloy Hole (i.e., one of the deepest known depressions in the Arctic Ocean), and across the Greenland continental margin.Three sampling sites close to the ice edge between 79°30'N and 80°00'N in the north-eastern Fram Strait and a supplementary site in a permanently ice-free area at 78°30'N in the eastern part of the strait complete the network.The Fram Strait is the only deep-water connection between the Nordic Seas and the central Arctic Ocean, with a sill depth of about 2600 m.The hydrography in the eastern part of the strait is characterized by the inflow of relatively warm and nutrient-rich Atlantic Water into the central Arctic Ocean (Beszczynska-Möller et al. 2012).Cooler and less-saline Polar Water exits the central Arctic Ocean as the Eastern Greenland Current in the western part of the Fram Strait (de Steur et al. 2009).Hydrographic patterns in the strait result in a variable sea-ice cover, with predominantly ice-covered areas in the west, permanently ice-free areas in the southeast, and seasonally varying ice conditions in the central and north-eastern parts.Details in sampling and methods are given by Soltwedel et al. (2005Soltwedel et al. ( , 2016)).
The Baltic Sea LTER site Seili (https://deims.org/9d4222a2-c50f-4fac-8b1d-3b685072b34d) is located in the Archipelago Sea in Finland (60°15'33"N; 21°57'39"E) (Fig. 1), which is part of the northern Baltic Sea, and hence a brackish water marine area.According to the Water Framework Directive of the European Union, this coastal area belongs to the Southwestern middle archipelago of the Finnish coastal waters.The area is characterized by a large number of skerries and islands; the median depth of this part of the Archipelago Sea is 40 m.The depth of the sampling site is 50.6 m.The salinity ranges between 5.1 and 7.0 (median value 5.9).The sampling area is described in detail by Kauppila (2007 and references therein).Samples for phytoplankton have been taken since 1991.Details in sampling and methods are given by HELCOM (2017).
The LTER North Sea Benthos site (https://deims.org/50946250-c0fa-41b0-a917-17d2a3992eee)constitutes a network of seven sampling areas in the North Sea.Box A, the area used for this study, is situated in the German Bight about 25 nautical miles northwest of the Island of Helgoland, in close proximity to the 30 m depth contour near the former glacial valley of the River Elbe (Fig. 1).The mean depth of this area is 40 m and the water column is generally well mixed due to tidal mixing throughout the year.The freshwater run-off from continental rivers results in low salinities and regionally high nutrient input in the area, which can cause enhanced primary production and high food supply for benthic fauna (Kröncke et al. 2004).The area is also characterized by strong seasonal variability in temperature and salinity as well as by high bottom stress due to strong tidal currents (Neumann et al. 2017).The sediments consist of more than 20% (<63 μm fraction) mud in the south-west corner gradually decreasing towards the north-east corner (0-5%).The time series started in 1998 and has since been sampled annually with the exception of 2013, when it was not sampled.Details in sampling and methods are given by Neumann et al. (2008).
The Mediterranean Sea LTER site Gulf of Olbia (IT14-002-M) (https://deims.org/3178d0fb-0789-4992-9c51-1ddb50b7e871) has belonged to the LTER-Italy network since 2006.It is situated on the eastern coast of Sardinia (Italy, Fig. 1).Morphologically, it is a ria, 5 km long, 2 km wide, with a total area of 6.5 km 2 , a mean depth of 5 m, and a maximum depth of 10 m.The Gulf receives an inflow from the Padrongianus River in its southern part.The town of Olbia, located in the inner part of the Gulf, represents one of the most important passenger harbors in the Mediterranean in addition to a commercial and an industrial harbor (www.olbiagolfoaranci.it).It is also the largest mussel-and clam-farming area in Sardinia (Sannio et al. 1996;Bazzoni et al. 2015) and one of the most important in Italy.Since the early 1990s, the Gulf has been monitored at three stations for phytoplankton abundance and composition as well as environmental variables.In particular, the Gulf of Olbia has been monitored with respect to its trophic state, phytoplankton abundance and composition, focusing on the presence of harmful algae, since the early 1990s (Sannio et al. 1996(Sannio et al. , 1997;;Lugliè et al. 2003a, b;Garcés et al. 2007;Satta et al. 2010;Bazzoni et al. 2015;Pulina et al. 2016).The first data were collected in 1987 in connection with an event of water discoloration due to a phytoplankton bloom (Sechi et al. 1987).
The collection of LTER data started in 1992 for phytoplankton.Details in sampling and methods are given by Pulina et al. (2016).

Sampling, data sets and data analysis
Table 1 provides an overview on the LTER study sites, including the sampling period, and the climate and biological data sets considered in the dynamic factor analysis (DFA).

Dynamic Factor Analysis (DFA)
DFA is a smoothing and dimension reduction technique to identify common trends in multivariate time-series and to determine the effect of explanatory variables.The "common trends" represent the underlying dynamic pattern over time in the considered system (Zuur et al. 2003a).Unlike other dimension-reduction techniques such as principal component analysis (PCA) or multidimensional scaling (MDS), DFA is especially designed for sequential time-series data.Another advantage of this method is that it can deal with missing data in time series, as they also occur in our time series (e.g.Olbia, Hausgarten, North Sea).
The time series were modeled as a function of a linear combination of common trends, a constant level parameter, two or more explanatory variables, and noise (Zuur et al. 2003a, b).A detailed mathematical description of the DFA algorithm is given in Zuur et al. (2003a) and Zuur et al. (2007).

Abundance of 16 epifauna species
Phytoplankton cell density as percentage composition on the base of contribution of 9 major groups DFA was used to analyze time series of phytoplankton (Baltic Sea, Western Mediterranean, Arctic Ocean), benthic fauna (North Sea, Arctic Ocean), benthic bacteria (Arctic Ocean) and temperate reef fish (Santa Barbara Channel).Water or air temperature and various climate indices (PDO, IL, NAO, NSE, WeMO) were used as explanatory variables.Several DFA models were tested for each time-series, ranging from the simplest (one common trend plus noise) to the most complex (two common trends, up to three explanatory variables plus noise).Models were fitted with both a diagonal covariance matrix and a symmetric non-diagonal covariance matrix and compared using the Akaike's information criterion (AIC) (Akaike 1974), which depends on the maximum likelihood estimates and is a measure of goodness-of-fit.Depending on the number of parameters, the model with the lowest AIC value was selected as the most appropriate one.Factor loadings were used to assess the relationship between response variables and the common trends (cut-off point of 0.2).The diagonal elements of er-ror covariance matrix obtained by the model were used as a measure of "misfits" of the fitted model.Large diagonal elements mean that the corresponding time series are not fitted well (Zuur and Pierce 2004).Significant (p < 0.05) relationships between the response and explanatory variables were assessed by using the estimated regression coefficients and their associated t-value, whereby t-values > 2 (in absolute sense) indicate strong relationships.All time series data (response variables) were square root transformed and both response and explanatory variables were standardized by subtracting their time-averaged mean and dividing by their temporal standard derivation (normalization).Computations were carried out using Brodgar computer software, developed by Highland Statistics LTD., UK (www.brodgar.com).

Results
The transformed and standardized time-series of the different biological components of marine communities and environmental parameters are presented in Fig. 2. The time series showed considerable variability with no obvious trends.To estimate the underlying common trends of the various biological components, different sets of DFA models, depending on the length of the time-series in relation to different climatic indices, were used for each of the LTER study sites (Fig. 3).The AIC values (Table 2) indicate the best model fit in each study case.
LTER Santa Barbara Channel: Sixteen dynamic factor models were calculated to estimate the underlying common trend of the fish time-series differing in the covariance matrix employed (diagonal versus non-diagonal), the number of trends (one or two) and the included explanatory variables (mean summer bottom temperature, PDO index or both).The AIC values indicated that the best model fit was obtained for a non-diagonal matrix with two common trends explained by bottom temperature and the PDO index (AIC = 415) (Table 2).
Both trends showed overall increases from 2001-2017 (Fig. 3).The first common trend was characterized by a sharp increase from 2001 to 2008 and a decrease until 2012 followed by a leveling out and slight increase.The second common trend showed a modest decrease from 2001 to 2006, followed by an increase afterwards (Fig. 3).Factor loadings showed that the first trend was related to the blue rockfish Sebastes mystinus (0.30).The second trend was correlated with the señorita Oxyjulis californica (0.41), the blacksmith Chromis punctipinnis (0.32), the pile surfperch Phanerodon furcatus (0.2) and the blackeye goby Rhinogobiops nicholsii (-0.23).Relatively low diagonal elements of the error covariance matrix (< 0.50) were obtained for all of these species except the pile surfperch (0.58) highlighting a good model fit of the individual series for these species (Table 3).
The estimated t-values for the explanatory variables bottom temperature and PDO are given in Table 4. Bottom temperature had a significant influence on the painted greenling Oxylebius pictus (-3.36), the señorita Oxyjulis californica (-2.48), and the kelp bass Paralabrax claturatus (2.44).Furthermore, the abundances of the blacksmith   Chromis punctipinnis (6.14), black surfperch Embiotica jacksoni (-2.90) and the kelp perch Brachyistius frenatus (-2.27) were significantly related to the PDO.LTER Observatory HAUSGARTEN: Thirty-two dynamic factor models, using SST as well as the intensity and positions of the IL as explanatory variables, revealed two trends in the planktonic community and underlying sediments.The AIC values (Table 2) indicated that the best model fit was obtained for a diagonal matrix with two common trends and bottom temperature, air pressure and latitude as explanatory variables (AIC = 424).The first common trend showed two major peaks in 2004/2005 and 2011/2012.The second common trend showed an overall decrease between the years 2000 and 2016 (Fig. 3).Factor loadings showed that the first trend was related to chlorophyll a (0.32), nanoflagellates (0.75), Dinophyceae (0.79), Phaeocystis pouchetii (-0.35),Coccolithophyceae (0.50), organic matter in the sediments (AFDW; 0.24), particulate proteins in the sediments (PROT; 0.57) and meiofauna (-0.98).The second trend was related with Bacillariophyceae (0.25) and the sediment-bound chloroplastic pigments (CPE; -0.31).Among the considered biological groups, the best model fits were found for the proportion of Bacillariophyceae, nanoflagellates, Coccolithophyceae and Phaeocystis pouchetii (Table 3).
t-values indicated a significant influence of SST on the proportion of Bacillariophyceae (-4.64) and Phaeocystis pouchetii (3.65) (Table 4).The air pressure anomaly of the IL strongly covaried with the proportion of Bacillariophyceae (4.08), Coccolithophyceae (-3.51), and protein concentrations (4.24) in the sediments (Table 4).The proportions of Coccolithophyceae (4.79) and Bacillariophyceae (-4.47) also showed a strong relationship to the latitudinal position of the IL (Table 4).Meiofauna density at HAUSGARTEN observatory was the only faunal parameter included in the DFA.Unfortunately, meiofauna data had to be restricted to results provided by Hoste et al. (2007) and Grzelak (2015) for the central HAUSGARTEN site (2500 m water depth), covering the years 2000 till 2009.Results from the DFA revealed no strong relationships with the explanatory variables chosen for the analysis; however, the model exhibited a fairly good fit with the observed values (Table 4).
LTER Baltic Sea Seili: Sixteen dynamic factor models were calculated to estimate the underlying common trend of the phytoplankton time-series differing in the covariance matrix employed (diagonal vs. non-diagonal), the number of trends (one or two) and the included explanatory variables, i.e.SST and NAO (Table 2).The AIC values indicated that the best model fit was obtained for a diagonal matrix with one common trend and no explanatory variables (AIC = 549).The common trend showed an increase from 1991 to 2005, and leveling out afterwards (Fig. 3).Based on the factor loadings, the trend was positively related to biomasses of Chlorophyta (0.32) and Cyanophyceae (0.29) and negatively to Cryptophyceae (-0.20).Observed and fitted series for the best correlated taxa indicated that the first two groups were fitted reasonably well, with measures of 0.41 and 0.51, respectively (Table 3).Because the best model fit did not contain any explanatory variables, relations of phytoplankton time series to explanatory variables could not be examined.LTER North Sea Benthos Observatory: Sixteen dynamic factor models were calculated to estimate the underlying common trend of the epibenthic time-series differing in the used covariance matrix (diagonal and non-diagonal), and the number of trends (one or two) and the included explanatory variable (SST or NSE or both).The AIC values indicated that the best model fit was obtained for a nondiagonal matrix with two common trends and SST and NSE as explanatory variables (AIC = 620) (Table 2).
The first common trend showed a sharp decrease from 1998 to 2007, followed by an increase afterwards (Fig. 3).The second common trend was characterized by an increase from 1998 to 2005 and a steady decrease from 2005 onwards.Factor loadings showed that the first trend was positively related to the brittle star Ophiura albida (0.27) and the brown shrimp Crangon crangon (0.20), and negatively related to the Auger shell Turritella communis (-0.20).The second trend was positively correlated with the shrimp Philocheras bispinosus bispinosus (0.27) and the swimming crab Liocarcinus holsatus (0.27) and negatively with the brittle star Ophiura ophiura (-0.20).Observed and fitted series for the best correlated species (Table 3) indicated a good model fit for the above species.In contrast, the individual model fit was low for the sea star Astropecten irregularis (0.82), the masked crab Corystes cassivelaunus (0.78) and the goby Pomatoschistus minutus (0.77).
Mediterranean Sea LTER Gulf of Olbia: The best fit of the sixteen model outputs for phytoplankton was one with a non-diagonal error covariance matrix, containing two common trends (Fig. 3), and with air temperature and WeMO as explanatory variables (the lowest AIC = 308) (Table 2).
The first common trend showed higher and similar values at the beginning (1996)(1997) and at the end of the considered time series (2013)(2014)(2015), with a central part of lower values (1998-2012, minimum value in 2005) (Fig. 3).The second common trend was slightly increasing at the beginning (1995)(1996)(1997)(1998)(1999)(2000)(2001), with a subsequent prominent decrease till the end of the study period, except for a small rise in 2011.
According to the largest factor loadings (Fig. 3), percentage contribution of Bacillariophyceae to total phytoplankton cell density followed strongly the first common trend (0.26) in a positive way.On the contrary, percentage contribution of nanoplankton to total phytoplankton cell density resulted in a negative factor loading (-0.27).The second common trend strongly correlated with the percentage contribution of Cryptophyceae to total phytoplankton cell density, with a negative sign (-0.25).Best fits were found for the time-series of nanoplankton (0.25), Bacillariophyceae (0.34), Chrysophyceae (0.41) and Cryptophyceae (0.49) (Table 3).
Considering the estimated t-values (Table 4), air temperature had a significant influence on percentage contribution to total phytoplankton cell density of Chrysophyceae (t value = -3.45),whereas WeMO influenced the percentage contribution of Chlorophyceae (t value = -2.87)and Dinophyceae (2.38).

Trends in ocean temperature and climate modes
Data from the five LTER sites considered in this study revealed an overall increase in temperature (SST, bottom and air temperature; data not shown).SST increased 0.8-3 °C at the northeastern Pacific, North Sea and Baltic Sea sites between 1995 and 2016.Air temperature at the LTER site in the Western Mediterranean Sea increased by 0.8-1.1 °C, while at the Arctic LTER site HAUSGARTEN SST increased by 0.06 °C y -1 for the time-period 1997-2010 (Beszczynska-Möller et al. 2012), but also revealed distinct interannual variability (Walczowski et al. 2017).The time-series of summer observations exhibited 5-6 year cycles with distinct warm-water anomalies (WWA).Recent WWAs in Fram Strait occurred between 2005 and 2007, in the years 2011/2012, and towards the year 2016.
Figure 4 shows the long-term variability of climate indices.Temperature increases in the HAUSGARTEN, in the North Sea and the Baltic Sea were related to increasing IL/NAO/NSE indices.The increase in air temperature at the LTER site in the Western Mediterranean Sea coincided with extreme negative WeMO values.There was a strong shift at the Santa Barbara Channel LTER site in the north-eastern Pacific from decreasing bottom temperature and negative PDO indices towards increasing SST of about 3 °C and positive PDO indices around 2014.The 2000 CRS is obvious as well, a similar pattern in all indices around 2010, which might be a hint of another CRS, as yet undetected.

Discussion
The comparison of climatic and long-term ecological data at the five LTER sites located in the north-eastern Pacific, the western Arctic Ocean, the northern Baltic Sea, the southern North Sea and the western Mediterranean Sea revealed a great similarity in common trends in marine systems in the northern hemisphere.

Common multiannual trends and biological shifts
At four of the five LTER sites considered in this study, the DFA revealed the presence of two common trends with temperature and climate indices as explanatory variables.Despite the different biological components and marine ecosystems analyzed, multiannual common trends within a site were quasi-synchronous, but usually in opposite directions.Common trends first crossed near the end of 1990s/early 2000s and then again around 2010, which coincided with the 2000 CRS and probably a new, yet to be described CRS around 2010.We found one inversion of tendencies of biological components at the Arctic LTER HAUSGARTEN and two inversions at sites in the north-eastern Pacific Ocean (Santa Barbara), the Western Mediterranean Sea (Gulf of Olbia) and the North Sea.A similar signal was also observed at the Baltic Sea site (Seili), where the single common trend reached its maximum in 2005, for which no explanatory variable was identified, similar to other long-term plankton studies (e.g. in the North Sea, van Beusekom et al. 2009).Three LTER data sets (Western Mediterranean Sea, North Sea, Baltic Sea) analyzed in this study started before 2000.The observed changes in the common trends and composition of the marine communities at these sites coincided with the 2000 CRS (Tsonis et al. 2007;Swanson and Tsonis 2009;Wang et al. 2009) and could provide further evidence for the presence of BRS, as already reported for other marine systems worldwide (Bond et al. 2003;Peterson and Schwing 2003;Collos et al. 2009;Pulina et al. 2011;Kröncke et al. 2013a;Beaugrand et al. 2014;Dippner and Kröncke 2015;Reid et al. 2016;Perretti et al. 2017).Beaugrand et al. (2015) reported quasi-synchronous regime shifts in the late 1980s and the mid-to late 1970s related to temperature and the Arctic circulation in multiple marine systems from two oceans and three regional seas in the Northern Hemisphere.Because changes monitored in the north-eastern Pacific Ocean since 2002 are similar to those from three other LTER sites considered in our study, they can probably also be related to the 2000 CRS.Turning points of the common trends at four of the sites in the mid to late 2000s might be explained as biological shifts between two CRS.Four of the five LTER sites revealed also biological shifts around 2010, which coincided with shifts in the climate indices and might be evidence for a new CRS.
The congruent trends in climate indices and CRS at the five LTER sites suggest that decadal fluctuations in atmospheric and ocean circulation are teleconnected between the Atlantic and Pacific Ocean regions as also found for previously reported CRS (Schwing et al. 2003).The congruent trends of climate indices and of common trends in biological variability provide evidence that a major part of inter-annual and inter-decadal biological variability can be attributed to physical forcing.

Biological responses to CRS and the role of temperature
We found consistent biological modifications to abrupt or continuous increases in sea water and air temperature and associated climatic indices at all five of the LTER sites analyzed in this study.
Temperature modulates a multitude of processes at the cellular, organismic and ecosystem levels of organization, and it can act as a mediator between organisms and climate (Kröncke et al. 2013a;Dippner and Kröncke 2015).Recent changes observed in the northern Atlantic communities caused by the 2000 CRS were similar to those caused by the CRSs in the 1920s and 1930s (Drinkwater 2006;Deser 2000).Ecosystem changes associated with warmer than normal SST and positive values of the NAO in the 1920s and 1930s included a general northward movement of cold-temperate species of fish and benthic invertebrates.These changes reversed during the CRS of 1938-1945(Tsonis et al. 2007;Swanson and Tsonis 2009;Wang et al. 2009), when temperature dropped again.Additional support for the important role of temperature in structuring marine communities comes from other studies that examined the relationship between population trends and species latitudinal distributions, which may, to some extent, act as a proxy for temperature niche of fish (Holbrook et al. 1997;Brooks et al. 2002;Heath 2005;Reed et al. 2016) and marine invertebrates (Neumann and Kröncke 2011;Birchenough et al. 2015;Hiddink et al. 2015;Gardmark et al. 2015: Bowler et al. 2017a, b).
In our study, at the south-eastern North Sea LTER site, the first occurrence of the epifaunal angular crab Goneplax rhomboides coincided with the 2000 CSR.The angular crab extended its distribution range from the north-eastern Atlantic to the North Sea, which was facilitated by an increase in water temperature (Neumann et al. 2013).The link between temperature and the occurrence of the angular crab is most likely given by effects on the survival of its larvae as it was the case for the Pacific Oyster Crassostrea gigas, when rising temperatures facilitated the successful survival of larvae and promoted the colonization of the entire German Wadden Sea (Brandt et al. 2008).Further, a shift in the functional composition of epifauna occurred in 2002, largely coinciding with the 2000 CSR (Neumann and Kröncke 2011).These results suggest that climate-induced variability of SST primarily affects the reproduction of epifaunal species rather than other functional traits.Drivers of this functional shift were the shrimps Crangon and C. allmanni as well as the hermit crab Pagurus bernhardus, whose abundances were strongly related to SST and NSE in this study.This corresponds to the results of Henderson et al. (2006), who also found that the recruitment of C. crangon is correlated to SST and the winter NAO index.This could be explained by the positive link between temperature and the duration of the larval stage as well as the frequency of breeding of shrimp species (Bergström 2000;Wear 1974).In contrast, Temming and Damm (2002) found higher recruitment of C. crangon after cold winters and postulated that predators get over-saturated in years with a pronounced peak of recruits and consequently more recruits would survive.Direct effects of temperature on key stages of reproduction are also known for P. bernhardus where timing and frequency of breeding is favored by cold water temperatures (Lancaster 1990).Unfortunately, the study by Neumann and Kröncke (2011) ended in 2008, but if the results of our analyses are generalizable, then we would expect that the 2010 CRS had similar effects on the functional composition of epifauna in the southern North Sea since these three species were positively related to the DFA trends.The epibenthic community of the North Sea LTER site was also severely affected by the cold winter of 1995/1996, resulting in the outbreak of the opportunistic brittlestar Ophiura albida.This outbreak was followed by characteristic post-disturbance succession stages from 1998 to 2000 (Neumann et al. 2008(Neumann et al. , 2009)), which is clearly reflected in trend 1 of the DFA analysis.The period from 2000 onwards coincided with the 2000 CRS and was characterized by a continuous decrease of O. albida and by an increase in diversity and secondary production due to an increase in SST.
The biomass and diversity of reef fish at the Santa Barbara Coastal LTER site also shifted towards warm-temperate species with warm water affinities (Reed et al. 2016).
The DFA results for fish from this site demonstrated that species composition responded strongly to temperature and climate mode, the PDO.Changes in fish communities were driven by common species that tend to have warm-water or cold-water affinities such as the blue rockfish Sebastes mystinus, the señorita Oxyjulis californica, the blacksmith Chromis punctipinnis and the blackeye goby Rhinogobiops nicholsii that, interestingly, are all partially or wholly planktivorous and depend on zooplankton including copepods and salps for their diet (Love and Ebeling 1978).This trophic dependency may make these species more responsive to climatic forcing.
It is widely recognized that temperature has also a strong influence on the physiology of planktonic organisms because it controls basic metabolic processes (Vidussi et al. 2011).The role of temperature in structuring phytoplankton is still strongly debated, in part due to the covariation of temperature, physical structure of water masses, nutrients and grazing pressure in marine ecosystems (Gardner et al. 2011;Boyce et al. 2015;Marañón et al. 2015;Sommer et al. 2017).Instead, there is growing evidence that the decline in cell size together with the shift of species' ranges toward higher latitudes and changes in phenology, could be considered as one of the major universal biotic responses to global warming, even if with no contradictory evidence (Daufresne et al. 2009;Sommer et al. 2017).Cascading consequences of these changes on marine pelagic ecosystem functioning are still to be evaluated and should be considered in the context of a revised paradigm of pelagic ecosystem structure and function (Verity et al. 2002).Previous studies determined a shift from larger to smaller species at Western Mediterranean, Baltic and Artic LTER sites around the mid-2000s (Suikkanen et al. 2013;Nöthig et al. 2015;Pulina et al. 2016).The results obtained in the present study showed that these common trends in phytoplankton composition found at the above mentioned sites could be connected to climate variability and especially to the temperature variation.For example, at the LTER Arctic HAUSGARTEN observatory, the WWA during the years 2005-2007 caused a major shift in the composition of unicellular plankton organisms with the affirmation of the coccolithophores Phaeocystis pouchetii, a significantly increased proportion of nanoflagellates between 2009 and 2011 and a reduction of diatoms (Nöthig et al. 2015).A similar change was observed at the LTER Mediterranean site, but slightly earlier during 2002-2004, in correspondence with the strongest variation of temperature.An intense decrement in summer phytoplankton abundance was observed due to an abrupt Bacillariophyceae decrease and an anomalous increase of Prasinophyceae and Chrysophyceae in 2002.There was also the beginning of a more lasting predominance of nanoplankton replacing larger species of Bacillariophyceae, Dinophyceae and Euglenophyceae (Pulina et al. 2016).The LTER Baltic site experienced a significant increase in Cyanobacteria biomass and, similarly to the Arctic site, in small Coccolithophyceae of the genus Chrysochromulina coupled to the increase in SST (Suikkanen et al. 2013).In the present study, results from the DFA confirmed these trends on the considered longer dataset, especially for the northernmost sites.Instead, at the Mediterranean site, further relevant modifications were detected in the phytoplankton composition from 2011, when further strong temperature instability occurred, bringing to a turnaround for Bacillariophyceae and to an increase for Cryptophyceae.
Quantitative and qualitative changes in phytoplankton have unavoidable consequences at the other trophic levels.Several studies of benthic communities revealed that their diversity and function depended on the amount and quality of organic matter produced in the water column, even at extreme depths (Rowe and Pariente 1992;Clare et al. 2017).The DFA's results of the Arctic HAUSGARTEN time-series (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) revealed strong relationships between the phytoplankton composition and concentrations of biogenic sediment compounds, indicating organic matter in surface sediments with the explanatory variables chosen for our analysis.Higher amounts of sedimentbound chloroplastic pigments seemed to provide additional food supply for meiofauna and megafauna, which resulted in increases in abundance with a time lag of 1-2 years in larger organism size classes, as also found in other deep-sea areas (Sibuet et al. 1989;Rowe and Pariente 1992;Kröncke et al. 2013b).

Conclusions
Time series data collected at the five widely distributed marine ILTER sites in the Northern Hemisphere revealed relatively synchronous changes in the abundance and species composition of marine biota across the study regions.The community changes coincided with the 2000 CRS and provide the first data on a possible additional CRS in 2010, highlighting the existence of teleconnections among the climate modes in the Northern Hemisphere.Although recent climate models predict further increases in global temperature, future CRS might cause unexpected changes in marine communities.The significant relationship described in this study between climate modes and marine communities mediated by sea water or air temperature in the Northern Hemisphere can aid in predicting future changes in these marine systems in response to climate variability and ocean warming.Such improvements in ecological forecasting should prove useful to environmental managers responsible for implementing measures aimed at mitigating adverse ecological effects of future changes in climate.Our findings highlight the importance of spatially distributed quantitative Long-Term Ecological Research in developing a predictive understanding of ecological responses to climate change in the world's oceans and they underpin the continued need for long-term research within the global and the national LTER networks, as already suggested in similar studies (Pugnetti et al. 2013;Mirtl et al. 2018;Morabito et al. 2018).Last but not least, ILTER facilitate collaborations among scientific groups with different scientific topics and this study, catching the opportunity of the Special Issue proposed by the LTER-Italy network, demonstrates it.

Figure 1 .
Figure 1.Locations of the LTER study sites.

Figure 2 .
Figure 2. Transformed and standardized time-series at the various LTER sites.AFDW = organic matter in the sediments as ash-free dry weight, CPE = sediment-bound chloroplastic pigments, PROT = particulate proteins in the sediments.

Figure 3 .
Figure 3. Common trends (left) and corresponding factor loadings (right) for the five LTER sites obtained by means of DFA.Only factor loadings above the cut-off of 0.2 in absolute value are shown.Common trends and factor loadings are untitled.Dashed line in graphs indicates the confidence interval of the DFA model.Bac = Bacillariophyceae, Chl = Chlorophyceae/Chlorophyta, Coc = Coccolithophyceae, Cry = Cryptophyceae, Cya = Cyanophyceae, Din = Dinophyceae, nano = nanoplankton/nanoflagellates, Ppou = Phaecystis pouchetii, CHLa = chlorophyll a, AFDW = organic matter in the sediments as ash-free dry weight, CPE = sediment-bound chloroplastic pigments, PROT = particulate proteins in the sediments, Meio = meiofauna.

Figure 4 .
Figure 4. Multiannual yearly mean variability of the different climate indices (PDO, IL, NAO, NSE, WeMO) considered in the study.The blue dotted lines indicated the Climate Regime Shifts (CRS).

Table 1 .
Location, depth, sampling periods and designs, climate and biological data sets used for the DFA at the five LTER sites.

Table 2 .
Model selection based on values of Akaike's information criterion (AIC).The optimal DynamicFactor Analysis model with one or two trends is given in bold.TEMP = summer bottom temperature, PDO = Pacific Decadal Oscillation index, SST = sea surface temperature, PRES = pressure, LAT = latitude, LONG = longitude, NAO = North Atlantic Oscillation index, NSE = North Sea Environmental index, AirTemp = air temperature, WeMO = Western Mediterranean Oscillation index.

Table 3 .
Measures of fit (diagonal elements of error covariance matrix) for the different time-series.Relatively low diagonal elements of the error covariance matrix (< 0.50) indicate good fit.AFDW = organic matter in the sediments as ash-free dry weight, CPE = sediment-bound chloroplastic pigments, PROT = particulate proteins in the sediments.

Table 4 .
Estimated t-values for the explanatory variables.Only values above 0.2 are shown.TEMP = summer bottom temperature, PDO = Pacific Decadal Oscillation index, SST = sea surface temperature, PRES = pressure, LAT = latitude, NSE = North Sea Environmental index, AirTemp = air temperature, WeMO = Western Mediterranean Oscillation index, PROT = particulate proteins in the sediments.