Disease epidemic and a marine heat wave are associated with the continental-scale collapse of a pivotal predator (Pycnopodia helianthoides).

Multihost infectious disease outbreaks have endangered wildlife, causing extinction of frogs and endemic birds, and widespread declines of bats, corals, and abalone. Since 2013, a sea star wasting disease has affected >20 sea star species from Mexico to Alaska. The common, predatory sunflower star (Pycnopodia helianthoides), shown to be highly susceptible to sea star wasting disease, has been extirpated across most of its range. Diver surveys conducted in shallow nearshore waters (n = 10,956; 2006-2017) from California to Alaska and deep offshore (55 to 1280 m) trawl surveys from California to Washington (n = 8968; 2004-2016) reveal 80 to 100% declines across a ~3000-km range. Furthermore, timing of peak declines in nearshore waters coincided with anomalously warm sea surface temperatures. The rapid, widespread decline of this pivotal subtidal predator threatens its persistence and may have large ecosystem-level consequences.

INTRODUCTION

Host-pathogen theory predicts that multihost pathogens can cause extreme population impacts, including extinction of susceptible species if they are continuously infected from reservoir species (, ). For example, introduced multihost pathogens such as avian malaria and avian pox have driven multiple native Hawaiian bird species to extinction (). Similarly, the Batrachochytrium dendrobatidis pandemic may have caused hundreds of species extinctions worldwide and decimated more than 38 amphibian species in Central America (, ). In another example, spillover of shared pathogens from domesticated bees to wild bumblebees (e.g., deformed wing virus and Nosema ceranae) is driving declines in the wild populations (). These pathogen-associated impacts are further exacerbated in a changing climate (–).

Since 2013, sea star wasting disease (SSWD) has caused massive, ongoing mortality from Mexico to Alaska (known as the Northeast Pacific SSWD event). In particular, the epidemic phase of the Northeast Pacific SSWD event (2013–2015) was notably different from previous events elsewhere in terms of its geographic extent, persistence, involvement of multiple species, symptoms in reproductive stars, and the extremely rapid progression of disease to death (–). More than 20 asteroid species have been affected in what is currently the largest documented epizootic of a noncommercial marine taxon (, , ). Diseased sea stars develop progressively worse dermal lesions (, ), arms detach from the central disc, gonads spilled from fully reproductive stars and individuals die, often leaving white piles of ossicles and disconnected limbs (fig. S1). Sea star mortalities during the first years of the Northeast Pacific SSWD event (2013–2015) were linked to a sea star–associated densovirus (SSaDV; family Parvoviridae), based on metagenomic analysis of bacteria and viruses in field samples, the experimental generation of disease in sea stars challenged with nonheated viral-sized material, the correlation between SSWD progression and SSaDV loads, and higher SSaDV prevalence in symptomatic stars (). Further support for densovirus involvement in the Northeast Pacific SSWD event includes experimental infection and morbidity in Pycnopodia helianthoides (), the relationship between asteroid densovirus load and SSWD in P. helianthoides, and experimental transmission of SSWD disease to asymptomatic individuals through exposure to a viral-sized agent (0.22 μm) from SSWD symptomatic stars (, ).

While the Northeast Pacific SSWD event caused severe reductions of the keystone intertidal ochre sea star across its entire west-coast range (, , –), the impacts on subtidal species are less well known. Because it was the most abundant subtidal star and also the most susceptible in early reports, initial studies quantified devastating declines of P. helianthoides, in Washington and British Columbia, linked with high prevalence of wasting disease (, ), but these observations are limited to the initial years, specific regions, and shallow depths accessible to divers. Hence, the current status of P. helianthoides in the Northeast Pacific and at all depths is unknown. This species is an important predator of sea urchins, and urchin populations released from top-down predatory control can expand and threaten kelp forests and biodiversity (, ). In many locations, P. helianthoides is the apex subtidal predator when other urchin predators are absent (, ). In these locations, its decline has triggered a trophic cascade, causing urchin populations to explode and kelp to rapidly diminish (). If P. helianthoides reductions occur throughout its range at the magnitude reported in local surveys, and at greater depths, this disease could not only threaten the long-term persistence of this species but also have wide-ranging cascading ecosystem effects.

