Changing Physical Conditions and Lower and Upper Trophic Level Responses on the US Northeast Shelf

Friedland, Kevin D.; Morse, Ryan E.; Shackell, Nancy; Tam, Jamie C.; Morano, Janelle L.; Moisan, John R.; Brady, Damian C.

doi:10.3389/fmars.2020.567445

ORIGINAL RESEARCH article

Front. Mar. Sci., 15 September 2020

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.567445

This article is part of the Research Topic Managing for the Future: Understanding the Relative Roles of Climate and Fishing on Structure and Dynamics of Marine Ecosystems View all 22 articles

Changing Physical Conditions and Lower and Upper Trophic Level Responses on the US Northeast Shelf

$\r\nKevin D. Friedland*$ Kevin D. Friedland^1*

Ryan E. Morse¹

Janelle L. Morano³

¹Narragansett Laboratory, National Marine Fisheries Service, Narragansett, RI, United States
²Fisheries and Oceans Canada, Bedford Institute of Oceanography, Dartmouth, NS, Canada
³Department of Natural Resources, Cornell University, Ithaca, NY, United States
⁴Wallops Flight Facility, National Aeronautics and Space Administration (NASA) Goddard Space Flight Center, Wallops Island, VA, United States
⁵School of Marine Sciences, University of Maine, Walpole, ME, United States

Sea surface temperature (SST), salinity, and chlorophyll concentration (CHL) have changed in the US Northeast Shelf ecosystem over recent decades. The changes in these parameters were distinctly marked by change points around the year 2012 resulting in a 0.83°C increase in SST, a 0.3 PSU increase in salinity, and decrease in CHL in excess of 0.4 mg m^–3. Where temperature and salinity shifted in mean level around their respective change points, CHL declined in a more monotonic fashion. Modeled data suggest that the shift in CHL resulted in a greater contribution of pico- and nanophytoplankton and a decreased contribution of microphytoplankton to overall CHL. Complementary estimates of the contribution of different phytoplankton functional types suggest a diminished contribution of diatoms to the phytoplankton community. Hence, not only is there evidence of a decline in the overall primary production capacity of the ecosystem, but also evidence of a fundamental change in the size and quality of phytoplankton supporting food webs. Two ecosystem responses to the observed changes in SST, salinity, and CHL were analyzed. Both length and weight at age have declined for a number of species, and both measures of growth appear to be negatively associated with temperature and positively associated with CHL. Biomass of fish and macroinvertebrates has declined in recent years, with a decrease in pelagic species associated with a decrease in CHL, while the decline in demersal species was associated with an increase in temperature. Collectively, these ecosystem changes appear to be the result of the complex interactions of both thermal effects and changes at the base of the food web.

Introduction

Changing climate conditions are expected to impact the base of marine food webs by changing the seasonal timing and intensity of phytoplankton blooms (Lotze et al., 2019). The main mechanism governing this anticipated change is the shoaling of mixed layer depths due to increasing temperature, which in turn restricts vertical mixing (Somavilla et al., 2017) and inhibits the redistribution of dissolved nutrients leading to lower overall chlorophyll biomass (Henson et al., 2013). However, projections that incorporate physical mechanisms such as stratification and biogeochemical interactions require validation. Examination of changes in phytoplankton production in marine ecosystems subject to recent, rapid warming can provide the context to improve the modeled response of lower trophic levels to projected change in climate. Rapid climate change is already occurring and the associated warming is arguably the most evident change observed thus far (Cheng et al., 2019). Therefore, it is prudent to examine the response of both lower and upper trophic levels as proxies for the impacts of longer-term climate change effects.

Remote sensing data sources make it possible to consider contemporary change in the upper water column thermal dynamics and phytoplankton biomass in marine ecosystems. These data sources have matured to time series in excess of two decades (Groom et al., 2019) and now provide the ability to evaluate trends and shifts in ecosystem conditions (Friedland et al., 2018). Temperature can influence marine ecosystem function at both the system and individual levels. Increasing temperature can cause the metabolism in individuals to increase (Dantas et al., 2019), alter the population structure and phenology of phytoplankton communities (Moisan et al., 2002; Conversi et al., 2010), and result in changes in trait variation within populations (Salo et al., 2020). Consequently, marine ecosystems are reorganizing as temperature influences both range size and species richness (Batt et al., 2017). Increasingly, we see that where once fishing had the greatest influence on marine ecosystem structure and productivity (Shackell et al., 2012), climate effects are exacerbating forcing factors (Mérillet et al., 2020). The observational and remote sensing time series of temperature data predates estimates of chlorophyll biomass; hence, our thinking and investigations have tended to be skewed toward long-standing hypotheses of thermal control of marine ecosystem function. However, the base of the food web has not been static in response to climate factors resulting in changes in productivity (Roxy et al., 2016), phytoplankton community structure (Dutkiewicz et al., 2019), and bloom phenology (Friedland et al., 2018). While it is anticipated that climate change will alter the vertical structure of the water column, it is worth emphasizing that stratification has already undergone substantial change (Yamaguchi and Suga, 2019).

