Doliolid (Tunicata, Thaliacea) Blooms in the Southeastern Gulf of Alaska as a Result of the Recent Marine Heat Wave of 2014–2016

The eastern North Pacific experienced a prolonged heat wave in 2014–2016 manifested by high sea surface temperature anomalies in the south-central Gulf of Alaska (GOA). The event provided a natural experiment on the response of the southern GOA ecosystem to a dramatic change in sea temperature. Spatial and temporal variability in zooplankton communities following the culmination of the heat wave was investigated as a part of the NOAA Eastern GOA Ecosystem Assessment program in 2016–2017. Here, for the first time in the GOA, we report consistent observations of doliolid (Dolioletta tritonis) swarms observed in the upper mixed layer beyond the shelf break during both years, with the maximal density of 3,847 ind m–3 recorded in August 2016 and coinciding with the location of an offshore cyclonic mesoscale eddy. Doliolid density was significantly lower on the shelf. The long-term Continuous Plankton Recorder (CPR) data indicated that doliolid blooms in the south-central GOA may have occurred in the past two decades during El-Nino events. Coincidentally, doliolids prevailed in the diets of juvenile sablefish collected along the eastern coast of GOA both during the 2014–2016 heat wave and during 1997–1998 El Nino. Thus, we speculate that warming trends may increase the importance of doliolids in the GOA pelagic food web.


