Differences in Rate and Direction of Shifts between Phytoplankton Size Structure and Sea Surface Temperature

Species distributions are changing with various rates and directions in response to recent global warming. The velocity of sea surface temperature (SST) has been used to predict species migration and persistence as an expectation of how species track their thermal niches; however, several studies have found that evidence for species shifts has deviated from the velocity of SST. This study investigated whether estimation of the velocity of shifts in phytoplankton size structure using remote sensing data could contribute to better prediction of species shifts. A chlorophyll-a (Chla) size distribution (CSD) model was developed by quantifying the relationships between the size structure of the phytoplankton community and the spectral features of the phytoplankton absorption coefficient (aph(λ)), based on the principal component analysis approach. Model validation demonstrated that the exponent of CSD (hereafter, CSD slope), which can describe the synoptic size structure of a phytoplankton community, was derived successfully with a relative root mean square error of 18.5%. The median velocity of CSD slope across the ocean was 485.2 km·decade−1, broadly similar to Chla (531.5 km·decade−1). These values were twice the velocity of SST, and the directions of shifts in CSD slope and Chla were quite different from that of SST. Because Chla is generally covariant with the size structure of a phytoplankton community, we believe that spatiotemporal changes in Chla can explain the variations of phytoplankton size structure. Obvious differences in both rate and direction of shifts were found between the phytoplankton size structure and SST, implying that shifts of phytoplankton size structure could be a powerful tool for assessing the distributional shifts of marine species. Our results will contribute to generate global and regional maps of expected species shifts in response to environmental forcing.


Introduction
Global warming has brought dramatic changes to the marine environment [1].Physical and chemical variations have strong effects on the behavior of marine species [2], and responses of species to changing environmental conditions are manifest at the population level as shifts in distribution [3].Empirical and theoretical studies have suggested that while species tend to migrate toward higher latitudes [4,5], up to 60% of species have not responded as expected, and that shifts in species distributions have shown various rates and directions [6].
Improved understanding of how species shift their distributions is crucial for effective conservation under conditions of climate change [7].Climate velocity is derived as the ratio of the long-term temperature trend to the two-dimensional spatial gradient in temperature, and it enables predictions regarding species migration and persistence as an expectation of how species track the location of their thermal niches [8].In marine ecosystem research, climate velocity is calculated from sea surface temperature (SST); hence, we refer here to climate velocity in terms of SST.Some recent studies have exposed a complex mosaic of velocity of SST, suggesting non-uniform shifts in the global oceans with a wide range of rates and directions that deviate from simple poleward migration [9].Therefore, the velocity of SST is expected to describe observed species shifts [10,11]; however, unexpected shifts that are inconsistent with predictions generated by the velocity of SST have been reported for some species and in certain regions [12,13].
Phytoplankton plays numerous roles in marine ecosystems and biogeochemical processes.Consequently, information on the spatiotemporal pattern of the community composition and size structure is necessary to understand the ecological and biogeochemical functions [14].One of the most important functions of phytoplankton is to determine the energy transfer efficiency, which depends strongly on its size structure [15], and the size composition of a phytoplankton community is expected to affect species distribution through changes in the food web structure [16].These facts suggest that the size structure of the phytoplankton community is one of the most crucial factors affecting species distribution; hence, in addition to temperature, monitoring phytoplankton size structure is also likely to contribute to better predictions of species shifts.Although the velocity of SST has been used widely as a new index for describing expectations regarding species shifts [12,13], the velocity of shifts in phytoplankton size structures has not been reported yet.Therefore, we hypothesized that the velocity of shifts in phytoplankton size structures could contribute to the improvement of our understanding of species shifts that deviate from the warming perspective.
In the last decade, various methods have been developed for global estimation of phytoplankton size structure using satellite remote sensing based on optical properties such as the phytoplankton absorption coefficient (a ph (λ)) and the particle backscattering coefficient (b bp (λ)) [17][18][19].Kostadinov et al. [17] proposed a novel method to estimate phytoplankton size structure based on the particle size distribution (PSD) using the slope of b bp (λ).As phytoplankton is the main particle component in the pelagic ocean [20], PSD can in turn be related to the phytoplankton size distribution.The principal advantage of this approach is the arbitrariness of the arrangement of the size range, where other methods generally adopt a fixed target group or size class.However, b bp (λ) is affected not only by phytoplankton but also by non-algal particles (NAPs) [21].To overcome this problem, Roy et al. [19] reconstructed the PSD model using a chlorophyll-a (Chla)-specific absorption coefficient (a ph * (λ)) at the wavelength of 676 nm.They successfully estimated the phytoplankton size structure; however, the retrieval of inherent optical properties (IOPs) at longer wavelengths (such as 676 nm) is quite difficult because of the strong absorption by the water itself [22].Some recent studies have demonstrated the utility of the principal component analysis (PCA) approach in deriving information on phytoplankton community structure [23][24][25].For example, Wang et al. [25] developed a method to estimate phytoplankton size structure using the spectral features of normalized a ph (λ) captured by PCA.This method relies on the fact that the oscillation of loading factors can reveal the major factors driving a ph (λ) and its applicability.Thus, we believe that a PCA approach based on PSD could constitute a powerful tool for the global estimation of the synoptic size structure of the phytoplankton community.
In this study, we aimed to clarify the similarities and differences between the changes in phytoplankton size structure and SST.To achieve this objective, we developed a model for retrieving phytoplankton size structure by quantifying the phytoplankton size distribution using spectral features of a ph (λ).Then, the developed model was applied to satellite remote sensing data to estimate the phytoplankton size distribution.Furthermore, climate velocities derived from spatiotemporal changes in phytoplankton size distribution and SST were compared with regard to their contribution to better prediction of species distributional shifts occurring now and in the future.