Change in the productivity of lower trophic levels can affect the stability and function of ecosystems and the production of surplus biomass to support fisheries. Ecosystem stability has been described as a function of energy flow (Huxel and McCann, 1998). For example, Ullah et al. (2018) described scenarios where trophic transfer is inhibited by thermally related factors and new equilibrium among phytoplankton functional groups are established. This can be particularly problematic when groups like cyanobacteria begin to dominate since they often tend to be refractory foods (Friedland et al., 2005) and some species pose a toxicity threat (Paerl, 2018). Perhaps of greater fundamental importance is how food web function changes if energy flow is modified in magnitude, timing, or becomes highly event driven. Terrestrial and aquatic species time their reproductive cycles to benefit from relatively constant seasonal production cycles. For example, cohort recruitment of fish species is affected by the timing and/or size of blooms (Asch et al., 2019). In much the same vein, subsequent growth and reproduction of cohorts will be impacted by the changing energy content and availability of forage species (Durant et al., 2019). The consequences of a loss of system stability and diversity all too often results in a concentration of fishing on lower trophic level target species, which tends to exacerbate the problem (Howarth et al., 2014). These perturbed ecosystems are often continental shelf large marine ecosystems, which play the dominant role in providing national and global food security compared to the relatively negligible role of high seas fisheries (Schiller et al., 2018). Hence, instability in production in marine systems can translate to instability in economic systems.

Climate change is expected to increase the frequency of extreme events, both transient and transformative in nature, that affect multiple aspects of marine ecosystems. Heatwaves have occurred worldwide and they are characterized as anomalously warm temperatures that persist and produce a myriad of follow-on effects (Holbrook et al., 2019; Pershing et al., 2019). Heatwaves can actuate change in community structure by stimulating emigration and causing regionalized mortality (Sanford et al., 2019). They have been associated with phytoplankton blooms and the development of hazardous blooms species that can cause fish mortality in other species (Roberts et al., 2019). Heatwaves can also impact fishery populations as they have been observed to play a role in the recruitment of both invertebrate (Chandrapavan et al., 2019) and vertebrate species (Smith et al., 2019).

The US Northeast Continental Shelf is a well-studied ecosystem that is already displaying warming induced changes in physical forcing and ecosystem response. It is a system that has experienced rapid temperature change (Pershing et al., 2015), has seen two prominent heatwaves in as many decades (Mills et al., 2013; Gawarkiewicz et al., 2019), and has experienced progressive warming in large measure due to change in basin scale circulation (Chen et al., 2020). Fish and macroinvertebrate species have undergone significant shifts in their distribution (Friedland et al., 2019), which is already impacting human populations (Rogers et al., 2019). Rising temperatures appears to have increased the available occupancy habitat for most species and concomitant with this expansion of habitat, there has been an increase in upper trophic level biomass across a range of functional groups (Friedland et al., 2020). The principal multi-species fisheries in the region were greatly reduced by over-fishing in the 1990s, and much of the economic space is now occupied by single species fisheries on invertebrate taxa (Goode et al., 2019; Wiedenmann and Jensen, 2019).

The aim of this study was to examine contemporary trends in sea surface temperature and chlorophyll concentration from remote sensing data sources, and salinity from observational data for the US Continental Shelf ecosystem. Furthermore, these trends are scrutinized for time series change points as an indication of potential regime change. The phytoplankton community was characterized by its size fraction composition, functional type contributions, and seasonal bloom patterns. Ecosystem response to both trend and events in the physical forcing and phytoplankton change was evaluated using growth and abundance data for fish and macroinvertebrates from a concomitant fisheries independent trawl survey.

Materials and Methods

Study Site

We studied physical and biological changes occurring in the Northeast US Continental Shelf ecosystem (NES). The shaded portion of the continental shelf shown in Figure 1 denotes the study area. We characterized conditions and responses of the ecosystem based on five subdivisions as described in Friedland et al. (2015a). The five subareas were Georges Bank (GBK), Gulf of Maine east and west (GOMe and GOMw, respectively), and Middle Atlantic Bight north and south (MABn and MABs, respectively). These five areas capture much of the variability of the system, especially as related to the distribution and dynamics of primary producers and the fish populations that rely on them.

FIGURE 1

Figure 1. US Northeast Shelf study area (NES) highlighted with shading and delimited to five subareas including: Georges Bank (GBK); Gulf of Maine east and west (GOMe and GOMw, respectively); and, Middle Atlantic Bight north and south (MABn and MABs, respectively), dashed line is the 100 m depth contour.

Sea Surface Temperature

High resolution sea surface temperature (SST) data for the ecosystem and subdivisions were sourced from the NOAA Optimum Interpolation (OISTT) 0.25 Degree Daily Sea Surface Temperature Analysis dataset (Reynolds et al., 2007). For the purposes of this study, we calculated monthly means from the daily data while retaining its spatial resolution (see “Data Availability Statement”).

Salinity

Salinity in practical salinity units (PSU) for the NES was sourced from shipboard data and estimated using an interpolation method described in Friedland et al. (2019). The main source of the data was ongoing resource and ecosystem surveys of the NES conducted by the Northeast Fisheries Science Center (NEFSC) of the National Marine Fisheries Service (see “Data Availability Statement”). Water column temperature and salinity have been collected contemporaneously with tows associated with the fall seasonal bottom trawl survey beginning in 1963 with the spring survey starting 5 years later (Desprespatanjo et al., 1988). In addition, data was also collected in two comprehensive ecosystem surveys over the study period including the Marine Resources Monitoring Assessment and Prediction program or MARMAP (1977–1987) and the Ecosystem Monitoring program or EcoMon (1992-present), both providing shelf-wide observations (Sherman et al., 1998; Kane, 2007). These data provide time series of spring and fall salinity in all the subareas. The spring and fall salinity time series were correlated (Spearman rank order correlation, p < 0.05); hence, we combined the seasonal time series to produce an annual PSU time series by taking the mean.

Chlorophyll Concentration