INTRODUCTION
Ecological impacts of climate change have become increasingly apparent affecting the organizational hierarchy of species and communities with marine ecosystems (Walther et al., 2002). In the recent decade, there has been an apparent rise in the frequency and persistence of marine heat waves (Oliver et al., 2018). The ecosystem impacts of these climate events can be severe resulting in shifts in geographical distributions and seasonal cycles and production of marine organisms, adversely affecting human activities (Uye, 2008;Purcell, 2012;Mills et al., 2015). Such extreme events can be viewed as climate change "stress-tests, " providing insight into possible changes in ecosystem structure, species' abundance and distribution, and the frequency and strength of ecosystem disruptions. In order to maintain resilient and sustainable marine ecosystems, it is imperative that management and conservation strategies are developed to mitigate and/or adapt to these impacts.
The North Pacific undergoes periodic influences of warmer periods typically resulting from large scale teleconnections of El Nino events. The 2014-2015 period was marked by unusually high sea level pressure which resulted in reduced wind mixing and advection and led to high near surface sea temperature anomalies in the south-central Gulf of Alaska (GOA) (termed "the warm blob") (Bond et al., 2015). The heat wave peaked in 2016 coinciding with the strong 2015-2016 El Niño (Newman et al., 2018;Walsh et al., 2018). After a brief cooling in 2017, temperatures increased again during 2018-2019, suggesting a prolonged climate event rather than a discrete heatwave (Litzow et al., 2020). The phenomenon attracted considerable attention as a potential future scenario for the changing climate. It has been speculated that the persistent elevated water temperatures may result in dramatic shifts in latitudinal distribution of marine organisms and structural reorganization of food webs along the North American Pacific coast, thus influencing commercial fisheries and impacting coastal communities. For example, warm temperatures in combination with other stressors have resulted in increase of abundance and production rates of gelatinous organisms such as hydrozoan and scyphozoan medusae, and pyrosomes (Purcell, 2012;Brodeur et al., 2019), and in decrease of commercially important fish such as Pacific Cod (Gadus macrocephalus) (Barbeaux et al., 2019).
The southeastern GOA coastal and bottom topography is defined by the Fairweather-Queen Charlotte fault system which forms the boundary between the convergent North American and Pacific tectonic plates (Brothers et al., 2018). The Queen Charlotte fault lies underwater in close vicinity to the narrow (5-10 km) continental shelf with depths<300 m and steep continental slope (Weingartner et al., 2009). The northern part of the boundary zone is adjacent to the submerged Yakutat Terrane which extends the shelf zone to ∼50 km offshore. The basin and continental slope circulation are dominated by a diffuse, weak eastern boundary Alaska Current with mean speeds of ∼5 cm s −1 and a width of several hundred kilometers (Weingartner et al., 2009; Figure 1). The Alaska Current diverges from the eastflowing North Pacific Current offshore from British Columbia and carries relatively warm and saline waters into the higher latitudes of the northern GOA. In contrast, the GOA shelf is characterized by a less saline, buoyancy driven Alaska Coastal Current flowing westward within 20-50 km of the shore (Royer, 1982;Weingartner et al., 2005). Mesoscale eddies are often generated along the eastern GOA continental slope leading to substantial exchange between the shelf and the basin, resulting in transport of heat, nutrients, microelements, and organisms across the slope which increase local productivity (Okkonen et al., 2003;Crawford et al., 2007;Ladd et al., 2007;Coyle et al., 2012).
Doliolids are planktonic tunicates with a complex life cycle, which includes polymorphic asexual zooid stages (called oozooids, trophozooids, and phorozooids) alternating with sexual solitary gonozooids (Deibel, 1998). The high asexual reproduction rate in combination with relatively short generation time governed primarily by temperature and food concentration facilitates rapid population growth under favorable conditions and often results in dense (thousands m −3 ) swarms of gonozooids typically restricted to the mixed layer above the pycnocline (Paffenhofer and Gibson, 1999;Deibel and Lowen, 2012;Takahashi et al., 2015). Doliolids are filter feeders capable of ingesting a wide size range of food Frontiers in Marine Science | www.frontiersin.org particles including microbes, packaging them into large fecal pellets and thus contributing to vertical carbon flux and a shunt of the microbial loop (Deibel and Paffenhofer, 2009). Their high individual feeding and excretion rates combined with explosive population growth may result in considerable impact on biogeochemical processes in the areas where the blooms occur. Doliolids seem to be adapted to event-scale fluctuations in the environment as opposed to seasonal-scale changes (Deibel and Lowen, 2012). The blooms are often associated with mesoscale physical events such as mid-shelf thermohaline fronts (Atkinson et al., 1978), eddies spinning off a strong boundary current and propagating along the adjacent shelf (Deibel and Paffenhofer, 2009), or meanders and filaments formed by convergent currents with contrasting physical properties (Takahashi et al., 2015). Such intrusions tend to bring nutrients into the euphotic zone fostering high rates of phytoplankton production, which, in turn, can be quickly exploited by doliolids (Paffenhofer and Lee, 1987;Deibel and Paffenhofer, 2009). The patchy and sporadic nature of such events makes doliolid blooms a challenging topic to study.
The sibling species Dolioletta tritonis (Herdman 1888) and Dolioletta gegenbauri (Uljanin 1884) occur in tropical and temperate waters of the Atlantic, Pacific, and Indian oceans and often penetrate as far north as the Faroe Islands following the warm summer transgression (Godeaux, 2003). However, until recently, dense aggregations were found only within the warmer tropical and subtropical core of their biogeographical range (Mackas et al., 1991;Deibel and Paffenhofer, 2009;Takahashi et al., 2015). In this paper, we report the first observation of intense D. tritonis bloom in the southeastern GOA during the 2016-2017 heat wave and utilize long-term data sets and historical data to link doliolid emergence in the southeastern GOA to climate variability.