Increasingly warm or anomalous temperatures are being shown to influence the prevalence and severity of marine infectious diseases [e.g., ()]. Experimental and field studies support a role for temperature in SSWD morbidity. Current evidence suggests that at warmer temperatures, Pisaster ochraceus have a higher risk of infection and progression to mortality (, , , ). If this relationship applies to other sea star species, then we expect that declines of P. helianthoides populations will be associated with warmer temperature exposure.

Here, we investigated the current status of P. helianthoides in both shallow nearshore and deep offshore waters, from California to British Columbia (~3000 km), using data from over a decade of complementary survey methods that cover from pre- to post-outbreak periods of the Northeast Pacific SSWD event. We report the rapid collapse of P. helianthoides populations along most of its range after the onset of the Northeast Pacific SSWD event, and at all depths, confirming the lack of a deep-water refuge for this species. Furthermore, to explore the hypothesis that warmer waters may be linked to the decline, we assessed the relationship between P. helianthoides abundance in shallow nearshore waters and sea surface temperature (SST). We detected a negative association between P. helianthoides abundance and anomalously warm SST. Potential ecosystem impacts of this decimation are discussed.

RESULTS

To assess P. helianthoides decline in deep offshore waters, we estimated the average yearly biomass (kg/10 ha) collected in 8968 bottom trawls (55 to 1280 m depth) conducted from California to Washington between 2004 and 2016. Deep-water trawl surveys showed fluctuating but constant P. helianthoides biomass up to 2011–2012 and an unprecedented decline after the onset of the Northeast Pacific SSWD event. In California and Oregon, the average biomass decreased 100% during 2013–2015 (from 2.78 and 1.73 kg/10 ha, respectively). In Washington, average biomass declined 99.2% (from 3.11 to 0.02 kg/10 ha) during this period. In 2016, no P. helianthoides were collected across the 1264-ha area covered by 692 trawl surveys. The collapse in biomass collection occurred 1 year earlier in California compared with other regions (2014; Fig. 1).

Fig. 1

Continental collapse of a pivotal predator: Deep offshore surveys.

Mean biomass of sunflower star in 8968 deep offshore trawls (55 to 1280 m) from (A) Washington, (B) Oregon, and (C) California from 2004 to 2016 with 95% confidence interval in light blue. Gray line marks the year 2013 for comparison of SSWD initiation across regions. Yellow circles depict the 2013–2016 trawl locations. The trawls per jurisdiction per year are shown in the top of each plot.

To quantify the impact of the Northeast Pacific SSWD event on P. helianthoides in shallow nearshore waters, we analyzed changes in abundance reported in 10,956 roving-diver surveys conducted between 2006 and 2017 from California to Alaska (sparse survey coverage for Alaska not plotted). Abundance categories (ACs; 0 to 4 corresponding to 0, 1, 2 to 10, 11 to 100, and >100 individuals, respectively) were analyzed within regional jurisdictions as an annual abundance score (). Furthermore, we report nearshore biomass of P. helianthoides (kg/10 m2) along the coast of central British Columbia using annual subtidal belt transect surveys conducted between 3 and 18 m depth in 2010–2011 and again between 2013 and 2017.

In the years before the Northeast Pacific SSWD event onset, ACs 2 (2 to 10 stars) and 3 (11 to 100 stars) were the most commonly reported (64 to 80% of the reports), but since 2014 (after the epidemic onset), at least 60% of surveys across the study area and up to 100% in California and Oregon report declines to ACs 0 and 1 (Fig. 2). Belt transect surveys in central British Columbia also showed that P. helianthoides biomass declined by ~96% (from 0.57 to 0.93 kg/10 m2 in 2010–2014, pre-Northeast Pacific SSWD event) to virtually zero (0.01 to 0.07 kg/10 m2) in 2015–2017. The abundance score from the shallow nearshore waters surveys revealed fluctuating but regionally constant P. helianthoides abundance during 2006–2013 and a consistent continental-wide decline after the onset of the Northeast Pacific SSWD event (Fig. 2).

Fig. 2

Continental collapse of a pivotal predator: Shallow nearshore surveys.

(A to D) Percentage of shallow nearshore ACs of sunflower star (P. helianthoides) reported in roving-diver surveys from southern California to southern British Columbia, Canada, from 2006 to 2017 (blue scale bars, right axis). Black line, annual abundance score (left axis); red line, annual mean of the maximum temperature anomaly 60 days before each survey (whiskers, 95% confidence interval; left axis). (A) British Columbia. (B) Washington. (C) Oregon. (D) California. (E) Mean biomass (kg/10 m2) in belt transect surveys in central British Columbia, with 95% confidence interval in light blue. Yellow circles depict the 2013–2017 locations. The red rectangle depicts the area where the belt transect surveys were conducted. The surveys per jurisdiction per year are shown in the top of each plot. For other details, see Fig. 1.