Chlorophyll concentration (CHL) data were extracted from a merged multi-sensor ocean color data product from the Hermes GlobColour (see “Data Availability Statement”). This product includes measurements made with the Sea-viewing Wide Field of View Sensor (SeaWiFS), Moderate Resolution Imaging Spectroradiometer on the Aqua satellite (MODIS), Medium Resolution Imaging Spectrometer (MERIS), and Visible and Infrared Imaging/Radiometer Suite (VIIRS) sensors during the period 1998–2019. The data were merged using the Garver, Siegel, Maritorena Model (GSM) algorithm (Maritorena et al., 2010). CHL was summarized as annual and seasonal means; the spring or first half of year seasonal was the average CHL from January through June; and, fall or second half of the year was the average CHL from July through December. The two halves of the year reflect the bloom patterns found in the NES that includes both spring and fall blooms (Friedland et al., 2016). Seasonal CHL time series were plotted as Z-scores (observation minus the mean and divided by the standard deviation) as were the phytoplankton size fraction and phytoplankton functional type data described below.

Time Series Change Points and Trend

We identified change in study time series as both discreet change points and trends, both of which may potentially indicate underlying processes such as a regime shift. We applied two parametric modeling approaches to identify potential change points in the annual means of SST, salinity, and CHL data. First, discrete change points or a change in signal level within a single time step were identified with the sequential averaging algorithm called STARS or “sequential t-test analysis of regime shifts” (Rodionov, 2004, 2006; Thomson and Emery, 2014). The STARS algorithm parameters were specified a priori to detect change in thermal, salinity, and chlorophyll regimes (alpha level α = 0.05; the length criteria was set to 10; and, the Huber weight was set to 1). To be consistent with the terminology associated with the STARS algorithm, the change in a parameter level before and after a change point is referred to as the change in regime means, which may have the units of°C, PSU, or mg m^–3. Second, we fit piecewise linear relationships using the R package “segmented” (version 1.1-0), which fits the data as segmented linear models while identifying time series break-points (Muggeo, 2003, 2017). The performance of these approaches was compared for an individual time series with Akaike information criterion (AIC) using the simplified form intended for model comparison.

A I C = 2 k + n l n (R S S)

Where RSS is the residual sum of squares, n is the number of observations, and k is the number of independent parameters (Burnham and Anderson, 2004). To evaluate trends in the data, we applied a non-parametric test of time series trend using the R package “zyp” (version 0.10–1.1). We used the Yue and Pilon method to estimate Theil-Sen slopes and performs an auto-correlation corrected Mann-Kendall test of trend (Yue et al., 2002).

Phytoplankton Size Fractions

Estimates of phytoplankton size fraction contribution to total CHL were based on the three-component model of phytoplankton size classes (Lamont et al., 2019) using the global parameters from Brewin et al. (2015). The model provides estimates of the contribution of microphytoplankton (> 20 μm), nanophytoplankton (2–20 μm), and picophytoplankton (< 2 μm) to total CHL. Model equations provide an estimate of the CHL fraction that is associated with combined contribution of pico- and nanoplankton:

C_{p, n} = C_{p, n}^{m} [1 - e x p (- \frac{D_{p}}{C_{p}^{m}}) C H L]

and picoplankton only:

C_{p} = C_{p}^{m} [1 - e x p (- \frac{D_{p}}{C_{p}^{m}}) C H L]

Which by difference yield the nano- and the microplankton fractions, Cn and Cm, respectively. The model was parameterized using the following global model estimates.

C_{p, n}^{m} = 0.77; C_{p}^{m} = 0.13; D_{p, n} = 0.94; D_{p} = 0.80

These parameters are similar to the average across regional and other global study estimates.

Phytoplankton Functional Types

Three complimentary approaches are considered for estimating Phytoplankton Functional Types (PFTs) to evaluate change in the quality of phytoplankton entering the food web of the NES. First, we used the method presented in Moisan et al. (2017) for studying spatial distributions of PFTs, their phenology, and ecotones along the U.S. east coast. That study utilized a total of 172 independent measurements of phytoplankton absorption spectra and High-Performance Liquid Chromatography (HPLC) pigment measurements to develop a satellite-based model for PFTs. Pigment-specific absorption spectra for 18 phytoplankton pigments were obtained by inverting the observed pigment concentrations and the total phytoplankton absorption at each wavelength using Singular Value Decomposition, SVD (Press et al., 1987). The resulting pigment-specific absorption spectra were then used with individual phytoplankton absorption spectra with the Non-Negative Least Squares, NNLS (Lawson and Hanson, 1995) inversion method to estimate the pigment concentrations. The individual phytoplankton absorption spectra values at each wavelength were modeled as a second order function with chlorophyll a as the independent variable.

In this study, the resulting SVD-derived pigment-specific absorption spectra and the wavelength-dependent, second order model coefficients for absorption spectra from Moisan et al. (2011) are used with the annual mean 1998–2019 chlorophyll a estimates for each of the five subareas to estimate the PFT time series. The chlorophyll a values are used with the second-order linear model to estimate the phytoplankton absorption spectra. These spectra are then used with the SVD-derived pigment-specific absorption spectra in the NNLS inverse process to estimate the 18 different phytoplankton pigments.

Once maps of the 18 phytoplankton pigments were derived, they were used to generate estimates of the various PFTs by using the estimation formulas outlined in Table 1 of Hirata et al. (2011) for diatoms, dinoflagellates, prymnesiophytes, prokaryotes, and green algae. The pigments necessary as inputs for these algorithms included: fucoxanthin, peridinin, chlorophyll-b, 19-butanoyloxyfucoxanthin, 19-hexanoyloxyfucoxanthin, alloxanthin, and zeaxanthin. A time series of these functional types were calculated for all five subareas.