Net Sampling
Samples were collected during NOAA Eastern GOA Ecosystem Assessment surveys in summer 2016 and 2017 (Figure 1). Sampling within the US Exclusive Economic Zone (EEZ) was done along 10 main parallel transects run perpendicular to the coast and supplemental stations spaced 10 to 18-37 km apart. An additional rectangular sampling grid consisted of offshore stations spaced 37 km apart. The nearshore surveys were conducted between July 4 and 28 followed by the offshore sampling continued from August 4 through 10. A total of 72 stations were sampled each year.
Conductivity-temperature-depth (CTD) measurements were taken with a Sea-Bird Electronics Inc. (SBE) 911 or SBE 25 CTD at each station to determine vertical temperature and conductivity profiles. Casts did not exceed 200 m (if bottom depths were>200 m) and were limited to the upper 50 m while sampling the offshore grid.
Seawater samples were collected at five depths (0, 10, 20, 40, and 50 m) to determine total chlorophyll-a (chl-a) concentrations (using Whatman GF/F filters). Filters were stored at -70 • C and analyzed with a Turner Designs (TD-700) laboratory fluorometer (Parsons et al., 1984). In situ fluorometric data were calibrated with discrete chl-asamples to estimate chl-a concentrations (r 2 = 0.678). The satellite image [Visible Infrared Imaging Radiometer Suite(VIIRS) on the Suomi-NPP] of the northern GOA in June 2016 was obtained from https://earthobservatory. nasa.gov/images/88238/bloom-in-the-gulf-of-alaska. No clear satellite images of the study area were available for the summer 2017 because of extensive cloudiness.
Doliolids were sampled with a 60 cm MARMAP-style bongo frame with a 505 µm mesh net. The net frame was equipped with an SBE 49 CTD transmitting real time tow data. Oblique tows from within 5-10 m off the bottom or from 200 m (at depths>200 m) to the surface were conducted primarily during daylight hours. Only the upper 50 m layer was sampled at the additional offshore grid. Volume filtered was measured with calibrated General Oceanics flowmeters mounted inside the net. All samples were preserved in 5% formalin, buffered with seawater for later processing. In the laboratory, the samples were split with a Folsom splitter and the subsamples were examined to identify and count doliolids. Collections from 2016 were sorted at the Polish Plankton Sorting and Identification Center (Szczecin, Poland) following NOAA protocols which omit identification of doliolid stages and mass (e.g., Eisner et al., 2018). Collections from 2017 were sorted at the University of Alaska (College of Fisheries and Ocean Sciences) in Juneau, AK, United States. Wet mass of doliolids were determined for 2017 samples. Blotted preserved wet mass was measured on a Cahn Electrobalance or Mettler top-loading balance, depending on the size of the doliolids. The wet mass was converted to carbon (C) content using 0.0072 coefficient established for Tunicata (Kiørboe, 2013). The data were uploaded to a Microsoft Access database, and analyses were performed using STATISTICA 6 and R software (R Core Team, 2015).

