Environmental drivers of mesozooplankton ... - Wiley Online Library

PUBLICATIONS Journal of Geophysical Research: Biogeosciences RESEARCH ARTICLE 10.1002/2016JG003544 Key Points: • Fluctuations in primary production are the main drivers of variability of mesozooplankton biomass at station ALOHA • NPGO is the main climate pattern modulating productivity at station ALOHA • MEI and PDO effects on mesozooplankton biomass were weaker and less consistent

Supporting Information: • Supporting Information S1 Correspondence to: B. Valencia, [email protected]

Citation: Valencia, B., M. R. Landry, M. Décima, and C. C. S. Hannides (2016), Environmental drivers of mesozooplankton biomass variability in the North Pacific Subtropical Gyre, J. Geophys. Res. Biogeosci., 121, 3131–3143, doi:10.1002/2016JG003544. Received 3 JUL 2016 Accepted 16 NOV 2016 Accepted article online 21 NOV 2016 Published online 28 DEC 2016

Environmental drivers of mesozooplankton biomass variability in the North Pacific Subtropical Gyre Bellineth Valencia1

, Michael R. Landry1, Moira Décima2, and Cecelia C. S. Hannides3

1

Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California, USA, 2National Institute of Water and Atmospheric Research, NIWA, Wellington, New Zealand, 3Department of Geology and Geophysics, University of Hawaiˈi at Mānoa, Honolulu, Hawaii, USA

Abstract

The environmental drivers of zooplankton variability are poorly explored for the central subtropical Pacific, where a direct bottom-up food-web connection is suggested by increasing trends in primary production and mesozooplankton biomass at station ALOHA (A Long-term Oligotrophic Habitat Assessment) over the past 20 years (1994–2013). Here we use generalized additive models (GAMs) to investigate how these trends relate to the major modes of North Pacific climate variability. A GAM based on monthly mean data explains 43% of the temporal variability in mesozooplankton biomass with significant influences from primary productivity (PP), sea surface temperature (SST), North Pacific Gyre Oscillation (NPGO), and El Niño. This result mainly reflects the seasonal plankton cycle at station ALOHA, in which increasing light and SST lead to enhanced nitrogen fixation, productivity, and zooplankton biomass during summertime. Based on annual mean data, GAMs for two variables suggest that PP and 3–4 year lagged NPGO individually account for ~40% of zooplankton variability. The full annual mean GAM explains 70% of variability of zooplankton biomass with significant influences from PP, 4 year lagged NPGO, and 4 year lagged Pacific Decadal Oscillation (PDO). The NPGO affects wind stress, sea surface height, and subtropical gyre circulation and has been linked to mideuphotic zone anomalies in salinity and PP at station ALOHA. Our study broadens the known impact of this climate mode on plankton dynamics in the North Pacific. While lagged transport effects are also evident for subtropical waters, our study highlights a strong coupling between zooplankton fluctuations and PP, which differs from the transport-dominated climate influences that have been found for North Pacific boundary currents.

1. Introduction Time series data have been increasingly used over the past decade to shed new light on the impacts of natural cycles of environmental variability on marine plankton communities and to assess their potential vulnerabilities to climate change [Beaugrand and Reid, 2003; Chiba et al., 2006, 2012; Piontkovski et al., 2006; Ji et al., 2010; Garcia-Comas et al., 2011]. Zooplankton time series are particularly valuable in this regard because zooplankton relate to important food-web functions such as export and trophic transfer [Steinberg et al., 2012] and because the week-to-month generation times and motility of pelagic animals integrate over the high-frequency and fine spatial scales of environmental variability [Mackas and Beaugrand, 2010]. In the North Pacific Ocean, investigations of the environmental drivers of zooplankton variability have principally focused on the eastern and western boundary currents [Di Lorenzo et al., 2013]. In the California Current, abrupt changes in zooplankton biomass and community structure on interannual scales are strongly linked to fluctuations of El Niño–Southern Oscillation (ENSO) [e.g., Peterson et al., 2002; Lavaniegos and Ohman, 2007], while decadal-scale variability appears to be mainly modulated by fluctuations in transport processes associated with the Pacific Decadal Oscillation (PDO) [Peterson and Keister, 2003; Hooff and Peterson, 2006; Keister et al., 2011; Di Lorenzo and Ohman, 2013]. In the western Pacific boundary currents, transport processes also account for observed anomalies in zooplankton community [Chiba et al., 2013], but these are driven by lagged responses to the North Pacific Gyre Oscillation (NPGO), a climate mode that reflects variations in North Pacific wind stress and sea surface height [Di Lorenzo et al., 2008].