The second approach used was an algorithm designed to estimate the dominant PFTs at a location. The method is called PHYSTWO and is described in Correa-Ramirez et al. (2018). The approach uses Empirical Orthogonal Function (EOF) decomposition to relate pigment signatures to water-leaving radiances. It applies normalized radiances at seven wavelengths including 412, 443, 469, 488, 531, 547, and 555 nm. The analysis was based on the monthly 25 km merged reflectance data retrieved from the same source as the CHL data. The PFT estimates were made for the period 2003 through 2019, which corresponds to the sensor collections for these wavelengths. Reflectance at 488 nm was estimated based on an interpolation of reflectance data for 469 and 490 nm. The PFT identified Nanoeucaryotes (typically chlorophytes and cryptophytes), Diatoms, Coccolithophorids, Phaeocystis-like (typically haptophytes), and two cyanobacteria, Prochloroccocus and Synechococcus. Following the approach used in Alvain et al. (2008), the percent contribution of a PFT was estimated as the ratio of the area of NES with each PFT to the total. Monthly means of the percent contributions of PFT were calculated over the course of each year.

Finally, as verification of modeled PFT data, we queried the abundance of total diatoms and dinoflagellates from the Continuous Plankton Recorder (CPR) Survey dataset. The CPR collects continuous measurements of zooplankton and phytoplankton taxa retained on the CPR mesh, hence it includes a partial sampling of the microphytoplankton (Batten et al., 2003). Monthly samples were collected in the NES during the period of 1998–2018, however, not all NES subareas were sampled over the entire period. The GBK and MABn subareas were sampled in all 21 years the time series and GOMe was sample in 20 years, however, GOMw and MABs were sampled in less than 7 years each. The average number of annual sample units was 138; a sample consisted of a total number of diatoms and dinoflagellates from 68 to 35 taxa, respectively (see “Data Availability Statement”).

Survey Growth Trends

Change in the size at age of fish in the NES was assessed using biological data collected during the Northeast Fisheries Science Center (NEFSC) bottom trawl survey (Desprespatanjo et al., 1988). Catch weight was collected during the full duration of the survey, however, individual weight measurements began in 1992 (see “Data Availability Statement”). From these data, changes in length and weight at age were examined. The size at age metrics were computed for a subset of survey species having sufficient age samples. Annual growth indices were computed using a general linear model (GLM) following the approaches applied to estimate catch per unit effort (Forrestal et al., 2019). The general form of the size at age GLM was:

G L M (L e n g t h / W e i g h t) \sim f a c t o r (S p e c i e s) + f a c t o r (A g e) + f a c t o r (Y e a r)

Moreover, it was also partitioned by season and sex. Year factor coefficients were used to represent the rate of change over time. The change in size at age was assessed with canonical correlation with temperature, salinity, and chlorophyll concentration as environmental covariates. The canonical correlation was fit using the CCA package in R (version 1.2).

Survey Biomass by Functional Groups

Seasonal time series trends in biota were represented by the CPUE for biomass of all taxa from the bottom trawl survey (see “Data Availability Statement”), and by assignment to functional groups based on their adult prey preferences and vertical distribution (Friedland et al., 2020). The functional groupings included benthivores, demersal piscivores, pelagic piscivores, and planktivores. Catches were standardized for various correction factors related to vessels and gears used in the time series (Miller et al., 2010). The spring and fall CPUE time series were correlated with Spearman rank order correlation reasoning that since many taxa are known to undertake seasonal migrations, correlated spring and fall signals could be considered an indication of more reliable abundance data. Similar to the change in size at age analysis, the change in biomass was assessed with canonical correlation with temperature, salinity, and chlorophyll concentration as environmental covariates.

Results

Trends in SST, PSU, and CHL

The main result of the trend analysis is that CHL declined in the Northeast Shelf as both SST and salinity increased, with change points in these data around the year 2012. We found that SST in the NES increased from ∼12 to 13°C in a stepwise fashion with a change point in 2012 (Figure 2A), which was also the highest SST in the time series. Salinity also increased in a stepwise fashion by ∼0.3 PSU, albeit in the year after the SST shift, noting there were no local maxima in salinity around the year 2012 (Figure 2B). The non-parametric trend test for NES SST yielded a significant positive Theil–Sen slope, however, the positive slope for salinity was non-significant at p = 0.1 (Table 1). In contrast, CHL decreased from ∼1.5 to 1.1 mg m^–3 in a segmented linear fashion, with a change point in 2012 (Figure 2C); and, the highest system-wide CHL level occurred in 2011. The Theil–Sen slope for NES CHL was negative and significant.

FIGURE 2

Figure 2. Mean annual SST (A,D,G,J,M,P) for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs, respectively), mean annual PSU (B,E,H,K,N,Q) and mean annual CHL (C,F,I,L,O,R) for the same areas, respectively. Lines are either STARS or segmented regression model fits; 2012 is highlighted with vertical dotted line.

TABLE 1

Table 1. Theil–Sen slope estimates expressed in decadal rate of change for environmental variables by NES and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs).

Similar patterns were found in the NES subareas time series of SST, salinity, and CHL. The pattern found for the NES as a whole was reflected in Georges Bank (Figures 2D–F) and MABn (Figures 2M–O). The SST change points in Gulf of Maine east (Figures 2G–I), and Gulf of Maine west (Figures 2J–L) subareas initiated earlier than the NES as a whole or in the neighboring Georges Bank area. In addition, the GOMe and GOMw salinity time series were modeled marginally better with a segment regression suggesting a less distinct transition in salinity. The segmented regression depiction of change in CHL in these areas was similar to the other parts of the ecosystem. The SST time series in the MABs was modeled best with segmented regression, and though 2012 was the warmest year in the time series, there was no evidence of a change in SST level as seen in the more northerly subareas (Figures 2P–R). The Theil–Sen slopes for the subareas follow the NES trend with the exception of Georges Bank SST, which had a positive trend that was non-significant. To summarize, the NES and its regional subareas experienced a shift in SST, salinity, and CHL around 2012. Both SST and salinity started increasing around this change point while CHL decreased, noting that the change in SST and salinity were generally more stepwise whereas the change in CHL was gradual.