In Situ Data
The dataset used in this study was obtained from a wide area following 23 cruises over a 16-year period (Figure 1, Table 1).At each station, size fractionated Chla (Chla size ), light absorption coefficient, and spectral radiation were measured, the detailed procedures of which are described below.Note that the spectral radiation was not observed at some stations.Seawater samples were collected from the sea surface (0 m) using a clean plastic bucket.Seventy percent of in situ data were used for model development and the remainder was reserved for model validation.We used various filters (Table 2) to obtain Chlasize during each cruise.Seawater samples (n = 364) were filtered onto polycarbonate membrane or nylon mesh filters (20, 10, 5, and 2 μm pore size) and GF/F filters (Whatman, ca.0.7 μm pore size) under low vacuum pressure (<0.013MPa).After filtration, the filter was immediately soaked in N,N-dimethylformamide, and pigments were

Absorption Coefficient of Phytoplankton
To measure phytoplankton absorption spectra, particles in water samples were collected onto a GF/F filter.The optical density (OD) of particles on the filter was measured immediately after sampling using a MPS-2400 or a MPS-2450 spectrophotometer (Shimadzu, Japan) at wavelengths of 350-750 nm, and the absorption coefficient of the total suspended particles (a p (λ)) was determined from the OD according to the quantitative filter technique [29].The OD of the NAPs was then measured after soaking in methanol [30] or sodium hypochlorite [31], and the absorption coefficient of the NAPs (a d (λ)) was quantified as well as a p (λ).Finally, a ph (λ) was obtained by subtracting a d (λ) from a p (λ).

Retrieval of IOPs from In Situ Spectral Radiation Data
In situ spectral distributions of the underwater radiation (n = 317) were acquired using a SMPR (Satlantic Inc., Halifax, NS, Canada), a HyperPro (Satlantic Inc., Halifax, NS, Canada) or a PRR-800/810 (Biospherical Instruments Inc., San Diego, CA, USA).The HyperPro measures underwater downward spectral irradiance (E d (λ, z)) and upward spectral radiance (L u (λ, z)) at 136 wavelengths between 350 and 800 nm.The SPMR and PRR-800/810 measures those radiations at 13 and 17 wavelengths between 400 and 705 nm and between 380 and 765 nm, respectively.These measurements were acquired at the same time as above-water measurements of downward spectral irradiance (E ds (λ)) were taken using a sensor mounted on the ship deck.

Satellite Remote Sensing Data
MODIS-Aqua Level-3 standard mapped daily images of OC3M Chla and R rs (λ) at 10 MODIS bands were downloaded from the Goddard Space Flight Center/National Aeronautics and Space Administration for 2003-2015.The temporal and spatial resolutions of these satellite data were daily and 9 km, respectively.The a ph (λ) at the MODIS bands were derived from R rs (λ) data using the QAA-v6.Daily Advanced Very High Resolution Radiometer Pathfinder Version 5.2 SST data with 4-km spatial resolution were obtained for 2003-2015 from the USA National Oceanographic Data Center/National Oceanic and Atmospheric Administration and the Group for High Resolution Sea Surface Temperature.The Pathfinder Version 5.2 data are an updated version of the Pathfinder Version 5.0 and 5.1 collections described in Casey et al. [34].These satellite data were resampled on a 1 • × 1 • grid and composited as monthly averages to minimize local and transient variabilities.

Model Description
Assuming that the PSD follows the Junge-type [35] power law size distribution [36], the number of particles per unit volume (N) of seawater normalized by the size bin diameter (D), N(D), can be expressed as follows: where ξ is the Junge slope of the PSD, and D 0 is a reference diameter at which N 0 = N(D 0 ).Therefore, the total number of particles in a given size range can be derived by integrating Equation (4) from the minimum diameter (D min ) to the maximum diameter (D max ); thus, enabling the PSD to be partitioned into distinct classes described by Given that the Chla is particulate, and thus the Chla size distribution (CSD) also follows the Junge-type power law distribution, the total Chla (Chla total ) and Chla size in a given size range from D 1 to D 2 can be expressed using Equations ( 6) and (7) as follows: where Chla 0 is the reference Chla at D 0 (here, 0.7 µm), and η is the exponent of the Chla size spectrum (hereafter, CSD slope).Larger magnitudes of the CSD slope indicate a large proportion of smaller phytoplankton, whereas smaller magnitudes suggest predominance of larger phytoplankton.In this study, we assumed D min and D max as the pore size of GF/F filter (i.e., nominal 0.7 µm) and 200 µm [37], respectively.The CSD slope is derived as the slope of the linear regression in log-space between the inverse log-transformed median diameters from log-transformed D 1 to D 2 and Chla size normalized by the size bin width.Note that the Chla fraction of an arbitrarily defined size range (F size ) can be derived using the CSD slope as follows: where the constants Chla 0 and D 0 no longer exist in Equation ( 8), such that only the CSD slope and diameter range are required for estimating each fraction of the phytoplankton size classes.