©2016. American Geophysical Union. All Rights Reserved.

VALENCIA ET AL.

Relative to boundary currents of the North Pacific, much less is known about plankton community responses to climate variability in the central oligotrophic region, the North Pacific Subtropical Gyre (NPSG). Station ALOHA (A Long-term Oligotrophic Habitat Assessment), sampled approximately monthly since 1988 by the

MESOZOOPLANKTON BIOMASS AT STATION ALOHA

3131

Journal of Geophysical Research: Biogeosciences

10.1002/2016JG003544

Hawaii Ocean Time-series (HOT) program, is the only regularly studied location in the vast NPSG region [Karl and Lukas, 1996]. The station occupies a position slightly north of the Hawaiian Islands (22.45°N, 158°W), where connections to ENSO, PDO, and NPGO modes of climate variability have been advanced to explain decadal-scale variability of plankton dynamics in the region. Corno et al. [2007] suggested that ENSO and PDO influences on upper ocean stratification and nutrient delivery may regulate primary production and phytoplankton community composition at station ALOHA. Bidigare et al. [2009] similarly argued based on a coupled physical-biological model that observed increases in productivity, photosynthetic eukaryotes, and zooplankton biomass from 1990 to 2004 were a consequence of shifts in PDO and ENSO states from positive to negative in the late 1990s, which reduced stratification and increased nutrient mixing. In contrast, Dave and Lozier [2010] did not find a relationship between stratification and primary production in the HOT data set, but they did note, as did Saba et al. [2010], a significant correlation between NPGO forcing and primary production, although weaker correlations with the PDO and the Multivariate ENSO Index (MEI) were also reported. Previous studies suggest a strong coupling of zooplankton biomass to primary production at station ALOHA [Sheridan and Landry, 2004; Chiaverano et al., 2013]. Seasonally, zooplankton biomass peaks are observed during the summer months of highest productivity [Landry et al., 2001]. A 9 year increase of mesozooplankton biomass from 1994 to 2002 [Sheridan and Landry, 2004] also coincided with the period of increasing productivity, as highlighted by Bidigare et al. [2009]. Based on these observations, we hypothesize that North Pacific climate modes drive mesozooplankton biomass variability in the NPSG principally through their bottom-up influences on primary productivity. In the present study, we develop generalized additive models (GAMs) to investigate this hypothesis using 20 years of monthly and annual mean zooplankton biomass data at station ALOHA. These analyses provide new insights on the contributions of climate variability to zooplankton fluctuations in the North Pacific.

2. Materials and Methods 2.1. Sampling and Laboratory Analysis Zooplankton samples were collected on 198 cruises from 1994 to 2013 at station ALOHA as part of the Hawaii Ocean Time-series (HOT) program. Oblique tows were taken through the euphotic zone (tow depth: 161 35 m, mean standard deviation (SD)) with a 1 m2 plankton net from February 1994 to August 2005 and with a 1 m diameter ring net from November 2005 to the present. Both nets used 202 μm Nitex mesh and were equipped with a General Oceanic flowmeter across the net mouth to measure volume-filtered and a time-depth recorder (Brancker XL-200 or Vemco logger) attached to the net frame to measure depth of tow. On most cruises, three samples were collected during midday (1000–1400) and three during the midnight (2200–0200). Onboard, a Folsom-split subsample (1/2, 1/4, or 1/8) from each tow was wet sieved through five mesh sizes (5, 2, 1, 0.5, and 0.2 mm), and each fraction was concentrated onto preweighed 200 μm Nitex filters, rinsed with isotonic ammonium formate to remove salt, and flash frozen in liquid nitrogen as described by Landry et al. [2001]. In the laboratory, dry weight (DW) biomass of each size fraction was determined (Denver Instrument analytical balance, 0.01 mg) after thawing and oven drying (60°C, at least 24 h). Weighed biomass was corrected for the volume of the subsample split, the volume of water filtered (m3), the tow depth (m), and is presented as areal biomass: g DW m2. Cruise DW means were calculated from the sums of the five size fractions in each tow and the average of day and night tows. Three cruises were excluded because both day and night samples were not collected. Carbon and nitrogen estimates were also determined as %DW for each size fraction on at least one day and one night tow per cruise by combusting a measured subsample of finely ground-dried zooplankton in a CHN elemental analyzer against known standards [Landry et al., 2001]. Here we use only the mean ratios to convert DW trends to approximate C equivalents. 2.2. Long-Term Trend Analysis Mesozooplankton biomass data (DW) were transformed (log10) prior to analysis to meet the assumption of normality (Shapiro-Wilk test, w = 0.99, p = 0.06). The biomass trend was evaluated by the generalized least squares (GLS) approach to account for positive autocorrelations of residuals, which in this case means removing the influence of past biomass on that measured at a given later time. In GLS, autocorrelation is