Seasonal CHL

Seasonal CHL generally followed the patterns established in the annual means suggesting that spring and fall blooms have been affected by a similar set of factors. CHL for the first half of the year followed a coherent pattern of decline after 2012, similar to the CHL annual trend for the NES (Figure 3A). The Theil–Sen slopes for the first half of the year CHL were all negative and significant (Table 1). The CHL for the second half of the year were also coherent among the subareas, but had features that varied from the first half of the year (Figure 3B). In particular, localized minima around 2004 and maxima around 2011 varied between the data series. The Theil–Sen slopes for the second half of the year CHL were all negative and significant with exception of the data for Gulf of Maine east.

FIGURE 3

Figure 3. Z-score of CHL for the first (A) and second (B) half of the year for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs). The mean for each size class shown in the panel title with ± 1 SD; the NES time series is shown as a b-spline; 2012 is highlighted with vertical dotted line.

Phytoplankton Size Classes

The decrease in CHL throughout the ecosystem after 2012 was associated with a change in the size-class composition of the phytoplankton community. The estimates of the percent contribution to total CHL by picophytoplankton shifted from ∼10 to 14% after 2012, or a 40% increase in the contribution (Figure 4A). Similarly, nanophytoplankton exhibited a shift in composition from ∼23 to 28%, or a 22% increase in the contribution (Figure 4B). The increase in the contribution of smaller phytoplankton appears to come at the expense of microphytoplankton, which shifted from a contribution of 67 to 58%, or a 13% decline in contribution (Figure 4C). The Theil–Sen slopes for the trends in pico- and nanophytoplankton were all positive and significant; whereas, the trends for microphytoplankton were negative and significant (Table 1). In summary, the shift in NES CHL likely also represents a shift in the size class composition of the phytoplankton community driven by the increasing availability of smaller cells vs. larger cell taxa.

FIGURE 4

Figure 4. Z-score of estimated phytoplankton size fractional contribution percentage of pico (A), nano (B), and micro (C) plankton to total chlorophyll for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs). The mean percentage for each size class shown in the panel title with ± 1 SD; the NES time series is shown as a b-spline; 2012 is highlighted with vertical dotted line.

Phytoplankton Functional Types

Analysis of the PFT time series using the Moisan method suggests phytoplankton populations were relatively stable from 1998 to 2012. In the NES region, diatoms, dinoflagellates, green algae, prymnesiophytes and prokaryotes comprised nearly 63, 11, 13, 2, and 12 percent of the phytoplankton biomass in terms of chlorophyll a, respectively (Figures 5A–E). In 2012, all of these populations of phytoplankton show a dramatic linear trend away from the prior relatively stable levels. Specifically, the diatom population shows a decrease while all of the other functional types show increases in their relative percent biomass. These changes were consistent with a shift away from high production coastal phytoplankton populations to a more open ocean, less productive domain. These estimated changes were supported by the observed nearly 50% decrease in the chlorophyll a levels. These changes in phytoplankton composition were further corroborated by the phytoplankton size class estimates which show stable size class population levels prior to 2012 and decreasing microplankton with increasing nanoplankton and picoplankton levels after 2012, indicating a population shift away from a eutrophic coastal productive population to a more oligotrophic and less productive population. The Theil–Sen slopes for PFT in the analysis were significant and were positive trends except for diatoms, which had negative trends (Table 1).

FIGURE 5

Figure 5. Z-score of estimated percent contribution of PFT diatoms (A), dinoflagellates (B), green algae (C), prymnesiophytes (D), and prokaryotes (E) for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs). The mean percentage for each PFT shown in the panel title with ± 1 SD; the NES time series is shown as a b-spline; 2012 is highlighted with vertical dotted line.

Analysis of Phytoplankton Functional Types using the PHYSTWO method largely corroborated the changes identified by the Moisan method and had the advantage of identifying changes in the dominant PFT by subarea. For the NES, Nanoeucaryotes and Diatoms were the most important functional types with average percentages of 43% in all areas (Figures 6A,B). Phaeocystis-like forms, Coccolithophorids, and cyanobacteria were low contributors with mean percentages of only 8, 3, and > 1–2%, respectively (Figures 6C–F). The only consistent trend observed was a decrease in the contribution of diatoms, however, this trend was only significant in the NES and the Gulf of Maine subareas (Table 1). The other dominant group, Nanoeucaryotes, had a significant positive trend in the Gulf of Mine east and a negative trend in the Middle Atlantic Bight south.

FIGURE 6

Figure 6. Z-score of percentage area where Nanoeucaryotes (A), Diatoms (B), Phaeocystis-like (C), Coccolithophorids (D), Synechococcus (E), Prochloroccocus (F) were the dominant PFT for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs). The mean percentage for each PFT shown in the panel title with ± 1 SD; the NES time series is shown as a b-spline; 2012 is highlighted with vertical dotted line.

The time series CPR abundances support the conclusion that there was a change in diatom populations, but were inconclusive in respect to any change in dinoflagellates. Diatom abundances on a CPR sample unit averaged over 89,000 cells and in all the time series with sufficient data trended downward over the study period (Figure 7A). The Theil–Sen slopes were all negative for diatom trends and significant for the GBK and NES areas (Table 1). There were no apparent trends in the dinoflagellate time series (Figure 7B) and Theil–Sen slopes for these time series were of mixed sign and all non-significant. In summary, using two model-based methods to characterize Phytoplankton Functional Types, a shift was observed in functional types associated with the change in CHL in 2012 and most notably evidence to suggest a decline in diatoms in the plankton community. The conclusion of a change in diatom population based on the model-based data is in part supported by in situ sampling data.