Quantification of CSD Slope Using Phytoplankton Absorption Spectra
To quantify the CSD slope using the spectral shape of a ph (λ), we applied PCA to the normalized a ph (λ), according to the method of Wang et al. [25].In brief, the normalized a ph (λ) (a std ph (λ)) was derived from the wavelength mean and standard deviation of a ph (λ) at 10 MODIS bands.The formula for a std ph (λ) is shown below: where mean(a ph (λ)) and std(a ph (λ)) are the wavelength mean and standard deviation of a ph (λ), respectively.Then, PCA was applied to a std ph (λ) to capture the spectral features of the phytoplankton absorption properties.The input values for the PCA comprised a matrix (m × N) constituted of a std ph (λ), where m and N are the numbers of the wavelengths and samples, respectively.The resulting principal component (PC) scores were assumed to correlate with the CSD slope; hence, the CSD slope was quantified as a logistic-type regression model using the i-th PC score (S i ) and regression coefficients between the CSD slope and PC scores (β 0 and β i ) as follows: (10) where k is the number of PCs.Here, S i can also be expressed as follows: where w i,j and a std ph (λ j ) are the loading factors for the i-th PC and a std ph (λ) value at wavelength j, respectively.Therefore, we obtained an equation by replacing S i in Equation (10) with Equation ( 11): Finally, the CSD slope is derived from Equation ( 12) using the model parameters of β 0 and C j .
Remote Sens. 2017, 9, 222 To avoid unrealistic values of the CSD slope, we defined the upper and lower limits of the CSD slope as 3.0 and 0.0, respectively.These limits mostly corresponded to 90% of pico-and microplankton dominance to total Chla.

Evaluation of Estimate Accuracy
The relative root mean square error (RRMSE) was adopted when we validated the agreement of two values, such as the measured and estimated CSD slope.The RRMSEs were computed as relative values to provide equal weights to all samples and then expressed as percentages [38,39], as described by where Meas i and Mod i are the i-th measured and modeled values, respectively.

Velocities of the Shift in Oceanographic Parameters
Velocities of the CSD slope, Chla, and SST were derived following the method of Burrows et al. [9].The velocity along angle θ (V θ ), with 0 • as North and 180 • as South, was calculated on a 1 • grid by dividing the temporal trend over the period 2003-2012 by the spatial gradient as follows: The long-term temporal trend was calculated as the slope of the linear regression between monthly anomalies and time.Note that the temporal trend was employed as an absolute value when we determined the velocity using Equation ( 15) to clarify the local speed of shifting contours.The spatial gradient was derived as the vector sum of the North-South and West-East spatial gradients of 3 × 3 neighborhoods, with the associated vector angle given by the direction of the gradient, as described by Spatial gradient = S NS cos θ + S EW sin θ (16) where S NS is the North-South spatial gradient and S EW is the East-West spatial gradient.The North-South spatial gradient was calculated as the difference in the temporal average for each northern and southern adjacent pair divided by the distance between them.The West-East spatial gradient was calculated similarly for each western and eastern adjacent pair divided by the distance between the cells in each pair.

Retrieval of the CSD Slope from In Situ Chla size
Various filters were used to obtain Chla size (Table 2).In general, phytoplankton is partitioned into three size classes: picoplankton (0.7-2.0 µm), nanoplankton (2.0-20.0µm) and microplankton (>20.0 µm).However, we have not used filters for dividing phytoplankton into the size classes during some cruises (e.g., UM0405, OS180 and KH-0904).This fact forces us to employ specific combinations of filters used for obtaining Chla size , resulting in only a part of our dataset can be used together.Thus, our dataset could not be used in its entirety with typical methods such as Brewin et al. [40] and Devred et al. [41] that need to fix the size range of each phytoplankton size class.In contrast, it is possible to derive the CSD slope from any combination of Chla size ; hence, we succeeded in developing a CSD model using the entire dataset.
To verify the robustness of the CSD slope values determined by various combinations of filters, the CSD slopes derived from different combination of Chla size were compared.For this comparison, the data obtained during the MR12-E03 cruise were employed because filters with various pore sizes were used during this cruise.Figure 2 shows the relationships between the CSD slope values derived from three typical Chla size (0.7-2.0, 2.0-20.0,and >20.0 µm) and different combination of Chla size .The resulting RRMSEs were 4.9% and 4.0%, and therefore we concluded that the differences in the CSD slope among the combinations of filters were negligible.
in developing a CSD model using the entire dataset.
To verify the robustness of the CSD slope values determined by various combinations of filters, the CSD slopes derived from different combination of Chlasize were compared.For this comparison, the data obtained during the MR12-E03 cruise were employed because filters with various pore sizes were used during this cruise.Figure 2 shows the relationships between the CSD slope values derived from three typical Chlasize (0.7-2.0, 2.0-20.0,and >20.0 μm) and different combination of Chlasize.The resulting RRMSEs were 4.9% and 4.0%, and therefore we concluded that the differences in the CSD slope among the combinations of filters were negligible.

Empirical Development of the CSD Model
The method described in Section 2.3 was applied to our dataset for developing the CSD model.The in situ CSD slope ranged from 0.417 to 2.427.Note that smaller CSD slope values indicate phytoplankton assemblages with larger sizes, and vice versa.Figure 3 shows the relationship between the CSD slope and spectral shape of (λ).It is clear that there are remarkable differences in (λ) accompanying the variation of CSD slope.For instance, smaller CSD slope values show remarkably higher values of (λ) between 400 and 430 nm than do larger values of the CSD slope.In addition, (λ) values associated with smaller CSD slope values are higher than those of larger CSD slope values between 660 and 690 nm.These differences in (λ) were mainly due to a coupling of the packaging effect and pigment composition, which is strongly affected by phytoplankton cell size and species composition [25].Overall, the results were generally consistent with Ciotti et al. [42] and Wang et al. [25].

