Changes in the Ecological Environment of the Marginal Seas along the Eurasian Continent from 2003 to 2014

: Based on time-series satellite-retrieved records of the marine ecological environment from Aqua/MODIS, we investigated changes in the sea surface temperature (SST), photosynthetically active radiation (PAR), Secchi Disc depth (SDD), and chlorophyll-a concentration (Chla) in 12 Eurasian marginal seas from 2003 to 2014. Results showed that the SST increased in all 12 marginal seas, with the enclosed marginal seas (i.e., Black Sea, Baltic Sea, Japan Sea, Mediterranean Sea, and Persian Gulf) exhibiting relatively higher rates of increase. The PAR generally decreased, except in the European marginal seas, though not signiﬁcantly. Similar to the changes in the SST, the SDD increased in all 12 marginal seas, with a maximum rate of 3.02%/year (or 0.25 m/year, p = 0.0003) found in the Persian Gulf. As expected, Chla generally decreased in the tropical marginal seas, but increased in the high-latitude marginal seas. The different relationships between SST and Chla changes indicate the complexity of global warming effects on marine phytoplankton in different marginal seas.


Introduction
The marginal seas along the Eurasian continent are major distribution regions of coral reefs and mangroves, as well as important seaweed production areas worldwide. The global coral reef area is estimated to cover 249,713 km 2 to 284,300 km 2 , with more than 90% located within the tropical areas of the Indian and West Pacific oceans [1,2]. The South China Sea, Java-Banda Sea, and coasts along the Indian Ocean are major distribution regions of coral reefs [3]. Moreover, the coasts along the South China Sea, Java-Banda Sea, and Bay of Bengal are among the most abundant regions for mangrove vegetation and seaweed production in global oceans [4,5]. Therefore, these marginal seas play key roles in the global oceanic ecosystem and carbon cycle.
Marine ecosystems are highly sensitive to environmental variation, especially long-term changes caused by climate change and human activities, which can markedly alter ecosystem health. Many studies have suggested that global phytoplankton biomass and primary production might be reduced with global warming, as inferred from satellite remote sensing records [6][7][8][9], coupled carbon-climate models [10], and century-long in situ water clarity data [11]. However, high spatial variability of such changes has also been observed, and different regions and time scales demonstrate inconsistencies and contrasting results [12][13][14][15][16], especially in the marginal seas where natural changes and anthropogenic effects can significantly affect the marine environment [17,18]. Thus, examining changes in the marine ecological environment of the marginal seas along the Eurasian continent is critical.
In the past two decades, based on time-series satellite ocean color remote sensing data, lots of studies have been carried out to investigate the changes in the chlorophyll concentrations and phytoplankton blooms in most of the marginal seas along the Eurasian continent, such as in the East China Sea [18], Japan Sea [19], Arabian Sea [20], Persian Gulf [21], Mediterranean Sea [22] and Baltic Sea [23]. However, most of them were focused on one or a few specific regions. In this study, we used time-series satellite-retrieved products from the Aqua/MODIS to investigate changes in the ecological environment of the marginal seas along the Eurasian continent from 2003 to 2014. Four environmental parameters related to marine ecosystems were retrieved from Aqua/MODIS and analyzed, and included sea surface temperature (SST), photosynthetically active radiation (PAR), Secchi Disc depth (SDD), and chlorophyll-a concentration (Chla). Moreover, the spatio-temporal variations and long-term changes in these parameters were compared among the marginal seas along the Eurasian continent.