FIGURE 7

Figure 7. Z-score of the abundance of diatoms (A) and dinoflagellates (B) in a CPR sample unit for the Northeast Shelf (NES) and NES subareas (NES, GBK, GOMe, GOMw, MABn, and MABs). The mean for each abundance shown in the panel title with ± 1 SD; the NES time series is shown as a b-spline; 2012 is highlighted with vertical dotted line.

Change in Fish Growth

Size at age has declined for a number of fish species in the NES. GLM coefficients indicate that length at age for males declined during spring (Figure 8A), however, the decline in female length at age was not significant (Table 2). Spring weight at age declined for both sexes (Figure 8B). Similar to spring, fall length at age for males showed a decline, whereas females did not decline significantly (Figure 8C). Fall weight at age declined for both sexes (Figure 8D). These GLM results are a composite trend of a number of important groundfish species. Two change points emerge from these data, one in the spring male length at age data in 2009 and in the fall male weight at age data in 2005. When we decomposed the analysis to the species level, we found that not all species had declining trends in growth. For example, two species, silver hake Merluccius bilinearis and butterfish Peprilus triacanthus, were consistently increasing in length and weight at age for both sexes and seasonal categories (Table 2). The canonical correlation analysis between environmental and growth variables for the composite fish sizes yielded two significant dimensions at p < 0.1 (Table 3). The growth variables were positively correlated with CHL and negatively related to temperature (Figure 9A); salinity produced weaker correlative relationships.

FIGURE 8

Figure 8. GLM model coefficients of length (A) and weight (B) for males and females of 15 species of fish during spring. GLM model coefficients of length (C) and weight (D) for males and females of 15 species of fish during fall. Lines are either STARS or segmented regression model fits; vertical dotted line marks 2012.

TABLE 2

Table 2. Theil–Sen slope estimates expressed in decadal rates of change for growth variables by season and sex for the results of GLM models and by species.

TABLE 3

Table 3. Tests of Canonical Dimensions for environmental and growth and biomass data.

FIGURE 9

Figure 9. Canonical correlation two-dimensional plot of environmental factors vs. growth (A) and biomass (B) variables. Environmental variable (shown in red) include temperature (TEMP), salinity (SAL), and CHL. Growth variables (shown in blue) are coded with first position as S for spring and F for fall, second position as M for males and F for females, and third positons as W for weight and L for length. Biomass variables (shown in blue) include benthivores (BEN), demersal piscivores (DEM), pelagic piscivores (PEL), planktivores (PLANK), and total biomass (TOTAL).

Change in Survey Biomass by Functional Groups

Time series of biomass constructed from scientific surveys are difficult to interpret since they are shaped by a number of different factors and often dependent upon the catch of migratory species. It would appear the total biomass for the system has increased over time (Figure 10E). The biomass of the functional groups, benthivores, pelagic piscivores, and planktivores, seem to suggest a recent decline in biomass (Figures 10A,C,D). In contrast, the biomass of demersal piscivores appears to have remained high, although the level of biomass has been highly variable from year to year (Figure 10B); this variability may be due to the influence of migratory species. Two functional groups, benthivores and planktivores, had spring and fall signals that were correlated suggesting the same populations were sampled between seasons. These time series show similar patterns, namely that biomass increased to high levels within the last decade and has been declining since around 2015. The benthivore time series was the only CPUE time series with the suggestion of a change point, which occurred in 2012. The analysis of overall trend suggests the increase in CPUE has been higher in the spring and that the increase in total biomass may be driven by changes in benthivore and demersal piscivore populations (Table 4). The canonical correlation analysis between environmental and biomass variables yielded two significant dimensions at p < 0.1 (Table 3). The biomass of benthivores and demersal piscivores were positively related to temperature whereas pelagic piscivores were correlated with CHL (Figure 9B); planktivores appear to be uncorrelated with these factors.

FIGURE 10

Figure 10. Biomass as catch per unit effort for benthivores (A), demersal piscivores (B), pelagic piscivores (C), planktivores (D), and total biomass (E) during spring and fall time frames. Dashed black line is mean of spring and fall curves. In each caption, r_s is the Spearman correlation with associated p-value. Dashed line are b-splines; solid line either STARS or segmented regression models; vertical dotted line marks 2012.

TABLE 4

Table 4. Theil–Sen slope estimates expressed in decadal rates of change for CPUE variables by functional groups.

Discussion

The US Northeast Shelf ecosystem has gone through a transformational change in thermal state, primary producer biomass levels, and functional group assemblages over the past two decades. The upper part of the water column increased in temperature around the year 2012 and has since remained at a temperature level nearly 1°C higher than prior to 2012. The year 2012 was an exceptionally warm year and is recognized as the warmest year in the historical record of SST measurements for the ecosystem (Chen et al., 2014). Salinity has also increased on the order of 0.3 PSU, indicating a change in advective source water in the upper water column and photic zone of the ecosystem (Austin et al., 2019). Coincident with this transition in temperature and salinity, chlorophyll biomass peaked in 2011, and then declined from a relatively stable pre-2011 level in excess of 1.5 mg m^–3 progressively to a low value of ∼1.0 mg m^–3. These broad scale ecosystem level changes were most evident in the northern segments of the ecosystem, specifically in the Gulf of Maine and on Georges Bank.