Empirical Development of the CSD Model
The method described in Section 2.3 was applied to our dataset for developing the CSD model.The in situ CSD slope ranged from 0.417 to 2.427.Note that smaller CSD slope values indicate phytoplankton assemblages with larger sizes, and vice versa.Figure 3 shows the relationship between the CSD slope and spectral shape of a std ph (λ).It is clear that there are remarkable differences in a std ph (λ) accompanying the variation of CSD slope.For instance, smaller CSD slope values show remarkably higher values of a std ph (λ) between 400 and 430 nm than do larger values of the CSD slope.In addition, a std ph (λ) values associated with smaller CSD slope values are higher than those of larger CSD slope values between 660 and 690 nm.These differences in a std ph (λ) were mainly due to a coupling of the packaging effect and pigment composition, which is strongly affected by phytoplankton cell size and species composition [25].Overall, the results were generally consistent with Ciotti et al. [42] and Wang et al. [25].
To capture the spectral differences in a std ph (λ), PCA was conducted using a std ph (λ) between 400 and 700 nm with 1-nm resolution.The first four PC modes explained >94.2% of the total variance in a std ph (λ) and the first mode captured 52.1%.The spectral shape of the first mode indicates negative loading factor between 400 and 430, and 660 and 690 nm and positive loading factor between 430 and 550 nm (Figure 4a).The second mode explained 23.2% of the total variance and the spectrum shows an obvious negative peak between 660 and 690 nm (Figure 4b).The third mode explained 12.3% of the total variance in a std ph (λ).This mode shows a monotonous decrease from 400 to 430 nm, followed by an increase above 500 nm (Figure 4c).The fourth mode explained only 3.5% of the total variance; several positive and negative peaks are evident in the spectral shape (Figure 4d).To capture the spectral differences in (λ), PCA was conducted using (λ) between 400 and 700 nm with 1-nm resolution.The first four PC modes explained >94.2% of the total variance in (λ) and the first mode captured 52.1%.The spectral shape of the first mode indicates negative loading factor between 400 and 430, and 660 and 690 nm and positive loading factor between 430 and 550 nm (Figure 4a).The second mode explained 23.2% of the total variance and the spectrum shows an obvious negative peak between 660 and 690 nm (Figure 4b).The third mode explained 12.3% of the total variance in (λ).This mode shows a monotonous decrease from 400 to 430 nm, followed by an increase above 500 nm (Figure 4c).The fourth mode explained only 3.5% of the total variance; several positive and negative peaks are evident in the spectral shape (Figure 4d).To capture the spectral differences in (λ), PCA was conducted using (λ) between 400 and 700 nm with 1-nm resolution.The first four PC modes explained >94.2% of the total variance in (λ) and the first mode captured 52.1%.The spectral shape of the first mode indicates negative loading factor between 400 and 430, and 660 and 690 nm and positive loading factor between 430 and 550 nm (Figure 4a).The second mode explained 23.2% of the total variance and the spectrum shows an obvious negative peak between 660 and 690 nm (Figure 4b).The third mode explained 12.3% of the total variance in (λ).This mode shows a monotonous decrease from 400 to 430 nm, followed by an increase above 500 nm (Figure 4c).The fourth mode explained only 3.5% of the total variance; several positive and negative peaks are evident in the spectral shape (Figure 4d).These spectral features derived from PCA can be interpreted as signatures of changes in a std ph (λ).However, ocean color sensors restrict us to using only specific numbers of wavelengths for PCA because of the limitation of available observation bands.Although additional wavelengths are better for a PCA approach, previous studies have reported that reduced wavelength information has little influence on model performance [23,25].Initially, we attempted to conduct PCA for a std ph (λ) resampled at 10 MODIS bands between 400 and 700 nm; however, the estimation accuracy of the QAA-v6 became quite poor at longer wavelengths such as 645, 667, and 678 nm with huge RRMSE values of 67.2%, 77.9%, and 58.1%, respectively (Table 3).Therefore, only a ph (λ) at the seven shorter wavelengths of 412, 443, 469, 488, 531, 547, and 555 nm were used for developing our CSD model.The first four PC scores derived from a std ph (λ) at the seven wavelengths and the CSD slope were used to fit the logistic-type function, following which the model parameters β 0 and C j were determined.Comparison of the CSD slopes derived from in situ Chla size and in situ a ph (λ) shows good agreement (RRMSE: 17.4%); however, data scattering for higher CSD slope values is apparent (Figure 5a).This might be because one combination of regression coefficients cannot capture the entire relationship between the phytoplankton size structure and the spectral shape of a ph (λ).When focusing on the average a std ph (λ) at the MODIS bands, the a std ph (λ) at 412 nm was higher than at 469 nm for smaller CSD slopes, and vice versa.These spectral features were also found in previous studies, and they were attributed to differences in pigment composition related to phytoplankton cell size [42].For such practical considerations, we divided the dataset for model development into two, according to the relationship between a std ph (λ) at 412 and 469 nm, and then determined the model parameters for each dataset.Thus, the CSD model is defined as summarized in Table 4.The resulting relationship between in situ and modeled CSD slopes shows better agreement (Figure 5b), for which the RRMSE was reduced significantly to 13.2%.Note that sampling locations of these two dataset were randomly distributed in latitude and longitude, suggesting geographic bias seemed negligible.We also tried to quantify CSD slope based on previous methods described in Kostadinov et al. [17], Brewin et al. [40], and Fujiwara et al. [18], for which the resulting RRMSEs were 31.6%,31.7%, and 31.8%,respectively.These results suggest that the PCA approach is superior to other methods in the quantification of CSD slope.