Continuous Plankton Recorder Sampling
The continuous plankton recorder (CPR) was towed behind commercial ships across the Gulf of Alaska on routes from Alaska to the North American West Coast, from March through October each year,and from the West Coast towards Asia year-around, from 2000 to 2017. Several vessels conducted the sampling over the time series but in each case the CPR was towed in the wake of the ship at a depth of about 7 m.
The details on CPR sampling are provided elsewhere (Batten et al., 2003(Batten et al., , 2018. Briefly, water and associated plankton entered the front of the CPR through a small square aperture (1.27 cm 2 ), and then through silk filtering mesh (with a mesh size of 270 µm) which retained the plankton and allowed the water to exit at the back of the machine. The movement of the CPR through the water turned an external propeller which, via a drive shaft and gearbox, moved the filtering mesh across the tunnel at a rate of approximately 10 cm per 18.5 km of tow. As the filtering mesh left the tunnel it was covered by a second band of mesh so that the plankton were sandwiched between these two layers. The mesh layers were then wound into a storage chamber containing buffered 40% formaldehyde preservative (which diluted in the seawater to a concentration of about 4%, sufficient to fix and preserve the plankton). The towed mesh was processed according to standard CPR protocols (Batten et al., 2003(Batten et al., , 2018. It was first cut into separate samples (each representing 18.5 km of tow and about 3 m 3 of seawater filtered) and every 4th sample was randomly apportioned amongst the analysts for plankton identification and counting. The ship's log was used to determine the mid-point latitude and longitude of each sample along with the date and time.

Data Analysis
The stations were split into two groups: (1) the shelf group with<300 m depth and (2) the offshore group with >300 m depth. Since the shelf slope is narrow and steep in the study area, the majority of the offshore stations were located at >1,000 m depths. Depth of the thermocline was computed for each station by locating the depth where dT/dz was maximum (T = temperature, • C; z = depth, m). The mean water-column temperatures above and below the pycnocline, as well as throughout the upper 50 m water column were then computed.
Since only the upper 50 m were sampled on the deep offshore stations, doliolid abundances estimated from other stations were adjusted by proportion of the total volume filtered in the upper 50 m based on the assumption that doliolid aggregations typically occur at shallower depths near the thermocline (e.g., Paffenhofer et al., 1991;Takahashi et al., 2015). Because of uneven spatial distribution of doliolids, both abundance and directly measured biomass data were fourth-root power transformed to stabilize the variance and reduce heteroscedasticity. Analysis of variance (ANOVA) was used to test for significant effects of location and year on physical and biological variables, and the distribution of residuals was analyzed to ensure it meets the normality assumption. Interpolation of physical and biological values over the sampling area was done on non-transformed data using kriging subroutines in SURFER 12 (Golden Software).
Station-based doliolid abundance data were used to examine how D. tritonis responded to biophysical environmental factors during summer in 2016-2017. Generalized Additive Models (GAMs; Hastie and Tibshirani, 1990;Wood, 2006) were then used to identify significant relationships between doliolid abundance and water properties that had significant effects on doliolid distribution: where g is a log link, µ is the expectation of observations, s is a smoothing function, ε is an error term, and k 1 -k 4 are the numbers of knots for smoothing functions. We used bottom depth, and the mean temperature, salinity, and nonfractioned chl-a concentration integrated through the upper 50 m at sampling locations. GAMs were run using the mgcv library in R, version 3.2.0 (R Core Team, 2015), with a quasi-Poisson family and cubic regression splines as the smoothing function of predictor variables. A subset of predictor variables that produced the best-fit models was selected using generalized cross validation (GCV) methods (Wood, 2006). Significance of the predictor variables was assessed by using Chi-square statistics calculated by R (R Core Team, 2015). In addition, partial regression plots showing the additive effect of each predictor variable on the doliolid abundance were examined for linearity, significance, and positive and negative effect. Residuals from best-fit models were checked for the assumption of independence and identical distributions. The resulting best-fit model were used to determine which predictor variables may have a significant effect on the doliolid abundance.
The time series of monthly average Pacific Decadal Oscillation (PDO) index anomalies available from the Joint Institute for the Study of the Atmosphere and Ocean 1 was used to define cold and warm phases on the Gulf of Alaska during the past two decades. We divided the years into positive, negative, or neutral using 0.5 or -0.5 annual mean PDO values from March through October as the determining values ( Table 1). The time series of monthly Ocean Nino index (ONI) was obtained from the National Climate Prediction Center. 2

RESULTS
In both 2016 and 2017, a two-layered water column structure marked by a pycnocline located at ∼15 m depth was observed throughout the study area (Figure 2 and Table 2). The upper While mean temperature below the pycnocline did not change significantly between years or regions, the upper layer was ∼1.5 • C warmer in 2016 than in 2017 (ANOVA, F = 70.8, p < 0.0001). On average, the entire upper 50 m layer over the slope and basin was ∼1 • C warmer in 2016 than in 2017 (ANOVA, F = 44.7, p < 0.0001), while no interannual difference was found on the shelf. The warm core was located offshore beyond the shelf break during both years (Figure 3). While offshore salinity did not differ between the years, a brackish riverine plume was observed during both years in vicinity of Alsek River, spreading along the coast and significantly lowering salinity above the pycnocline by 1.5 PSU (ANOVA, F = 18.8, p < 0.0001). Similar to temperature, mean chl-a concentration was also significantly higher offshore in 2016 than in 2017 (ANOVA, F = 7.4, p < 0.01). A VIIRS satellite image taken three weeks before the 2016 cruise revealed a mesoscale anti-cyclonic eddy ∼400 km offshore in the study area marked by elevated chlorophyll concentrations at the eddy periphery (Figure 4) indicating higher primary production and likely better feeding conditions for zooplankton grazers in 2016. Dolioletta tritonis were significantly more abundant at offshore stations in the basin than over the shelf (ANOVA, F = 60.3, p < 0.0001) ( Table 2). No significant interannual differences in doliolid abundance were observed within each station group. Dense (∼100-800 ind m −3 ) aggregations of doliolids were observed well beyond the shelf break at deep offshore stations (>3,000 m) in both years (Figure 5). The maximum abundance of 3,847 ind m −3 was recorded in 2016 coinciding with the position of the offshore eddy observed a month earlier. In 2017, doliolid biomass was significantly higher offshore than over the shelf (Table 2), ranging from 0.02 to 2,429 mg m −3 (0.14-17,489 µg C m −3 ) offshore and from 0 to 260 mg m −3 (0-1,872 µg C m −3 ) over the shelf. 1 -shelf stations with<300 m depth, 2 -offshore stations with>300 m depth. While doliolid developmental stages were not identified in 2016, gonozooids were the most abundant zooid stage in 2017 making up 93.5% of the population with the rest being nurse stages ranging from 1.5 to 20 mm in total length. Our relatively coarse nets probably severely undersampled smaller zooids making analysis of generation distribution patterns impossible.
The mean individual gonozooid wet mass at stations, where the bloom was developed in 2017, was 4.15 ± 0.64(SE) mg. This value corresponds to 6.6 mm total length calculated using the length-mass regression for D. tritonis gonozooids from the western subarctic Pacific (Nakamura et al., 2017). This size corresponds to 36.43 µg C according to the length-C mass relationship obtained for laboratory-reared D. gegenbauri gonozooids (Gibson and Paffenhofer, 2000). The GAM model with four predictor variables explained 58.5% of the deviance in doliolid abundance. Depth (F = 5.7, p < 0.001), salinity (F = 2.4, p < 0.001), and chl-a (F = 0.925, p < 0.01) had non-linear effects, while temperature did not have any significant effect on doliolid abundance (Figure 6). Abundance sharply increased with depths increasing to 1,000 m and it showed another increase at stations over 3,000 m depth. Salinity had positive effect on abundance when it had exceeded 31.25. The effect of chl-a on doliolid abundance was negative when chl-a concentration exceeded ∼0.7 mg m −3 .
The CPR surveys covered a two-decade period with contrasting "warm" vs "cold" stanzas in the North Pacific highly correlated to El-Nino events (Figure 7). The interannual changes in doliolid abundance in the GOA during 2000-2017 showed a distinctive pattern apparently related to climate variability (Figure 8). During strong positive PDO years, doliolids were substantially more abundant and they occurred in deep waters in both the eastern and western GOA. In contrast, doliolids were absent from CPR samples during negative PDO years. During neutral PDO years, doliolids appeared to be restricted to the southeastern GOA where they occurred in moderate quantities.
In addition, Tunicata [reported as "Ascidia sp. (tunicate)" [sic] and "Salpa sp. (pelagic salp)" [sic]] contributed up to 60% by weight to overall juvenile sablefish (Anoplopoma fimbria) stomach contents collected in the northern GOA during a strong 1997-1998 El-Nino, while during other years they were virtually absent from the diets (Sigler et al., 2001 ; Figure 9). Juvenile sablefish are pelagic predators and feed almost exclusively near the surface and in the upper layers (e.g., Brodeur and Pearcy, 1992). Therefore, it is plausible that the reported "Ascidia sp.", a distinctly sessile organisms morphologically somewhat resembling doliolids, were misidentified, and the report in fact refers to doliolids. Since the contribution of salps to total fish prey weight on average was<1% (Sigler et al., 2001), the increase of tunicates in the sablefish diets during El-Nino years was probably facilitated by doliolids.

DISCUSSION
For the first time in the southeastern GOA, we recorded dense (>100 ind m −3 ) doliolid aggregations during summer 2016-2017, typically classified as "blooms" in other world regions    (Crocker et al., 1991;Paffenhofer, 2013;Takahashi et al., 2015). The aggregations occurred at offshore stations influenced by the saline Alaska Current which connects the GOA to southern regions of the Pacific Ocean. Since the Alaska Current originates from the lower latitude (35-50 • N) North Pacific Current, it appears that the source population of D. tritonis is located within the northern reaches of the North Pacific subtropical gyre and the seed organisms are advected into the GOA. Observations from the US Northeast Pacific GLOBEC program (Weingartner et al., 2002) indicated that doliolids regularly occurred, albeit in extremely low numbers (<0.05 ind m −3 ), in the northern GOA on the outermost stations along the Seward Line since 1997 (Pinchuk, unpublished), which implies that the Alaska Current is the major conduit for doliolid dispersal along the Alaskan shelf. In contrast, few doliolids were found closer to the shore suggesting limited shoreward advection and/or unfavorable environment on the shelf affected by the low salinity, run-off driven Alaska Coastal Current.
The high abundance of doliolids offshore, where mean chl-a concentration was higher in both years, indicates strong dependence of the thaliacean bloom on phytoplankton availability. Doliolids can feed efficiently on a large variety of  2003and 2014, PDO neutral years were 2002, 2006, and 2010and PDO negative years were 2000, 2008, and 2011 prey from<5 µm bacteria to>100 µm diatoms (Crocker et al., 1991), allowing them to take advantage of intense diatom blooms typically occurring during spring and summer in the GOA. While large diatom spring blooms on the GOA shelf are triggered by the onset of thermal stratification starting in wind-protected fjord systems in March and extending over the shelf through April and May (Strom et al., 2001(Strom et al., , 2006, summer blooms are often associated with mesoscale eddies which propagate along the shelf edge and bring iron-rich coastal water offshore (Coyle et al., 2012(Coyle et al., , 2019Aguilar-Islas et al., 2016). The relatively slow eddy propagation speed (1.5-4 km day −1 ) and their extended longevity (6-24 months) (Okkonen et al., 2001(Okkonen et al., , 2003 provide persistent conditions for a doliolid population to develop. The physical interactions between propagating eddies and adjacent water masses in the northern GOA are similar to those observed along the North American Atlantic coast (Okkonen et al., 2003). Studies in the South Atlantic Bight indicated that doliolid blooms on the broad continental shelf during summer are associated with formation of mesoscale eddies spinning off the Gulf Stream at the shelf break, and their consequential propagation along the shelf (Deibel and Paffenhofer, 2009). However, in contrast to the GOA, the Gulf Stream eddies bring nutrients from the offshore onto the shallow shelf and foster high rates of phytoplankton production, resulting in favorable feeding environment for doliolids nearshore (Paffenhofer and Lee, 1987;Deibel and Paffenhofer, 2009).
Despite the doliolid's apparent preference to the offshore water rich in phytoplankton, the GAM showed negative effect of high chl-a concentration on the doliolid abundance. Similarly, the location of dense doliolid patches closely coincided with areas of low chl-a concentration in the Oyashio-Kuroshio region (Nakamura, 1998;Takahashi et al., 2015). This is not surprising considering the high filtration ability of doliolids. Laboratory studies showed that D. gegenbauri gonozooid clearance rates measured at phytoplankton concentrations ranging from 20 to 160 µg C L −1 were similar depending mainly on temperature and individual size, and that a gonozooid of 6.5 mm total length (=35 µg C weight) incubated at 16.5 • C cleared approximately 275 mL of seawater per day ingesting ∼15 µg C L −1 (Gibson and Paffenhofer, 2000) (see their Figures 5, 6). With the exception of one station in 2016, phytoplankton concentrations in our study area never exceeded 100 µg C L −1 , assuming C mg:Chl-a mg ratio of 51 (Lomas et al., 2012). Thus, assuming the same laboratory-derived constant feeding rate for D. tritonis and applying Q 10 = 3 (Haskell et al., 1999) to correct for the mean 14.25 • C observed in 2016-2017, the aggregated gonozooids at ∼3,000 ind m −3 can sweep clear their resident water volume in ∼47 h ingesting 0.045 g C m −3 day −1 . Similarly, less dense aggregations of doliolids in the Oyashio-Kuroshio transitional region cleared their residential volume in ∼10 days (Ishak et al., 2020). While no direct measurements of primary productivity were obtained during our cruises, model simulations run for 2000-2013 using the Regional Ocean Modeling System (ROMS) with an embedded nutrient-phytoplankton-zooplankton (NPZ) estimated late summer (September-October) mean primary production at 30.5 g C m −2 in the upper 30 m layer for the offshore region "P5" from the Canadian boundary to 147 • W (Coyle et al., 2019), which converts to a daily value of 0.0169 g C m −3 day −1 . Therefore, dense aggregations of doliolids can successfully control primary production in the southeast GOA.
At 20 • C, it would take ∼20 days for a newly released D. gegenbauri gonozooid to grow to that size (6.5 mm total length or 35 µg C weight) (Deibel, 1982). Using this rate and applying the same assumptions, the minimal longevity of the D. tritonis bloom in the GOA can be estimated at ∼35 days. The presence of old nurse stages in 2017 also indicate intensive asexual reproduction which must have been persisting for a considerable period. Thus, doliolid blooms in the Gulf of Alaska can substantially impact pelagic food web by quickly reducing FIGURE 9 | Interannual variability in diet composition (% by weight) of 0-age sable fish in the Gulf of Alaska (from Sigler et al., 2001, modified). local resources available for other planktonic grazers such as copepods and euphausiids.
Gelatinous zooplankton have been considered a trophic culde-sac (Dubischar et al., 2012) due to high water content and low nutritional value per unit of volume (Verity and Smetacek, 1996;Sommer et al., 2002). However, recent studies suggest that gelatinous zooplankton has been underestimated as food (Cardona et al., 2012;González Carman et al., 2014). For instance, generally low protein and lipid weight-specific content (∼5% each, e.g., Dubischar et al., 2006) in Thaliaceans may be compensated by their relatively large body size and high density (Arai, 1986;Larson, 1986), which also make them more visible to visual predators. In addition, Thaliacean passive mobility requires comparably little effort from pelagic predators to consume large quantities (Arai, 1988). Juvenile sablefish (Anoplopoma fimbria), a commercially important species, have been documented consuming Thaliaceans in large proportions during their first summer (Brodeur and Pearcy, 1992;Sigler et al., 2001;Strasburger et al., 2018). During our surveys in 2016, Thalieaceans on average comprised ∼40% of juvenile sablefish diets by weight after examined in a fresh state at sea, while at locations where the availability of other prey (juvenile rockfish) was minimal they contributed up to ∼90% by weight (Strasburger et al., 2018;Strasburger, unpublished data). In the mid-1980s, juvenile sablefish consumed Thaliaceans up to 76% by weight off the coastal Washington and Oregon (Brodeur and Pearcy, 1992). Both studies showed an increase in the proportion of Thaliacea when juvenile rockfish are not available (Brodeur and Pearcy, 1992;Strasburger et al., 2018). This suggests that while Thaliaceans may not be preferred food, they can sustain juvenile sablefish during transit to more optimal prey patches, or to areas with higher prey densities, thus enhancing fish survival.
During the mid to late 1990s, pelagic Tunicata [sic] (assumed doliolids) were mostly absent from juvenile sablefish diet from the GOA (Sigler et al., 2001), even though their contribution was likely underestimated as samples had been preserved and examined in the lab at a much later date, when many gelatinous prey items may no longer be recognizable (Arai, 1988). The noticeable shift from predominant euphausiids to doliolids in the sablefish diet occurred in 1997-1998 coinciding with a strong El-Nino-Southern Oscillation (ENSO) event. ENSO events in the tropical Pacific affect the upper ocean coastal circulation in the GOA via coastal Kelvin waves and atmospheric teleconnections. In the GOA, the mixed layer temperatures were substantially (∼1.5 • C) warmer, and salinity and nitrate concentrations were lower in 1998 compared to 1999 (Weingartner et al., 2002). These differences were likely facilitated by increased runoff and strong downwelling over the coastal North Pacific in 1997-1998 (Weingartner et al., 2005). ENSO events also tend to destabilize the Alaska Current enhancing the formation of the multiple mesoscale (∼200 km) eddies along the coast, which slowly propagate into the GOA where they can survive for more than 1 year (Melsom et al., 1999). Thus, oceanographic conditions in 1997-1998 in the GOA might have promoted doliolid blooms, which were reflected in the sablefish diets.
Our long-term CPR time series clearly indicate elevated abundance of doliolids during warm years identified with positive PDO and ONI indices. The relatively small CPR aperture limits sampling efficiency of larger organisms at low densities and precludes direct comparison of CPR estimates of doliolid abundance to those obtained with the plankton net tows. However, the CPR is a consistent sampling methodology and the consistent occurrence of doliolids in the samples during positive and neutral PDO years and the lack of records during negative PDO years must be indicative of conditions favorable for doliolid development. The doliolid blooms occurred mostly in the eastern GOA off the Canadian coast never extending north of 55 • N. The area is known for frequent eddy formation contributing nutrient enhanced shelf waters to the deeper ocean (Crawford, 2002;Whitney and Robert, 2002), which may have promoted doliolid blooms during positive PDO years.
The most prominent feature of the climate in the GOA is the Aleutian Low-Pressure System (AL), for which the average geographic location changes periodically during the winter (Mundy and Olsson, 2005). The position and strength of the Aleutian Low Pressure System (AL) affects oceanographic conditions in the GOA via differential climate patterns (Sugimoto and Hanawa, 2009). It has a profound effect on marine ecosystems via changes in sub-polar gyres (Ishi and Hanawa, 2005) and sea surface temperatures (Latif and Barnett, 1996;Mantua et al., 1997). It also has a strong link to larger climate patterns (Pacific North American Pattern, PNA) that drive temperature and precipitation over large spatial scales (Trenberth and Hurrell, 1994). During positive PDO years, the increased intensity of the wind-driven subarctic Alaska Gyre (Sugimoto and Hanawa, 2009) results in higher concentration of upwelled nutrients inside the gyre, while intensification of the Alaska Current enhances eddy formation along the GOA coast (Melsom et al., 1999;Ladd et al., 2007). During negative PDO years, a weaker low pressure facilitates reverse wind patterns resulting in relaxation of the Alaska Gyre and Alaska Current, in tandem with an increase of southerly flow into the California Current (Hollowed and Wooster, 1992), and a reduced eddy formation. With respect to doliolids, a stronger AL would provide ideal conditions for their bloom formation.
A combination of several factors resulted in the anomalous marine heat wave in the GOA in 2016. Unusually high atmospheric pressure reduced heat loss to the atmosphere in the northeastern Pacific south of the GOA and initiated the formation of the "warm blob" as early as in January 2014 (Bond et al., 2015). The warm temperatures were observed in the coastal GOA where they remained confined to the upper 100 m until the end of 2014 spreading into deep layers in 2015 and "preconditioning" the water (Walsh et al., 2018). A strong AL positioned over the Alaska Peninsula combined with a strong El Niño with a positive PDO (warm) phase led to enhanced atmospheric advection northward, apparently bringing extra warmth to Alaska in early 2016, while other factors potentially related to increased greenhouse emissions (e.g.,changing albedo from reduced snow cover over the terrestrial biomes) may have also contributed (Walsh et al., 2017). Climate model simulations indicate that warmth of this magnitude becomes the norm by the 2050s if greenhouse gas emissions follow their present scenario (Walsh et al., 2017). In that case, dense doliolid blooms in the GOA would become persistent and may profoundly change the GOA marine food web structure with potential repercussions to the ocean-based economy.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AP perceived the general concept and wrote the manuscript together with SB and WS. SB contributed CPR data and did retrospective analysis. WS and AP collected and analyzed the oceanographic and biological data in 2016-2017. All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
This publication is the result in part of research sponsored by the Cooperative Institute for Alaska Research with funds from the National Oceanic and Atmospheric Administration under the cooperative agreement NA13OAR4320056 with the University of Alaska. The NOAA Eastern GOA Ecosystem Assessment survey was completed aboard the chartered vessel Northwest Explorer, we thank the captain and crew for their expertise and in aiding a successful survey. Pacific CPR data collection is supported by a consortium for the North Pacific CPR survey coordinated by the North Pacific Marine Science Organization (PICES) and comprising the North Pacific Research Board (NPRB), Exxon Valdez Oil Spill Trustee Council through Gulf Watch Alaska, Canadian Department of Fisheries and Oceans (DFO) and the Marine Biological Association of the UK.