The link between climate-induced warming and CHL is extremely complex. If warming was primarily driven by internal heating processes, then we might expect that phytoplankton communities will shift to a warm tolerant assemblage over time (Barton et al., 2016). However, warming in the NES is often the result of changing offshore water masses entering the system and these water masses can have significantly different nutrient concentrations and ratios. For example, Townsend et al. (2015) documented more frequent incursions of low-nutrient Scotian shelf waters into the Gulf of Maine and Georges Bank over a ten year period. Shifts in the nature of incoming water masses may exert strong and unpredictable controls on NES productivity. However, stratification, another important driver of phytoplankton production, is generally expected to increase due to increases in freshwater inflow from the Arctic and associated warming (Greene, 2012; Li et al., 2015). Enhanced stratification can potentially reduce the depth of the photic zone, thereby increasing the light exposure of phytoplankton in the surface mixed layer, however, it can also disconnect surface phytoplankton from bottom waters limiting the availability of nutrients, particularly nitrate. The importance of vertical mixing for the delivery of bottom water nitrate for phytoplankton production is well documented in the NES (Townsend, 1998) and enhanced stratification induced uncoupling of the water column may be a reason for the 0.4 mg m^–3 reduction in CHL observed in this analysis since 2011. However, Li et al. (2015) points out that stratification dynamics in the Northwest Atlantic shelf systems can vary considerably with more haline control of stratification in the Gulf of Maine and thermal control in the Mid-Atlantic bight. A recent model analysis by Shin and Alexander (2019) using downscaled Global Climate Models contends that the rate of bottom water warming may outstrip surface water warming. Ultimately, a comprehensive analysis of stratification in recent high resolution downscaled models of the region coupled to the biogeochemical processes controlled by vertical mixing will be necessary step in projecting future climate induced changes in the productivity of the NES.

In addition to the change in chlorophyll biomass observed in recent years, evidence suggest there has been a dramatic change in the phytoplankton structure of the NES that portends ecosystem impacts (Flombaum et al., 2020). The size structure of the NES phytoplankton community has shifted to smaller sized cells, which due to functional trait-environment interactions and the saturation of productivity of these smaller cells would be expected to have a negative effect on primary production (Chen et al., 2019). The shift in cell size may be related to the pulsed uptake of resources during winter months just prior to the development of the spring bloom that favor small cells relative to larger celled phytoplankton (Lin et al., 2020). As systems become more oligotrophic there tends to be a shift to dominance by picophytoplankton that produces a size-dependent change in food quality associated with low grazing pressure (Branco et al., 2020). For example, it is well known that important secondary producers in the NES, such as Calanus finmarchicus and other large copepods, do not graze on cells less than 10 μm (Marshall and Orr, 1955; Frost, 1972; Bundy et al., 1998). In addition to Calanus finmarchicus (a keystone species in the Gulf of Maine), the filter feeding by commercially harvested bivalves in the NES declines when feeding on particle sizes below 10 microns as seen in blue mussels Mytilus edulis (Strohmeier et al., 2012), bay scallops Argopecten irradians and eastern oyster Crassostrea virginica (Palmer and Williams, 1980), and sea scallops Placopecten magellanicus (Brillant and MacDonald, 2000). The opposite occurs in nutrient-rich waters that favor larger phytoplankton species and higher grazing pressure. In the Southern Ocean, warm, stratified conditions in the surface waters are dominated by picophytoplankton and a zooplankton community dominated by small bodied crustaceans (Venkataramana et al., 2019). Aquatic food webs utilize greater than 50% of their primary production owning to the higher turnover rates associated with phytoplankton, making variation in phytoplankton production that much more critical to the production of secondary consumers and higher trophic level organisms (Barbier and Loreau, 2019).

In addition to an apparent shift in cell size fractions of phytoplankton, the dominant functional types of phytoplankton appear to have shifted, most notably a putative decline in diatoms. Diatoms play a pivotal role in marine food-webs shaping the productivity of ecosystems and are an important link in biogeochemical cycling (Harvey et al., 2019). Since diatoms are often chain forming and large-celled, a decline in diatoms may be synonymous with the size fraction shifts posited for the NES. And as with the factors associated with the development of size fraction patterns, vertical mixing is associated with the presence of larger phytoplankton like diatoms (Fragoso et al., 2019). The source phytoplankton that generates the flux of organic matter to the benthos is typically reflective of the dominant functional types in a given area on a global basis (Durkin et al., 2016), but since the NES phytoplankton community consists largely of diatoms, the contribution of diatoms to the POC is that much more critical. And since diatoms are passively buoyant, unlike other actively mobile flagellated phytoplankton, diatom blooms represent active transport of fixed carbon to the benthos (Gemmell et al., 2016). This effect is most pronounced with the fall bloom, which is known to produce enhanced rates of particulate organic carbon flux. Zooplankton grazing communities are in a transitional state in fall and can leave portions of aggressive blooms underutilized (Fujiwara et al., 2018). The fall period is also associated with changeover in the phytoplankton community leading to functional types that more readily tend to settle to the benthos (Kemp et al., 2000). We would expect that the Gulf of Maine and Georges Bank would be most impacted by the change in diatom populations since they more reliably have fall bloom activity (Friedland et al., 2015b). As with any set of modeling results, it is critical be circumspect about the reliability of the model fits and data output; hence, we also stress that the model results are at lease consistent with the observational counts of diatoms from the CPR samples.