Validation of the CSD Model
The performance of the CSD model was examined by comparing measured CSD slope values determined from in situ Chlasize and modeled CSD slope values estimated from in situ Rrs(λ).The resulting RRMSE between the measured and modeled CSD slope values was 18.5% (Figure 6).As the CSD model relies on the spectrum shape of aph(λ), the accuracy of the modeled CSD slope values

Validation of the CSD Model
The performance of the CSD model was examined by comparing measured CSD slope values determined from in situ Chla size and modeled CSD slope values estimated from in situ R rs (λ).
The resulting RRMSE between the measured and modeled CSD slope values was 18.5% (Figure 6).As the CSD model relies on the spectrum shape of a ph (λ), the accuracy of the modeled CSD slope values depends strongly on the robustness of the IOP algorithms.Therefore, the performance of the CSD model will be improved when IOP algorithms that are more robust are developed.However, the CSD model successfully retrieved the CSD slope within ±0.25 ranges for 75 out of 106 (70.8%) validated data.Thus, we considered that the CSD model performed with sufficient accuracy for assessing the CSD slope; hence, it could be used to derive global estimations of phytoplankton size structure.

Validation of the CSD Model
The performance of the CSD model was examined by comparing measured CSD slope values determined from in situ Chlasize and modeled CSD slope values estimated from in situ Rrs(λ).The resulting RRMSE between the measured and modeled CSD slope values was 18.5% (Figure 6).As the CSD model relies on the spectrum shape of aph(λ), the accuracy of the modeled CSD slope values depends strongly on the robustness of the IOP algorithms.Therefore, the performance of the CSD model will be improved when IOP algorithms that are more robust are developed.However, the CSD model successfully retrieved the CSD slope within ±0.25 ranges for 75 out of 106 (70.8%) validated data.Thus, we considered that the CSD model performed with sufficient accuracy for assessing the CSD slope; hence, it could be used to derive global estimations of phytoplankton size structure.12)) for CSD slope estimation.

Global Distribution of CSD Slope and Chla
Figure 7 shows the average CSD slope and Chla values during the period from January 2003 to December 2015.It can be seen that smaller CSD slope values are present in coastal upwelling regions such as in the areas of the Benguela and Humboldt currents, whereas larger CSD slope values are found in oligotrophic regions, especially the equatorial regions.In general, these global patterns of CSD slope exhibit an inverse relationship with Chla, suggesting areas of high and low productivity are associated with the relative dominance of larger and smaller phytoplankton assemblages, respectively.However, smaller CSD slope values are also evident in areas of low productivity, such as subtropical gyres.This is because the total absorption coefficient was close to that of pure seawater and the estimated IOPs include large errors, resulting in huge uncertainties regarding the performance of the CSD model.
are associated with the relative of larger and smaller phytoplankton assemblages, respectively.However, smaller CSD slope values are also evident in areas of low productivity, such as subtropical gyres.This is because the total absorption coefficient was close to that of pure seawater and the estimated IOPs include large errors, resulting in huge uncertainties regarding the performance of the CSD model.

Velocities of Shift in CSD Slope and Chla
The temporal trend of CSD slope values was derived as a simple linear regression slope of the monthly anomaly of CSD slope against time (Figure 8a).Overall, 20.1% of the global oceans had significant increases or decreases (p < 0.05) during 2003-2015.An increasing trend in CSD slope, indicating a change toward a smaller phytoplankton community, was found for 57.4% of these regions.A significant increasing or decreasing trend in Chla was found in 54.8% of the global oceans for 2003-2015 (Figure 9a), with 63.6% of these regions exhibiting a decreasing trend.The spatial gradient of CSD slope was relatively large along coastal regions and subtropical gyres result from dynamic changes in CSD slope with modest spatial displacements (Figure 8b).This global pattern of spatial gradient was consistent with that of Chla (Figure 9b).However, it should be noted that the large estimation error of CSD slope was substantial in areas of low productivity.As a clear temporal trend in Chla was found in these regions, the biogeochemical and ecological processes within these regions are likely to have experienced climate forcing.Therefore, accurate estimation of the phytoplankton size structure would contribute to improved understanding of the changes in these processes by providing information that is more detailed than provided by the spatiotemporal variability of Chla.
It appeared that the velocities of CSD slope and Chla depended strongly on their spatial gradients rather than on their temporal trends because the velocities showed clear inverse distributions with the spatial gradients.When the rates of temporal change and spatial gradient of the CSD slope were combined, the resulting median velocity of the CSD slope across the ocean was 485.2 km•decade −1 (Figure 8c).This value is largely consistent with that of 531.5 km•decade −1 for Chla (Figure 9c).These velocities are equivalent to the observed distribution shift of the phytoplankton community (469.9 ± 115.3 km•decade −1 ) reported by Poloczanska et al. [12].The directions of shifts both in CSD slope (Figure 8d) and in Chla (Figure 9d) show a mosaic distribution, suggesting that spatiotemporal variations of CSD slope and Chla are non-uniform across the oceans.As species response to climate forcing is specific to the regional scale [6], a complex mosaic of the spatiotemporal variations of CSD slope and Chla has the potential for providing assessments of shifts in species distributions that are more accurate.

Comparison of Velocities Derived from CSD Slope and Chla with That of SST
The SST shows significant temporal change in 34.9% of the global oceans, of which regions of warming and cooling account for 78.4% and 21.6%, respectively (Figure 10a).The global median velocity of SST was 292.8 km•decade −1 and this is approximately half that of the CSD slope and Chla (Figure 10c), suggesting a modest movement of the isotherms compared with the shifts in CSD slope and Chla.Note that the velocity of SST in this study was more than 10 times that reported in Burrows et al. [9] (21.7 km•decade −1 ).This discrepancy was caused by differences in methods; i.e.,