Satellite Data Acquisition
We obtained monthly satellite products from MODIS onboard the Aqua satellite from 2003 to 2014, including SST, PAR, Chla, as well as remote sensing reflectance, water absorption coefficients, and particulate backscattering coefficients at 443 nm. These products were obtained from the NASA ocean color website (http://oceandata.sci.gsfc.nasa.gov). For the products of the water absorption and particulate backscattering coefficients, they were retrieved by the quasi-analytical algorithm developed by Lee et al. (2002). All products had a spatial resolution of~9 km with global coverage.

Satellite Retrieval of the SDD
The SDD is a good indicator of water transparency. Here, the SDD data were retrieved using the semi-analytic algorithm [24][25][26], as follows: where a and b b are the water absorption and backscattering coefficients, respectively; α and β are the refractive and reflective effects of the air-sea interface, respectively; ρ d is the reflectance of the Secchi disk in water with a value of 0.82 [27]; C e is the threshold contrast of the human eye-brain system with a value of 0.02 [27]; R is the irradiance reflectance just beneath the water surface. In Equation (1), only a, b b and R are variables, and R can be calculated from the satellite-retrieved remote sensing reflectance as follows: where Q is the bidirectional factor of the upwelling radiance just beneath the water surface, and R rs is the remote sensing reflectance. Based on Equations (1) and (2), we retrieved the SDD by satellite-derived remote sensing reflectance, water absorption coefficients, and particulate backscattering coefficients at 443 nm. The global SDD retrieved by the above semi-analytic algorithm has been validated previously using globally in situ samples which were located in the Northwestern Pacific Ocean, off the coasts of America, the North Sea, the Arctic Ocean and the Black Sea [26]. The validation results showed a correlation coefficient between satellite-retrieved and in situ SDD of 0.93, with a mean absolute relative error of 27.3% (for SDD ≥ 1.0 m) [26].

Statistical Analyses
To quantitatively analyze the changes in parameters from different regions, we partitioned the marginal seas around the Eurasian continent into 12 regions according to geographical locations, Sustainability 2018, 10, 635 3 of 15 as shown in Figure 1. These marginal seas included the Japan Sea, East China Sea, South China Sea, Java-Banda Sea, Bay of Bengal, Arabian Sea, Persian Gulf, Red Sea, Mediterranean Sea, Black Sea, North Sea, and Baltic Sea, with the marked numbers from 1 to 12 in Figure 1, respectively. Due to significant ice influence, the Arctic Ocean and Sea of Okhotsk were not included in this study.
Based on the monthly products of the SST, PAR, SDD and Chla over 2003-2014, we generated the climatology (or 12-years averaged) distribution maps for each parameter in the whole region, covering the 12 regions. Moreover, the mean values of the SST, PAR, SDD and Chla in each region were calculated from the climatology maps.
In addition, we also generated 12-year (2003-2014) averaged maps for each month to avoid contamination by abnormal values in certain months, and the difference between the maximal and minimal values was defined as the seasonal variation amplitude.
To analyze the changing rates in each region, we calculated the regional mean values of SST, PAR, SDD, and Chla in each month from 2003 to 2014. Based on the time series of the monthly regional mean values, we analyzed the absolute changing rates (the slope of the linear regression) for each region. Then, we obtained the relative changing rate (ratio of the slope of the linear regression to the regional mean value) for each parameter in each marginal sea.

Statistical Analyses
To quantitatively analyze the changes in parameters from different regions, we partitioned the marginal seas around the Eurasian continent into 12 regions according to geographical locations, as shown in Figure 1. These marginal seas included the Japan Sea, East China Sea, South China Sea, Java-Banda Sea, Bay of Bengal, Arabian Sea, Persian Gulf, Red Sea, Mediterranean Sea, Black Sea, North Sea, and Baltic Sea, with the marked numbers from 1 to 12 in Figure 1, respectively. Due to significant ice influence, the Arctic Ocean and Sea of Okhotsk were not included in this study.
Based on the monthly products of the SST, PAR, SDD and Chla over 2003-2014, we generated the climatology (or 12-years averaged) distribution maps for each parameter in the whole region, covering the 12 regions. Moreover, the mean values of the SST, PAR, SDD and Chla in each region were calculated from the climatology maps.
In addition, we also generated 12-year (2003-2014) averaged maps for each month to avoid contamination by abnormal values in certain months, and the difference between the maximal and minimal values was defined as the seasonal variation amplitude.
To analyze the changing rates in each region, we calculated the regional mean values of SST, PAR, SDD, and Chla in each month from 2003 to 2014. Based on the time series of the monthly regional mean values, we analyzed the absolute changing rates (the slope of the linear regression) for each region. Then, we obtained the relative changing rate (ratio of the slope of the linear regression to the regional mean value) for each parameter in each marginal sea.  Figure 2a shows the spatial distribution of the mean SST from 2003 to 2014. Influenced by the latitudinal variation of solar radiation, the SST generally decreased with increasing latitude. In some coastal regions, however, the SST was significantly lower than that in the same latitude oceans. For example, the mean SST in the East China Sea was about 5°C lower than that in the same latitude  Figure 2a shows the spatial distribution of the mean SST from 2003 to 2014. Influenced by the latitudinal variation of solar radiation, the SST generally decreased with increasing latitude. In some coastal regions, however, the SST was significantly lower than that in the same latitude oceans. For example, the mean SST in the East China Sea was about 5 • C lower than that in the same latitude oceans, likely due to southward coastal currents with cool waters in winter. In upwelling regions, the mean SST was slightly lower, such as observed in the Western Arabian Sea. oceans, likely due to southward coastal currents with cool waters in winter. In upwelling regions, the mean SST was slightly lower, such as observed in the Western Arabian Sea. Figure 2b shows the comparisons between mean SSTs in the 12 marginal seas. The mean SST values were 12.9 °C, 19   Japan Sea, and were larger than 10 • C in the Baltic Sea, North Sea, Mediterranean Sea, Black Sea, and Persian Gulf. Over 12 years (2003Over 12 years ( -2014, the SST increased in the 12 marginal seas, indicating oceanic warming around the Eurasian continent ( Figure 3b). However, the increasing rates showed significant regional inhomogeneity. The Black Sea exhibited the highest warming speed, with a rate of 1.23%/year (or 0.19 • C/year, p = 0.25), followed by the Baltic Sea, with a rate of 1.13%/year (or 0.10 • C/year, p = 0.56). The Japan Sea and East China Sea also demonstrated high warming speeds of 0.68%/year (or 0.09 • C/year, p = 0.57) and 0.47%/year (or 0.09 • C/year, p = 0.71), respectively. In the low latitude or tropical marginal seas, including the South China Sea, Java-Banda Sea, Bay of Bengal, and Arabian Sea, the warming speeds were very low, with rates of less than 0.1%/year (or 0.03 • C/year). It should be noted that although all significance values (p-values) in the 12 marginal seas were greater than 0.05, the rates still revealed net SST increases or decreases from 2003 to 2014. In contrast to the distribution of the mean SST, seasonal variation amplitudes in the SST showed high values in the high-latitude seas (Figure 3a), which approached 20 °C in the East China Sea and Japan Sea, and were larger than 10 °C in the Baltic Sea, North Sea, Mediterranean Sea, Black Sea, and Persian Gulf. Over 12 years (2003Over 12 years ( -2014, the SST increased in the 12 marginal seas, indicating oceanic warming around the Eurasian continent ( Figure 3b). However, the increasing rates showed significant regional inhomogeneity. The Black Sea exhibited the highest warming speed, with a rate of 1.23%/year (or 0.19 °C/year, p = 0.25), followed by the Baltic Sea, with a rate of 1.13%/year (or 0.10 °C/year, p = 0.56). The Japan Sea and East China Sea also demonstrated high warming speeds of 0.68%/year (or 0.09 °C/year, p = 0.57) and 0.47%/year (or 0.09 °C/year, p = 0.71), respectively. In the low latitude or tropical marginal seas, including the South China Sea, Java-Banda Sea, Bay of Bengal, and Arabian Sea, the warming speeds were very low, with rates of less than 0.1%/year (or 0.03 °C/year). It should be noted that although all significance values (p-values) in the 12 marginal seas were greater than 0.05, the rates still revealed net SST increases or decreases from 2003 to 2014.

Changes in PAR
Similar to the spatial distribution of the SST, the PAR also showed latitudinal variations (Figure 4a). Generally, the PAR decreased with increasing latitude; however, there were some regions in which the PAR values were lower or higher than those in the same latitude oceans, which may be caused by cloud coverage or aerosols. For example, the PAR was significantly low in the East China Sea, but relatively high in the Arabian Sea and Red Sea. Among the 12 marginal seas, the Red Sea had the highest PAR, whereas the North Sea had the lowest PAR, with a value about half that found for the Red Sea (Figure 4b).

Changes in PAR
Similar to the spatial distribution of the SST, the PAR also showed latitudinal variations (Figure 4a). Generally, the PAR decreased with increasing latitude; however, there were some regions in which the PAR values were lower or higher than those in the same latitude oceans, which may be caused by cloud coverage or aerosols. For example, the PAR was significantly low in the East China Sea, but relatively high in the Arabian Sea and Red Sea. Among the 12 marginal seas, the Red Sea had the highest PAR, whereas the North Sea had the lowest PAR, with a value about half that found for the Red Sea (Figure 4b). The high-latitude marginal seas had larger seasonal variation amplitudes than those in the lowlatitude seas (Figure 5a), which is likely related to the significant seasonal variation in solar radiation  The high-latitude marginal seas had larger seasonal variation amplitudes than those in the low-latitude seas (Figure 5a), which is likely related to the significant seasonal variation in solar radiation in high-latitude oceans. Among the 12 marginal seas, the Baltic Sea and North Sea demonstrated maximal seasonal variation amplitudes. Interestingly, the mean PAR in the Mediterranean Sea was significantly higher than that in the Japan Sea and East China Sea, despite their similar latitudes.
From 2003 to 2014, the PAR values decreased slightly in most marginal seas, except for a weak increase in the European marginal seas (Figure 5b). The tropical marginal seas (South China Sea, Bay of Bengal, and Java-Banda Sea) had high decreasing rates, with the maximum rate being −0.32%/year (or −0.14 [E/(m 2 d)]/year, p = 0.31) found in the South China Sea. The increases in PAR in the marginal seas around Europe were very weak. Overall, the changes in PAR were not significant for any marginal sea.

Changes in SDD
Compared to the latitudinal variations in SST and PAR, the spatial variations in SDD were relatively complex (Figure 6a) (Figure 6b). Due to the shallow water depths, the Baltic Sea had the lowest SDD, followed by the Persian Gulf. Conversely, the Mediterranean Sea had the highest mean SDD (more than 30 m), followed by the South China Sea and Bay of Bengal.  Figure 7 shows the seasonal variation amplitudes of the SDD across the whole region, covering the 12 regions. Notably, Figure 7 has much finer details than the maps of the mean SDD (Figure 6a). The reason for is that the mean value will smooth the image, while the seasonal variation amplitudes are the differences between maximal and minimal values, and will enhancement the image texture.   Figure 7 shows the seasonal variation amplitudes of the SDD across the whole region, covering the 12 regions. Notably, Figure 7 has much finer details than the maps of the mean SDD (Figure 6a). The reason for is that the mean value will smooth the image, while the seasonal variation amplitudes are the differences between maximal and minimal values, and will enhancement the image texture. Overall, the seasonal variation magnitudes of SDD in the nearshore regions are much smaller due to low SDD throughout the year. Similarly, the seasonal variation magnitudes of SDD in the oligotrophic basins are also small due to high SDD throughout the year. The high seasonal variation magnitudes of SDD mostly occur in the shelves or along the fronts between the shelves and basins, where the seasonal blooms and horizontal displacement of water masses are significant. Overall, the seasonal variation magnitudes of SDD in the nearshore regions are much smaller due to low SDD throughout the year. Similarly, the seasonal variation magnitudes of SDD in the oligotrophic basins are also small due to high SDD throughout the year. The high seasonal variation magnitudes of SDD mostly occur in the shelves or along the fronts between the shelves and basins, where the seasonal blooms and horizontal displacement of water masses are significant. For the seasonal variation amplitude in SDD, the Arabian Sea had the highest value among the 12 marginal seas, approaching 30 m (Figure 7a). With the effect of the seasonal monsoons, the mean SDD in the Arabian Sea can be up to 40 m in summer, while it decreases to around 10 m in winter, due to large-scale phytoplankton blooms [28]. The SDD in the Mediterranean Sea, Red Sea, Japan Sea, For the seasonal variation amplitude in SDD, the Arabian Sea had the highest value among the 12 marginal seas, approaching 30 m (Figure 7a). With the effect of the seasonal monsoons, the mean SDD in the Arabian Sea can be up to 40 m in summer, while it decreases to around 10 m in winter, due to large-scale phytoplankton blooms [28]. The SDD in the Mediterranean Sea, Red Sea, Japan Sea, Northern South China Sea, and Western Bay of Bengal also showed large seasonal variation amplitudes, with values up to 20 m. For the East China Sea, Baltic Sea, North Sea, Black Sea, and Persian Gulf, their small seasonal variation amplitudes of less than 5 m were due to the low SDD values for the whole year. In addition, in the Southern basin of the South China Sea and the Western basin of the Bay of Bengal, the seasonal variation amplitudes were also small, due to the high SDD values for the whole year.
From 2003 to 2014, the SDD in all 12 marginal seas decreased at different rates (Figure 7b). The Persian Gulf had the highest increase, with a rate of 3.02%/year (or 0.25 m/year, p = 0.0003). The Java-Banda Sea and Arabian Sea also had high increasing rates of around 1.3%/year (or 0.35 m/yearr, p < 0.007). The other nine marginal seas had increasing rates of less than 1.0%/year, with the South China Sea demonstrating the lowest rate of 0.18%/year (or 0.05 m/year, p = 0.49). The peak times of Chla differed among the 12 marginal seas, as shown in Figure 9. In high-latitude marginal seas, such as the Japan Sea, East China Sea, Baltic Sea and North Sea, the peak times generally occurred in April, corresponding to the spring bloom periods in these regions. However, in the tropical marginal seas, such as the South China Sea, Java-Banda Sea, Bay of Bengal and Persian Gulf, the peak Chla generally occurred in winter, due to cooling and strong monsoon-enhanced vertical mixing. Notably, the peak Chla occurred in February and August in the Arabian Sea and in August in the Black Sea. Figure 10a shows the spatial distribution of the seasonal variation amplitude of the Chla. The Baltic Sea, Arabian Sea, East China Sea, and Japan Sea showed large amplitudes. In the basins of the South China Sea, Bay of Bengal, and Mediterranean Sea, however, seasonal variation amplitudes were relatively low due to the low Chla throughout the whole year.

Discussion
In this study, we investigated the changes in SST, PAR, SDD, and Chla in 12 marginal seas along the Eurasian continent from 2003 to 2014. Results showed that SST increased in all 12 marginal seas, especially in the Black Sea, Baltic Sea, Japan Sea, and Mediterranean Sea. These marginal seas are all enclosed, with relatively long residence times and poor exchange with open ocean water. Furthermore, heat from global warming is more easily accumulated in these seas, and they are more easily affected by human activities [18]. Thus, compared to the open and semi-closed marginal seas, these enclosed marginal seas showed higher warming rates (≥0.5%/year or 0.09 • C/year).
Interestingly, results demonstrated opposite changes between PAR and SST in the 12 marginal seas. Changes in SST and PAR are generally connected. Enhanced PAR can increase water temperature by solar heating; however, increased SST can increase cloud coverage and reduce PAR. The general decrease in PAR with increasing SST indicates that the changes observed in SST were not primarily induced by changes in PAR, at least for the tropical marginal seas with decreasing PAR. However, in the high-latitude seas, such as the Baltic Sea and North Sea, whether the slight increase in PAR contributes to SST warming is worth future investigation.
From 2003 to 2014, the waters in the marginal seas along the Eurasian continent have become increasing clear, as indicated by the increasing SDD in all 12 marginal seas. In the tropical marginal seas, including the South China Sea, Bay of Bengal, Arabian Sea, and Red Sea, the increasing SDD should be related to the decreasing Chla [26]. However, in the enclosed marginal seas at high latitude, including the Baltic Sea, North Sea, Japan Sea, and East China Sea, the increasing SDD should be associated with decreasing terrestrial material. For example, sediment input from the Changjiang River into the East China Sea has been significantly reduced since the construction of the Three-Gorges Dam [29], whereas discharge and river plumes have been less impacted [30,31]. In fact, many river deltas worldwide have been retreating due to decreases in sediment supply from rivers [32].
Changes in Chla are more complex than changes in SDD as multiple factors can influence phytoplankton biomass and structure. Phytoplankton growth depends on an adequate supply of nutrients (e.g., nitrogen, phosphorus, silicon, and iron) and light in the euphotic layer. In tropical marginal seas, phytoplankton growth is often limited by low nutrient availability, especially in oligotrophic basins. Increasing SST values are generally associated with stronger stratification, which reduces upward nutrient transport from deeper layers and increases daily integrated exposure of phytoplankton cells within the upper mixed layer to ultraviolet (UV) and visible solar radiation. Thus, Chla will decrease with increasing SST in tropical marginal seas, such as that observed in the Arabian Sea, Java-Banda Sea, Red Sea, and Persian Gulf. The Chla decline in the tropical marginal seas could influence sustainable fisheries, as the Java-Banda Sea, South China Sea, and Bay of Bengal are among the richest regions for fish species worldwide [33].
In contrast to the overall decreasing Chla with warmer SST in the tropical marginal seas, Chla generally increased with increasing SST in the high-latitude seas. Phytoplankton growth is typically limited more by low light exposure than by low nutrient availability in high-latitude oceans. Warmer SST values can enhance phytoplankton biomass and primary production in a number of ways: (1) warming enhances enzymatic reactions, such as carboxylation catalysis by RuBisCO [34]; (2) warming encourages shoaling of the upper mixed layer, which, in turn, increases light availability for cells within the layer [35]; and (3) warming promotes earlier phytoplankton blooms and extended growing seasons. These combined effects were demonstrated in a 1958-2002 in situ time series from the Northeast Atlantic, wherein sea surface warming was accompanied by increasing phytoplankton abundance in cooler (annual mean SST < 12 • C) high-latitude regions and decreasing phytoplankton abundance in warmer (annual mean SST > 12 • C) low-latitude regions [36]. In addition, simulation results from atmosphere-ocean circulation models also show that sea surface warming is generally associated with increased marine primary production in high-latitude oceans because of increased light availability and a longer growing period [10]. In addition to the SST effects, the high-latitude marginal seas in the present study are mostly enclosed, with terrestrial nutrient input easily accumulated and phytoplankton growth consequently enhanced. For example, we previously observed double-bloom intensity in the Yellow Sea and Bohai Sea in the East China Sea from 1998 to 2011, due to a significant increase in nutrients in that region [18].

Conclusions
Based on the time series satellite products from the Aqua/MODIS, we investigated changes in the SST, PAR, SDD, and Chla in 12 marginal seas along the Eurasian continent from 2003 to 2014. Results showed that SST increased in all 12 marginal seas, though with different warming speeds. Among them, the enclosed marginal seas at high latitude, including the Black Sea, Baltic Sea, Japan Sea, and Mediterranean Sea, exhibited higher increasing rates than those observed for the tropical marginal seas. In contrast, PAR decreased slightly in most marginal seas, except for those in Europe. The increase in SST with decreasing PAR indicates that changes in the PAR did not induce the SST warming.
Similar to SST warming, the SDD in all 12 marginal seas increased from 2003 to 2014, indicating that the waters along the Eurasian continent have become increasingly clear. The decreasing SDD in the tropical marginal seas should be related to the decreasing Chla. In the enclosed marginal seas at high latitude, the increase in SDD is expected to be related to decreasing terrestrial material input.
As expected, the Chla generally declined in the tropical marginal seas, which can be linked to SST warming-enhanced stratification. However, Chla increased with warming SST in the high-latitude marginal seas, indicating complex effects of warming SST on phytoplankton growth in the different seas. Moreover, the accumulation of terrestrial nutrients in the enclosed marginal seas can also enhance phytoplankton growth.