The new levels of CHL and the associated potential decline in POC deposition in the NES may represent a stressor or critical production limitation for this ecosystem, and in particular, ecosystem services related to resource species. Though not well studied, it would appear large marine ecosystems may have a threshold of CHL associated with fishery yields. In LMEs with an average annual CHL over 1.0 mg m^–3, fishery total catch or fisheries production on an areal basis tend to be over twice the level of ecosystems with annual CHL less than 1.0 mg m^–3 (Friedland et al., 2012). Much of the NES is now at or below this threshold, which we would anticipate would result in lower productivity of higher trophic level organisms. In contemporary research, the most demonstrative effects of phytoplankton biomass and production relate to the feeding and condition of larval fish and its effect on year class strength. Over many decades in the North Sea, primary production measurements suggest a gradual decline in the base of the food chain and is implicated in the decline of secondary producers, notably small copepods including the genera Temora, Acartia, Pseudocalanus, and Paracalanus (Capuzzo et al., 2018). Small copepods often play a pivotal role in the first feeding of larval fish (Buckley and Durbin, 2006), hence the concern that recruitment levels have eroded as a consequence. Over shorter time series and with the benefit of remote sensing datasets, interannual variation in blooms and primary production levels have been associated with recruitment success of Pacific herring, Clupea pallasi, which seem to benefit from the proximity of the spring bloom to adult spawning (Boldt et al., 2019) and pre-metamorphic growth of rockfish (Sebastes spp.) larvae (Wheeler et al., 2017). Finally, in walleye pollock (Gadus chalcogrammus), whose recruitment is governed by thermally mediated predation effects (Uchiyama et al., 2020), bloom activity can serve as an event level forcing factor affecting year class strength (Gann et al., 2016).

It is difficult to predict the specific ecological response to any change in habitat (Friedland et al., 2020) due to the complexity of interactions and pressures in the region. Our observations indicate that there is an expansion of fish biomass in the NES, but at the same time, growth is slowing for many species. There are clear examples of the inverse relationship between fish growth and temperature, particularly at range edges (von Biela et al., 2015). Thus, while expanded thermal habitat is associated with higher recruitment, under the same temperature regime there is an expectation of a loss of larger fish as oxygen supply restricts increases in body size (Neuheimer et al., 2011). Long-term impacts of declining growth rates are also likely to reduce fecundity and egg quality, which varies with female fish size, and increases the risk of predation and starvation. For species such as silver hake Merluccius bilinearis that has shown a steady increase in growth over time in the NES, there is some evidence that they are able to supplement their diet with invertebrates such as Northern shrimp (Pandalus borealis) and avoid starvation (Link and Idoine, 2009). Thus, defining winners and losers under a thermal habitat regime shift, especially those coupled to changes in primary production, depends on more than just an individual species’ thermal tolerance; other ecological considerations need to be considered to fully understand and manage species under such changes to the habitat.

While the concept of regime change or shift in marine ecosystems has been extensively reviewed (Collie et al., 2004; Jiao, 2009), the definition of what constitutes a regime change varies by the ecosystem function under investigation. Regime change has been used in both formal and informal contexts. In the informal context, many practitioners simply associate a step change in conditions to be a regime change; and in many cases, the change in conditions is represented by a shift in a single factor (van Putten et al., 2019). However, in a more formal application of the term, a change of regime is considered a demonstration of a change in not only a suite of environmental indicators, but also evidence of a change in ecosystem function and productivity (Mollmann et al., 2015). In exploring the notion of high magnitude regime shifts, Scheffer and Carpenter (2003) considered what might cause a catastrophic shift in ecosystem state and conditions. They made the observation that a large change, or one with hysteresis, may be the result of accumulated gradual change that arrives at a threshold level and then the shift occurs. Another important distinction can be made between consumer and producer effects on transformative change (Connell et al., 2011). The producer level effects on the NES ecosystem state would appear to be dominant since the physical forcing and production level biology are associated with concomitant change points. Consumer effects would appear to be moderated by multiyear time lags in the sense the impacts of a change in growth and recruitment persists over the life cycle of the organism. This sharpens our focus on the shifts in temperature and salinity. The proximity of change on the NES to basin scale circulation (Kwon et al., 2019) would suggest the accumulated effects of climate perturbation resulted in a change point in ocean dynamics, which will likely be seen elsewhere beyond the confines of the NES. What should not be lost in the detail of the proximate response to change in physical and biological factors is that this climate change event on the NES simultaneously impacted the niche space of many species. What we do not know as of yet is whether the extent of niche change has or will actuate a change in biodiversity and fisheries productivity (Trisos et al., 2020).

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The data used in this analysis is available from the following sources: Sea surface temperature is available from the OISST website https://www.ncdc.noaa.gov/oisst; NEFSC survey salinity data is available from the National Center for Environmental Information https://www.nodc.noaa.gov/oads/stewardship/data_assets.html; CHL and water leaving radiance data is available from the Hermes GlobColour website http://hermes.acri.fr/index.php; NEFSC Survey length and weight at age and biomass is from the InPort NMFS Data Management Program https://inport.nmfs.noaa.gov/inport/; CPR data is from the Marine Biological survey site https://www.cprsurvey.org/

Author Contributions

KF led the analysis and drafting of the manuscript. RM and JRM assisted with data analysis and along with NS, JT, DB, and JLM contributed to the drafting of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

Funding for JM’s effort was provided by the NASA Ecosystem Forecasting Program, The Gordon and Betty Moore Foundation under Grant No. 3292, and the NASA Applied Information Systems Technology Program under NNH18ZDA001N-AIST.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank M. Correa-Ramirez and C. Melrose for help in estimating phytoplankton functional type data. J. Langan for comments on an early draft of the manuscript. P. Kostovick for assistance with data.

References

Alvain, S., Moulin, C., Dandonneau, Y., and Loisel, H. (2008). Seasonal distribution and succession of dominant phytoplankton groups in the global ocean: a satellite view. Glob. Biogeochem. Cycles 22:GB3001. doi: 10.1029/2007GB003154