Comparison of Velocities Derived from CSD Slope and Chla with That of SST
The SST shows significant temporal change in 34.9% of the global oceans, of which regions of warming and cooling account for 78.4% and 21.6%, respectively (Figure 10a).The global median velocity of SST was 292.8 km•decade −1 and this is approximately half that of the CSD slope and Chla (Figure 10c), suggesting a modest movement of the isotherms compared with the shifts in CSD slope and Chla.Note that the velocity of SST in this study was more than 10 times that reported in Burrows et al. [9] (21.7 km•decade −1 ).This discrepancy was caused by differences in methods; i.e., this study used only statistically significant trends (p < 0.05), and non-significant trends that result in low velocities were excluded.In fact, the median velocity of SST in this study was 35.3 km•decade −1 when both significant and non-significant trends were included.
velocities among the CSD slope, Chla, and SST (Figure 11).For instance, the directions of shifts in CSD slope and Chla are largely consistent (Figure 11a), whereas the direction of shift in SST deviates from those of CSD slope (Figure 11b) and Chla (Figure 11c).The direction of shift in SST shows a patchy pattern in which some areas are comparable and others dissimilar to those of both CSD slope and Chla, and specific patterns were not found in this study.Thus, the rates and directions of shifts in CSD slope and Chla could explain the variations in species distribution that deviate from that of SST [12,13].Although the direction of shift in SST also shows a complex pattern, which clearly differs from a simple poleward migration (Figure 9d), there are clear regional differences in the directions of velocities among the CSD slope, Chla, and SST (Figure 11).For instance, the directions of shifts in CSD slope and Chla are largely consistent (Figure 11a), whereas the direction of shift in SST deviates from those of CSD slope (Figure 11b) and Chla (Figure 11c).The direction of shift in SST shows a patchy pattern in which some areas are comparable and others dissimilar to those of both CSD slope and Chla, and specific patterns were not found in this study.Thus, the rates and directions of shifts in CSD slope and Chla could explain the variations in species distribution that deviate from that of SST [12,13].

Discussion
Satellite-based estimation of spatiotemporal variation in phytoplankton size structure can contribute to better understanding of the biogeochemical and ecological responses to recent climate forcing [14].In this study, we developed a CSD model based on empirical relationships between the size distribution of Chla and the spectral shape of aph(λ).This method enabled us to derive with sufficient accuracy the synoptic size structure of the phytoplankton community using satellite remote sensing, based on our dataset obtained from a wide area of the North Pacific, Southern