To assess the relationship between P. helianthoides declines at shallow nearshore waters and SST, we modeled the ACs reported for this species in the roving-diver surveys as a function of a satellite-derived SST anomaly metric and days since the SST metric was observed. Between 2013 and 2015, SST anomalies (departures from the climatologically expected temperature calculated from 1985 to 2012) were warm at all locations, but their magnitude and duration varied with latitude (Fig. 3). In California, 2014 was anomalously warm, increasing to extreme warming throughout 2015 with a peak anomaly of 4°C. In Washington, 2013 was anomalously warm during July to October but otherwise followed the seasonal climatology. In July 2014, a prolonged warm anomaly began, increasing to an extreme 2.5°C anomaly through 2015, coinciding with the long residence of the heat wave in the northeast Pacific Ocean (). Central British Columbia followed a similar pattern seen in Washington, with the arrival of the anomalously warm waters in the fall of 2015 (). Oregon was more variable than other locations, likely because of periodic cold upwelling events.

Fig. 3

Ocean temperature anomaly averaged over the roving-diver survey locations for the three initial years of the epidemic.

Blue, 2013; green, 2014; red, 2015. BC, British Columbia; WA, Washington; OR, Oregon; CA, California.

On the basis of ordinal regression models of P. helianthoides abundance and biologically relevant SST anomaly metrics, we provide evidence that P. helianthoides declines in shallow nearshore waters were associated with the maximum temperature anomaly exposure from within 60 days before each survey (tables S1 and S2). Our selected model indicates that with every 1°C increase in the maximum temperature anomaly, we would expect a 6% increase in the log odds of observing a low AC compared with all higher ACs (AC 0 versus 1 to 4, AC 1 versus 2 to 4, etc.), when all other variables are held constant (table S2). We evaluated the goodness of fit between the model and data using Nagelkerke’s pseudo R2 () and estimated that our model explained 68.5% of the variance in sea star ACs compared to a null model. To investigate the role of the maximum temperature anomaly exposure from within 60 days before the survey, we also calculated a pseudo R2 for a model that included all covariates other than this variable. This model explained 30.5% of the variance in P. helianthoides ACs, suggesting that this SST anomaly metric alone explains ~38% of the variance.

DISCUSSION

Using longitudinal data from complementary survey methods, we show (i) consistent or slightly increasing populations of P. helianthoides across most of its natural range in the decade before the Northeast Pacific SSWD event in shallow nearshore waters; (ii) consistent populations of P. helianthoides across most of its natural range in the decade before the Northeast Pacific SSWD event in deep offshore waters; (iii) a continental-scale collapse of this major ocean predator in these habitats associated temporally and spatially to the Northeast Pacific SSWD event and, therefore, to the multihost SSaDV; and (iv) an association between P. helianthoides declines in shallow nearshore waters and period of anomalously warm nearshore waters.

Our statistical analysis provides evidence for anomalous temperature as a key facilitator of the disease-related declines in the shallow nearshore waters, explaining more than a third of the variance by itself. These results align with field and experimental evidence of the interacting roles of temperature and SSWD in P. ochraceus morbidity and mortality, where infected stars exposed to warmer temperatures died at a faster rate (, , , ). Previous studies suggest that high water temperatures are associated with lower coelomic fluid volumes, higher metabolic demands, and metabolic stress in asteroids (, –), making the case for how a viral epidemic could be exacerbated in an invertebrate with a limited immune response capability (). Although warming waters likely accelerated and increased the scale of disease-induced morbidity, P. helianthoides mortality still occurred at high levels in colder temperatures of British Columbia. This finding supports a facilitating role of anomalously warm temperature for disease morbidity as previously reported in intertidal P. ochraceus during the Northeast Pacific SSWD event ().