VALENCIA ET AL.


3132


10.1002/2016JG003544

modeled by including an autocorrelation structure that removes the independence restriction on residuals [Zuur et al., 2007, 2009]. The model was fit as log10 ðDWÞ ¼ a þ bX þ εm ;

(1)

where a and b are the model coefficients, intercept, and slope, respectively; X is the covariate; and εm are the residuals at time m modeled with an autoregressive process (AR1) as εm ¼ ρεm1 þ η;

(2)

where ρ is the correlation parameter between lagged residuals and η is an error term. Results of a GLS without (linear regression model) and with an autocorrelation structure AR1 were compared based on the Akaike information criterion (AIC), which measures goodness of fit and model complexity [Zuur et al., 2009]. The GLS model gave a lower AIC value than the linear regression model. Model validation was performed by evaluating the assumptions of homogeneity of variances, normality, and independence of data based on graphical analyses of the normalized residuals (Figure S1 in the supporting information). Fitted values and coefficients presented in figures and tables correspond to detransformed values. Data analysis was performed using R and the function “gls” in the nlme package [Pinheiro et al., 2015]. In addition to the trend analysis of mesozooplankton biomass, anomalies were calculated in order to evaluate 0

deviations from the long-term mean. For this, monthly biomass anomalies (Am ) were computed as 0 Am ¼ log10 Am log10 Ai ;

(3)