Discussion
Satellite-based estimation of spatiotemporal variation in phytoplankton size structure can contribute to better understanding of the biogeochemical and ecological responses to recent climate forcing [14].In this study, we developed a CSD model based on empirical relationships between the size distribution of Chla and the spectral shape of a ph (λ).This method enabled us to derive with sufficient accuracy the synoptic size structure of the phytoplankton community using satellite remote sensing, based on our dataset obtained from a wide area of the North Pacific, Southern Ocean and Western Arctic oceans (Figure 1), despite approximately one third the samples being obtained from the optically complex waters of the Western Arctic Ocean [43,44].We should note that the global distribution of CSD slope includes uncertainties in clear waters due to the estimation error in a ph (λ) by QAA.Although other IOP algorithms for a ph (λ) estimation such as the Garver-Siegel-Maritorena model [45] and the Constrained Liner Matrix model [46] might perform well, even in the clearest waters, the fixed shape estimation of a ph (λ) is not applicable to our CSD model.Because the model relies strongly on the estimated spectrum of a ph (λ), which is derived from the water-leaving radiance through an IOP algorithm [25], the performance of the CSD model, especially in clear waters, will be improved when IOP algorithms that are more robust are developed.
Compared with existing models, our CSD model has several advantages.First, this model is applicable to both coastal and open waters.Although some models [47,48] use the spectral features of R rs (λ), which are affected strongly by non-algal components [20] and hence restricted to application in open waters, our CSD model has been developed using a ph (λ), which is separated from the absorption by colored dissolved organic matter and NAPs [22].Second, the CSD model can avoid the estimation errors of a ph (λ) because our method requires a ph (λ) at shorter wavelengths (i.e., between 412 and 555 mm), while the model proposed by Roy et al. [19] depends on the accuracy of the a ph (676) estimation, which generally contains significant error because of strong absorption by pure seawater [22].Third, the use of the spectral shape of a ph (λ) enables us to reduce the influence of non-algal components, whereas the b bp (λ)-based approach [17] can be affected not only by phytoplankton but also by all other suspended particles [21].In addition, the fraction of user-defined phytoplankton size classes can be derived with the CSD model according to Equation (8), while most other models can retrieve only predefined size classes.Furthermore, the CSD model has been developed using Chla size , which was determined by size-fractionated filtration.Although diagnostic pigment analysis (DPA) is frequently applied to classify phytoplankton size classes based on the photosynthetic pigment composition [19,25,39,40,49], Brewin et al. [50] reported that the DPA method had significant uncertainties.Therefore, we believe that the size-fractionated filtration-based method could be useful for obtaining phytoplankton size structures that are more accurate.Finally, the CSD model can represent the synoptic size structure of the phytoplankton community with one component, and hence the CSD slope is easily applicable to ocean biogeochemical models for retrieving phytoplankton size structures, which provide a useful index of the trophic status of the water mass, efficiency of the biological pump, and productivity of the marine ecosystem.
It has been well documented that the relationship between the size structure and abundance of the phytoplankton community shows strong covariance [39,40,49]; i.e., the fractional contribution of microplankton to Chla increases monotonically with increasing Chla, and vice versa.Indeed, the regional distribution of CSD slope showed a clear inverse pattern to that of Chla.Additionally, our demonstration of the temporal trend in Chla showed regional patterns of increase and decrease for the period 2003-2015, as reported in previous studies [51,52], and the distribution of CSD slope was generally the inverse.These results suggest that Chla and phytoplankton size structure are strongly covariant, and hence we consider that Chla could become an index of phytoplankton size structure.However, the regions where CSD slope and Chla showed statistically significant trends were quite limited.This is because the statistically significant regions of CSD slope were scarce in comparison with Chla.For instance, Chla showed significant increases or decreases in large parts of oligotrophic regions, where valid CSD slopes could not be retrieved because of the low estimation accuracy of the QAA-derived a ph (λ).Therefore, a robust a ph (λ) estimation in clear waters is an urgent requirement for the CSD model, as are other IOP-based approaches, for improved understanding of ecological and biogeochemical changes in these regions.
Species distributional shifts have been reported for a variety of organisms such as mollusks, fish, birds, and mammals [3,53], and their shifts have often been linked to temperature increases as expectations include simple poleward distribution shifts [54,55].However, in terms of a warming perspective, unexpected equatorward and shallower migrations have been observed [56,57]; thus, the complex mosaic of local velocity of SST (Figure 10c,d) is believed to represent a powerful tool for evaluating the variability in species distribution [9].Recent studies have applied the velocity of SST to the expectation of species range shifts [12,13].They found that some distributional shifts could be explained clearly by the velocity of SST, whereas others showed independent changes, attributable to the different sensitivities and strategies of marine organisms [11].VanDerWal et al. [10] reported that the expectation of shifts in species distributions driven by temperature alone is likely to underestimate the actual shifts.Moreover, Sunday et al. [58] revealed that the velocity of SST explained some of the variation in the range of shifts; however, including species traits more than doubled the explained variation.These facts suggest that adopting an approach that is more practical than the velocity of SST is required for better prediction of species shifts.
Trophic transfer efficiency describes the efficiency with which energy is transferred from one trophic level to the next [15] and it might control species distribution through changes in the food web structure [16].As the size structure of the phytoplankton community can have considerable effect in determining trophic transfer efficiency [15], the spatiotemporal variations of CSD slope and Chla are considered as important factors that could potentially cause species distribution shifts.In this study, both the rates and directions of the velocities of CSD slope and Chla showed clear regional differences with those of SST.These findings imply that the velocities of CSD slope and Chla could constitute powerful tools for predicting species distribution shifts that deviate from isothermal movement.Furthermore, an approach that combines the velocity of phytoplankton size structure with SST might contribute to predictions of species distribution shifts that are more accurate than existing approaches, which consider only variations in species thermal niches.
Numerous time series are currently being stored with the aim of monitoring the effects of climate change.Although species responses to climate change have been extremely variable across ocean regions and taxonomic groups [12,59], understanding the pattern of variation and identifying where and when species will respond to climate forcing is essential if we are to manage proactively changes in fisheries resources and to meet conservation goals [58].In consequence, as large proportions of species within marine ecosystems undergo distributional shifts [60], the provision of a cost-effective and relatively rapid integrated assessment of such shifts will be required [61].Therefore, techniques for assessing species distribution shifts are currently garnering considerable attention [62,63].In addition to the challenges involved in predicting species shifts, large gaps remain in our knowledge of species responses to climate change.This study provides compelling evidence regarding the potential for using phytoplankton size structure in assessing the persistence and migration of species distributions, which could contribute to the generation of global and regional maps of the expected rates and directions of species shifts that have improved accuracy.

Conclusions
This study proposed a CSD model to derive the synoptic size structure of the phytoplankton community using satellite remote sensing data.The CSD model is based on the relationship between phytoplankton size structure and its spectral absorption properties.The validation results showed that the CSD model performed with sufficient accuracy.Although satellite-derived CSD slope contains uncertainties, especially in clear waters, the CSD model has been demonstrated as a powerful tool for monitoring marine ecosystems through the synoptic size structure of the phytoplankton community.Recent evident deviations from uniform species shifts appear related to various environmental drivers, suggesting that a multifaceted approach, rather than using the velocity of SST, is probably required to assess species migration and persistence properly.Therefore, a new approach combining the velocity of phytoplankton size structure with other factors, including the velocity of SST, might contribute toward producing predictions of species distribution shifts with improved accuracy.Overall, our results could help generate global and regional maps of the expected rates and directions of species shifts.

Figure 1 .
Figure 1.Location of the sampling stations where in situ data obtained for this study through 22 cruises over the period 2000-2015.2.1.1.Size Fractionated Chla Concentration

Figure 1 .
Figure 1.Location of the sampling stations where in situ data obtained for this study through 23 cruises over the period 2000-2015.

Figure 2 .
Figure 2. The comparison of CSD slope derived from three typical Chlasize which were divided by the filters of 0.7, 2.0 and 20 μm pore size, and two different combinations of Chlasize: (a) Chlasize divided by the filters of 0.7, 2.0 and 10 μm pore size; and (b) 0.7, 5.0 and 20 μm pore size.Solid and dashed lined represent the 1:1 agreements and regression lines, respectively.