Our data document the widespread decline of P. helianthoides at depths beyond shallow nearshore waters, confirming the lack of a deep-water refuge for this species. Available data support that the Northeast Pacific SSWD event is the most parsimonious explanation for this collapse. This species was identified as the most susceptible to SSWD (, ), the observed widespread declines occurred right after the onset of this event, and they followed the timing and spatial pattern of the declines observed in nearshore P. helianthoides populations and intertidal P. ochraceus populations (). Moreover, necrotic stars with autotomizing arms were observed by one of the coauthors (K. Bosley) in the trawls at the initiation of the Northeast Pacific SSWD event. The limitations of deep offshore water temperature data prevented an analysis of the role of water temperature in exacerbating P. helianthoides declines at deep offshore waters, as this information is collected once annually. However, new findings show that anomalous warm temperatures associated with recent marine heat waves were broadly detected at about 140 m depth in the Northeast Pacific beginning in mid-2014 (), which may have contributed to SSWD morbidity at these depths.

The abrupt decline of P. helianthoides after the onset of the Northeast Pacific SSWD event occurred in shallow and deep waters regardless of the abundance before the collapse. Similarly, Miner et al. () did not detect a relationship between the degree of population decline and pre-outbreak intertidal P. ochraceus density. The current population status of P. helianthoides raises important questions about the potential for recovery and persistence of this species. With SSWD, reservoir species could be a continuous source of SSaDV to remaining P. helianthoides given that asymptomatic star species within the range of P. helianthoides have tested positive for SSaDV genetic material (). In its current status and with new bouts of mortality recorded by divers in August 2018, continuous infection from reservoir populations or small stochastic disturbances could cause the restricted remnant populations of P. helianthoides to vanish (, , , ).

Cascading effects of the P. helianthoides loss are expected across its range and will likely change the shallow water seascape in some locations and threaten biodiversity through the indirect loss of kelp (, , , ). P. helianthoides was the highest biomass subtidal asteroid across most of its range before the Northeast Pacific SSWD event (). Loss or absence of this major predator has already been associated with elevated densities of green (Strongylocentrotus droebachiensis), red (Mesocentrotus franciscanus), and purple urchins (Strongylocentrotus purpuratus) across their range (, , , ), even in regions with multiple urchin predators (). Associated kelp reductions have been reported following the outbreak (, ). Examples from other widespread marine diseases—the near extirpation of the intertidal sea star Heliaster kubiniji from Gulf of California (), the mass mortality of the urchin Diadema antillarum from Caribbean reefs (), and the withering syndrome–influenced endangerment of multiple California abalone species ()—demonstrate how the loss of key species can drive community effects that influence marine ecosystem processes.

SSWD, the anomalously warm water, P. helianthoides declines, and subsequent urchin explosions (fig. S2) have been described as the “perfect storm.” This “storm” could result not only in trophic cascades and reduced kelp beds () but also in abalone and urchin starvation (). We encourage scientific review of P. helianthoides recovery, assessment of the probability for endangerment, and, more generally, expanded surveillance of the consequences stemming from complex interactions that emerging infectious diseases and warming ocean temperatures can have in shaping ecosystems, even in the deep ocean.

MATERIALS AND METHODS

Study area

This study includes data from across the northeast Pacific, encompassing most of the natural range of the endemic sunflower sea star P. helianthoides. We analyzed data from shallow nearshore roving-diver surveys collected in California, Oregon, and Washington (USA) and British Columbia (Canada), in addition to deep offshore trawl surveys conducted from California to Washington.

Deep offshore surveys

The National Marine Fisheries Service [National Oceanic and Atmospheric Administration (NOAA)] conducted 8968 individual bottom trawls along the coasts of California, Oregon, and Washington between 2004 and 2016. These bottom trawls quantify the biomass of P. helianthoides (kg) as well as the total area swept by the trawl (ha). We estimated the mean kg/10 ha per state per year using bottom trawls conducted between 55 and 1280 m depth. The 95% confidence intervals of the yearly means were calculated.

Shallow nearshore surveys

Trained and tested recreational scuba divers searching for P. helianthoides on the coast of Washington, Oregon, and California (USA) and northern British Columbia (Canada) conducted 10,956 roving-diver surveys between 2006 and 2017. After the surveys were completed, the abundance of P. helianthoides was ranked between 1 and 4 corresponding to estimated ACs of 1, 2 to 10, 11 to 100, or >100 individuals sighted, respectively. This information was submitted to the Reef Environmental Education Foundation (REEF) Volunteer Fish Survey Project Database (). With this information, we calculated P. helianthoides abundance score as previously explained (, ) per state and per year. Those surveys that did not report P. helianthoides presence were assigned 0 abundance (category 0). Furthermore, we calculated the proportion that each AC, including the category 0, was reported per state and per year. The surveys were conducted in kelp forest, rock/shale reefs, open ocean, sea grass beds (Phyllospadix and Zostera spp.), pinnacles, bull kelp beds (Nereocystis sp.), cobblestone/boulder fields, and walls over 10 feet high.