where Ām is the biomass for each month sampled (e.g., January 1998) and Ām is the long-term average biomass for month i (e.g., mean of all Januarys). Annual anomalies were calculated as the means of the monthly anomalies for each year [Mackas et al., 2001; OˈBrien et al., 2013]. 2.3. Environmental Data From the set of environmental data collected at station ALOHA (HOT website: http://hahana.soest.hawaii. edu/hot/hot-dogs/index.html), we used chlorophyll a concentrations from the upper euphotic zone (mg m2; integrated: 0–50 m) and the deep chlorophyll maximum (mg m2; integrated: 75–125 m), primary production (g C m2 d1; integrated: 0–150 m), and sea surface temperature (°C; mean: 0–10 m) as the locally measured variables most likely to be associated with temporal changes in mesozooplankton biomass. As in Saba et al. [2010], primary production was calculated as light minus dark bottle values from 1994 to 2000 and subsequently by subtracting 5% from the light bottle values. For basin-wide indices of ocean-atmosphere variability in the NPSG, we used the Multivariate ENSO Index (MEI), the North Pacific Gyre Oscillation (NPGO), and the Pacific Decadal Oscillation (PDO). The MEI monitors El Niño–Southern Oscillation based on changes over the tropical Pacific in sea level pressure, zonal and meridional components of surface wind, sea surface temperature, surface air temperature, and total cloudiness fraction. Positive MEI values indicate El Niño events, and negative values indicate La Niña events. Bimonthly mean values of the MEI were obtained from http://www.esrl.noaa.gov/psd/enso/mei/. The NPGO is defined as the second dominant mode of variability of sea surface height in the northeast Pacific and reflects changes in the intensity of the North Pacific gyre circulation [Di Lorenzo et al., 2008]. Positive values of the NPGO indicate stronger gyre circulation. Monthly mean values of the NPGO were obtained from http://www.o3d.org/npgo/. The PDO is defined as the first mode of North Pacific sea surface temperature variability [Mantua et al., 1997]. Monthly mean values of the PDO were obtained from http://research.jisao. washington.edu/pdo/. 2.4. Mesozooplankton Biomass and Environmental Factors We used generalized additive models (GAMs) to evaluate what factors might be driving the high- and low-frequency variabilities of mesozooplankton biomass at monthly and annual scales. The GAM is a nonlinear regression technique that fits smooth functions through data to model the relationships between response variables and covariates [Wood, 2006; Zuur et al., 2009]. Because GAMs do not assume a particular type of response function, they represent an effective modeling approach for assessing the responses of plankton communities to environmental factors [e.g., Irwin and Finkel, 2008; Llope et al., 2009; Otto et al., 2014].

VALENCIA ET AL.


3133


10.1002/2016JG003544

Table 1. Summary of Mesozooplankton Biomass, Primary Production, and Sea Surface Temperature at Station ALOHA a From 1994 to 2013 2

Mesooplankton biomass (g DW m ; integrated: 0 to ~160 m) 2 1 Primary production (g C m d ; integrated: 0–150 m) Sea surface temperature (°C; mean: 0–10 m)

b

n

Mean

SD

Min

Max

195 188 195

0.97 0.51 24.90

0.36 0.13 1.18

0.64 (1997) 0.41 (1998) 24.21 (1999)

1.39 (2011) 0.61 (2000) 25.54 (2004)

a

Data of environmental factors include only the months in which zooplankton samples were collected. Parenthesis indicate the years in which the minimum (min) and maximum (max) values were found. b SD denotes standard deviation.

Prior to GAM runs for the HOT data, collinearity between covariates was identified from the variance inflation factor (VIF) using a value of less than 3 as a cutoff for inclusion of covariates in the models. Although upper water-column stratification, computed as the change in potential density between 0–10 m and 140–150 m [Lavaniegos and Ohman, 2007], was initially considered as a covariate, it was highly collinear with SST and therefore was not included in the analyses. For monthly mean data, we modeled the effects of primary production (log10PP; g C m2 d1), sea surface temperature (SST; °C), MEI, NPGO, and PDO on mesozooplankton biomass (log10 DW) as log10 ðDWÞ ¼ a þ bð log10 PPÞ þ f 1 ðSSTÞ þ f 2 ðMEI:ml Þ þ f 3 ðNPGO:ml Þ þ f 4 ðPDO:ml Þ þ ε;

(4)

where a is the intercept, b is a linear function, fi are the thin-plate regression spline functions describing the effects of environmental factors on biomass, and the error term ε ≈ N(0, σ 2). The lagged response of the mesozooplankton biomass to the climate patterns denoted as l, was evaluated by running the models with no lag (m0) and with lags of 1 (m1), 2 (m2), or 3 (m3) months (Table S1 in the supporting information). Because mesozooplankton samples were not collected for all 12 months in most years, for each lag, the three indices were matched with the respective month in which zooplankton data were available. Annual means were calculated for mesozooplankton biomass and environmental factors to describe lowfrequency fluctuations. However, primary production and 1 year lagged NPGO (NPGO. y1) were highly collinear (VIF > 3) and thus not included in the same run (Table S1). Likewise, high collinearity was found between MEI and PDO of the same year (VIF > 3), and thus, two separate sets of models were run: log10 ðDWÞ ¼ a þ bðSSTÞ þ f 1 ð log10 PPÞ þ f 2 ðNPGO:yl Þ þ f 3 ðMEI:yl Þ þ ε;

(5)

log10 ðDWÞ ¼ a þ bðSSTÞ þ f 1 ð log10 PPÞ þ f 2 ðNPGO:yl Þ þ f 3 ðPDO:yl Þ þ ε;

(6)

Similar to the monthly mean analyses, the annual mean models were run with no lag (y0) and with lags of 1 (y1), 2 (y2), 3 (y3), or 4 (y4) years (Table S1). Although chlorophyll a concentration was included initially in both monthly and annual mean models, it was not significant in either case and is therefore left out of the results reported here. In addition to the full models for monthly and annual means, the effect of each covariate on mesozooplankton biomass variability was evaluated by running GAMs for each covariate separately. For both monthly and annual mean models, overfitting of the smooth functions was reduced by restricting the effective degrees of freedom (edf ≤ 4). For cases where edf was equal to 1 (linear relationship), models were rerun including the covariate as a linear term, as for example with primary production in the monthly mean models and sea surface temperature in the annual mean models. Full models were run for each combination of lagged indices, and a stepwise backward approach was applied in which nonsignificant covariates (p > 0.05) were removed. After the ensemble of models was obtained, the best monthly and annual models were selected based on minimizing the generalized cross-validation (GCV) criterion that measures the degree of penalization during fitting [Wood, 2006]. Model validation was done by graphical analyses of the residuals to evaluate the assumptions of homogeneity of variances, normality, and independence, as well as the model fit to observed values (Figures S2 and S3). Influential observations were evaluated by the Cookˈs distance (>0.5). These analyses were performed using R and the function “gam” in the mgcv package [Wood, 2006]. VALENCIA ET AL.


3134


10.1002/2016JG003544

3. Results 3.1. Temporal Variability in Mesozooplankton Biomass

2

Figure 1. Mesozooplankton biomass (g DW m ) at station ALOHA from years 1994 to 2013. (a) Mean of day-night samples per cruise for the combined size fractions. Long-term increase in biomass is described by generalized least squares; the straight line corresponds to the 2 detransformed predicted values (n = 195; intercept = 729 mg DW m ; 2 1 2 1 slope = 20 mg DW m yr ; 95% CI = 20–21 mg DW m yr ; p = 0.002). Curved fit is a three-point moving average. (b) Biomass trend as described by generalized additive models; the solid line is the model fit, and the dashed lines are the 95% confidence intervals. (c) Biomass anomalies obtained as annual means of the monthly anomalies using transformed (log10) dry weights.

Mesozooplankton biomass at station ALOHA averaged 0.97 g DW m2 (0.36 SD) and varied over 1 order of magnitude during the 20 year period analyzed (Table 1). Mean mesozooplankton biomass was 157% higher in 2013 (1.12 0.41 g DW m2; mean SD) than at the beginning of the time series in 1994 (0.65 0.24 g DW m2). Based on the GLS analysis, the long-term trend was significant with a mean rate of increase of 20 mg DW m2 yr1 (Figure 1a). To test if biomass has risen continuously over the study period, we ran a GAM of mesozooplankton biomass versus time that included an autocorrelation function AR1. Based on this analysis, mesozooplankton biomass increased rapidly during the first 10 years of the time series but has remained relatively constant since 2004 (Figure 1b). Annual biomass anomalies show negative values from 1994 to 1999 and mainly positive values since 2001, except during 2006 and 2009–2010. Strong negative anomalies like those registered from 1994 to 1997 have notably not occurred in the later time series data to date (Figure 1c). 3.2. Monthly Mean Model

The single function (two-variable GAM) of the monthly mean models that relates primary production, SST, NPGO, MEI, and PDO explained 22%, 14%, 11%, 8%, and 2% of mesozooplankton biomass variability, respectively (Table 2). Zooplankton biomass was not significantly associated with chlorophyll (Chl a) in the upper euphotic zone (0–50 m) or the deep chlorophyll maximum (Table 2). Based on the full GAM ensemble, the NPGO was the only climate index significantly related to mesozooplankton biomass, independent of a time lag up to 3 months (Table S2). The PDO was not significantly related to mesozooplankton biomass in any of the models, and the MEI was only significant when considered without lag (Table S2). The selected GAM for monthly mean data based on the GCV (and AIC > 2) indicates that local environmental factors, primary production, and SST, as well as basin-wide climate forcing, NPGO, and MEI (both indices without lag), are the most important environmental predictors of monthly changes in mesozooplankton biomass at station ALOHA (Table 3). Of these environmental factors, primary production accounts for about half of the variability explained by the model (Table 4). However, the final model explains less than 50% of total biomass variability (Table 3). The relationship between mesozooplankton biomass and environmental factors in the final GAM for monthly mean data range from linear and positive for primary production (Figure 2a) and SST (Figure 2b) to nonlinear for the NPGO and the MEI, showing a saturated response for the NPGO (Figure 2c) and a complex response

VALENCIA ET AL.


3135


10.1002/2016JG003544

Table 2. Two-Variable Generalized Additive Models With Monthly and Annual Mean Mesozooplankton Biomass at a Station ALOHA Over the Period of 1994–2013 Mean Monthly Model 2

Mean Annual Model

edf

R

DE (%)

GCV

P Value

Chl.0–50m Chl.max logPP SST PDO.m0 PDO.m1 PDO.m2 PDO.m3

1.9 1.0 1.0 1.2 1.0 1.0 1.0 2.8

0.01 0.00 0.22 0.14 0.02 0.00 0.00 0.02

1.8 0.7 22.0 14.7 2.7 0.6 0.1 3.2

0.0243 0.0244 0.0194 0.0209 0.0238 0.0244 0.0245 0.0242

0.299 0.239