Figure 2 .
Figure 2. The comparison of CSD slope derived from three typical Chla size which were divided by the filters of 0.7, 2.0 and 20 µm pore size, and two different combinations of Chla size : (a) Chla size divided by the filters of 0.7, 2.0 and 10 µm pore size; and (b) 0.7, 5.0 and 20 µm pore size.Solid and dashed lined represent the 1:1 agreements and regression lines, respectively.

Figure 3 .
Figure 3.The relationship between (λ) spectra and CSD slope.Colored solid lines represent average (λ) for each range of CSD slope values; cool colors indicate larger size structure of phytoplankton, while warm colors indicate smaller size structure of phytoplankton.Black dashed lines represent the 7 MODIS bands between 412 and 555 nm, colored dotted lines indicate simple connecting lines between (λ) at 412 and 469 nm for each range of CSD slope values.

Figure 3 .
Figure 3.The relationship between a std ph (λ) spectra and CSD slope.Colored solid lines represent average a std ph (λ) for each range of CSD slope values; cool colors indicate larger size structure of phytoplankton, while warm colors indicate smaller size structure of phytoplankton.Black dashed lines represent the 7 MODIS bands between 412 and 555 nm, colored dotted lines indicate simple connecting lines between a std ph (λ) at 412 and 469 nm for each range of CSD slope values.

Figure 3 .
Figure 3.The relationship between (λ) spectra and CSD slope.Colored solid lines represent average (λ) for each range of CSD slope values; cool colors indicate larger size structure of phytoplankton, while warm colors indicate smaller size structure of phytoplankton.Black dashed lines represent the 7 MODIS bands between 412 and 555 nm, colored dotted lines indicate simple connecting lines between (λ) at 412 and 469 nm for each range of CSD slope values.

Figure 4 .
Figure 4. Loading factors of first four modes derived from (λ) spectra.CRs indicate cumulative proportion of variance for each mode.

Figure 4 .
Figure 4. (a-d) Loading factors of first four modes derived from a std ph (λ) spectra.CRs indicate cumulative proportion of variance for each mode.

Figure 5 .
Figure 5.Comparison between in situ and modeled CSD slope using the CSD slope: (a) before; and (b) after involving a condition branch.Solid line represents regression line, black and gray dashed lines indicate the 1:1 agreement and ±0.25 ranges, respectively.

Figure 5 .
Figure 5.Comparison between in situ and modeled CSD slope using the CSD slope: (a) before; and (b) after involving a condition branch.Solid line represents regression line, black and gray dashed lines indicate the 1:1 agreement and ±0.25 ranges, respectively.

Figure 6 .Figure 6 .
Figure 6.Comparison of in situ CSD slope determined from in situ Chlasize and modeled CSD slope derived from estimated (λ) through QAA-v6 using in situ PRR data.Solid line represents Figure 6.Comparison of in situ CSD slope determined from in situ Chla size and modeled CSD slope derived from estimated a std ph (λ) through QAA-v6 using in situ PRR data.Solid line represents regression line, and black and gray dashed lines indicate the 1:1 agreement and ±0.25 ranges, respectively.

Figure 7 .
Figure 7. Global distribution of: (a) CSD slope; and (b) Chla over 2003-2015.Cross-hatching shows areas of climatological aph(443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived aph(λ).

Figure 7 .
Figure 7. Global distribution of: (a) CSD slope; and (b) Chla over 2003-2015.Cross-hatching shows areas of climatological a ph (443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived a ph (λ).

23 Figure 8 .
Figure 8. Global distribution of: (a) temporal (b) spatial gradient; (c) velocity; and (d) direction of CSD slope.White shows the area with non-significant (p > 0.05) temporal trend for the period 2003-2015.Cross-hatching shows areas of climatological aph(443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived aph(λ).

Figure 8 .
Figure 8. Global distribution of: (a) temporal trend; (b) spatial gradient; (c) velocity; and (d) direction of CSD slope.White shows the area with non-significant (p > 0.05) temporal trend for the period 2003-2015.Cross-hatching shows areas of climatological a ph (443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived a ph (λ).

Figure 10 .
Figure 10.Global distribution of: (a) temporal trend; (b) spatial gradient; (c) velocity; and (d) direction of SST.White shows the area with non-significant (p > 0.05) temporal trend for the period 2003-2015.

Figure 10 .
Figure 10.Global distribution of: (a) temporal trend; (b) spatial gradient; (c) velocity; and (d) direction of SST.White shows the area with non-significant (p > 0.05) temporal trend for the period 2003-2015.

Figure 11 .
Figure 11.Distribution patterns of differences in direction of velocity between: (a) CSD slope and Chla; (b) CSD slope and SST; and (c) Chla and SST.White and gray areas show the area with non-significant (p > 0.05) trend in both and either of two velocities, respectively.Cross-hatching shows areas of climatological aph(443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived aph(λ).

Figure 11 .
Figure 11.Distribution patterns of differences in direction of velocity between: (a) CSD slope and Chla; (b) CSD slope and SST; and (c) Chla and SST.White and gray areas show the area with non-significant (p > 0.05) trend in both and either of two velocities, respectively.Cross-hatching shows areas of climatological a ph (443) = 0.01, where the CSD slope values contain uncertainties because of the low estimation accuracy of the QAA-derived a ph (λ).

Table 1 .
Cruises information and number of data obtained during each cruises.

Table 2 .
Summary of pore size of filters that were used to obtain Chla size during each cruises.

Table 3 .
Summary of performances of the QAA-v6 at 10 MODIS bands.R 2 and RRMSE were calculated in log10 space using full dataset (n = 317) after omitting minus values.

Table 4 .
Model parameters (β 0 and C j in Equation (