The abundance and size of P. helianthoides were recorded on annual scuba surveys using belt transects at 11 rocky reef sites located on the central coast of British Columbia between 2010 and 2017 (3 to 15 m depth). Belt transects (30 m × 2 m, n = 6 per site) were conducted at all 11 sites in 2013–2017, whereas 8 and 4 sites were surveyed in 2011 and 2010, respectively (belt transects, 10 m × 2 m; n = 3 to 9 per site). P. helianthoides biomass was calculated using a length-to-biomass regression () and summed across transects at a site to yield kg/10 m2. P. helianthoides biomass for the central coast region (range from 51°24.612′N to 52°4.242′N) was calculated as the mean across all site-year combinations. The 95% confidence intervals of these yearly means were calculated.

Sea surface temperature

Satellite SST data were obtained at 0.05° (~5 km) daily resolution from the CoralTemp product by NOAA Coral Reef Watch (NOAA Coral Reef Watch 2018). Time series of SST were acquired at each location of the roving-diver surveys or at the nearest neighboring satellite pixel near the coastal boundary. SST anomalies describe the variation in temperature from expected values at a given time of year and location and were determined using the Coral Reef Watch approach based on monthly climatologies (). Where SST is warmer (cooler) than expected, the SST anomaly is positive (negative). For each survey, maximum values of SST and SST anomaly from several periods immediately prior were extracted (30, 60, 90, 180, and 360 days) for comparison with survey data. Jurisdictional (California, Oregon, Washington, and British Columbia) SST anomaly summaries for each year were calculated by spatially averaging 60-day-prior maximum values from survey locations for that year (Fig. 2). Jurisdictional SST and SST anomaly time series (Fig. 3 and fig. S3) were averaged across all survey locations.

Statistical analysis

We estimated the relationship between sea star abundance reported in the shallow nearshore roving-diver surveys and SST by fitting a hierarchical ordinal regression model with a probit link functionwhere the cumulative probability of the ith survey falling in the jth AC or below is modeled as a function of the following: θ threshold parameters across ACs, which provide a separate intercept for each category j, an SST anomaly metric (SSTmetric), days since the SST anomaly metric was observed (Days.SSTmetric), year (Year2007–2017), month (Month), and latitude (Latitude). Year was included as a fixed effect to determine the year when sea star abundances collapsed (years that were statistically significantly different from the 2006 baseline). Month and latitude were included as random effects to account for additional variation over time and space. σ2μ and σ2γ corresponded to the variance of the distribution of month and latitude random effects, respectively. Our data met all model assumptions: (i) the response variable was measured on an ordinal scale; (ii) the predictor variables were continuous or categorical; (iii) there was no multicollinearity among predictor variables, which we assessed with correlation tests for correlations between two predictors and visually for correlations among three predictors; and (iv) there were proportional odds between each AC as indicated by nearly identical effects among generalized logistic regression models comparing each AC split individually (slopes < 2). We fit 10 candidate models that included the year, latitude, and month covariates and one of the following SST metrics: the maximum SST in the 30, 60, 90, 180, or 360 days prior to each roving diver survey; or the maximum anomalous SST in the 30, 60, 90, 180, or 360 days prior to each roving-diver survey. We compared the AIC value of the candidate models with and without the covariate “days since the SST metric was observed,” and then selected the model with the lowest AIC value (tables S1 and S2). We assessed convergence of models by inspecting the maximum absolute gradient of the log-likelihood function and the magnitude of the Hessian. Each model was empirically identifiable by ensuring that the condition number of the Hessian measure was no larger than 104 (). We evaluated variance explained by the final model using Nagelkerke’s pseudo R2 (). Nagelkerke’s pseudo R2 is a commonly used statistic to measure goodness of fit that is calculated by comparing likelihood ratios between a full model and an intercept model. We conducted this analysis in R statistical software v3.4.3 () using the clmm function of the “ordinal” package () for the ordinal regression model and the nagelkerke function in the “rcompanion” package to calculate pseudo R2 values ().