Marine Heatwaves Modulate Food Webs and Carbon Transport Processes

doi:10.21203/rs.3.rs-5256564/v1

Download PDF

Article

Marine Heatwaves Modulate Food Webs and Carbon Transport Processes

https://doi.org/10.21203/rs.3.rs-5256564/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

The impact of marine heatwaves (MHWs) on ecosystem functions and services remains poorly constrained due in part to the limited availability of time-resolved datasets integrating physical, chemical and biological parameters at relevant operating scales. Here we paired over a decade of autonomous Biogeochemical (BGC)-Argo profiling float data with water column plankton community composition profiles to investigate the impacts of MHWs on carbon cycling dynamics in the Northeastern subarctic Pacific Ocean (NESAP), resolving temperature-induced changes in production and transport of particulate organic carbon (POC) within the water column. POC concentrations were extremely high during MHWs in Spring 2015 and 2019, associated with detritus enhancement and changes in plankton community structure. Instead of sinking, particles < 100 µm accumulated in mesopelagic waters where remineralization happened slowly over the year, thus decreasing the rapid deep particle export and carbon sequestration potential of thermally-impacted waters.

Earth and environmental sciences/Biogeochemistry

Earth and environmental sciences/Ecology/Molecular ecology

Earth and environmental sciences/Environmental sciences/Environmental impact

Earth and environmental sciences/Ecology/Climate-change ecology

Earth and environmental sciences/Ocean sciences/Marine chemistry

biological carbon pump

BGC-Argo profiling floats

carbon cycle

carbon export

climate change

vertical migration

DNA metabarcoding

marine heatwaves

marine microbes

particulate organic carbon

Marine heatwaves manifest as prolonged warming events persisting from weeks to years. The resulting thermal stress can lead to increased mortality and biodiversity loss with concomitant changes in food web structure and biogeochemical cycling^1,‎2. Such impacts can overcome ecological resilience thresholds resulting in state changes that negatively impact ecosystem functions and services³. Ocean observations and models suggest that MHWs have been expanding in size and intensifying over the past few decades as a function of climate change⁴. Even a reduction in anthropogenic carbon emissions proposed by current climate action goals is unlikely to decelerate spatial expansion and duration of these events for many years to come⁴.

From the standpoint of carbon cycling dynamics, elevated ocean temperature promotes water column stratification and decreases vertical mixing with a resulting impact on surface nutrient renewal and primary production⁵. These changes percolate through marine food webs with the potential to alter the size and composition of particulate organic carbon (POC) and the efficiency of the biological carbon pump (BCP), an integral global mechanism for carbon export and deep ocean sequestration^6,‎7,8. Despite the magnitude of these effects, integrated analysis of plankton communities at the base of marine food webs and POC production and export during MHWs is limiting⁹.

Two successive and persistent MHWs impacted NESAP waters in 2013-15^{‎5,10‎,11} and 2019-20^‎11,12. The first event called “The Blob” was catalyzed by suppressed wind stress during the winters of 2013-14. Air-sea teleconnections with the equatorial Pacific during the development of an extreme El-Niño added persistence to the warm event and increased upper ocean stratification^10‎,11. Water column impacts included a decline in surface dissolved iron concentrations^‎13, and changes in primary productivity^‎5,‎14 and microbial community structure including heterotrophs‎¹⁵ and phytoplankton^‎16. The second MHW formed due to the weakening of the North Pacific high-pressure system and imprinted a low-density, low-salinity feature in the euphotic zone^‎11‎,12. High productivity during the summer of 2019 was associated with a volcanic eruption and land fires^{14,‎17,‎18}. Biogenic CaCO₃ export was insignificant in comparison to non-MHW years‎¹⁹, indicating potential changes in phytoplankton community composition and the BCP. During both events higher trophic levels were also affected, leading to the collapse of Bering Sea snow crab populations^‎20 and low birth rates among whales^‎21.

Here, we combine high-resolution BGC-Argo profiling float data collected along the Line P transect in conjunction with the NASA EXPORTS program²²between 2010-22 with pigment-based measurements of phytoplankton functional groups and DNA metabarcoding of plankton community structure including microscopic prokaryotic (bacteria and archaea) and eukaryotic constituents to explore the impacts of MHWs on food webs and carbon transport processes in the NESAP. BGC-Argo profiling floats equipped with temperature, salinity, nitrate, optical backscatter (b_bp, a proxy for POC) and fluorescence (a proxy for chlorophyll a) sensors profiled the area every 5 to 10 days, allowing high-resolution observations of net community production, particle production and export. Pigments and metabarcoding data collected near the float’s location provided biological context to link anomalous POC accumulation under thermal stress to plankton community structure. During both MHWs, the vertical transfer of small particles to a preferential mesopelagic depth indicated an intensification of zooplankton vertical migration, a well-known key export driver in the region^23,24. The subsequent slow remineralization of non-sinking particles <100 µm over the following months decreased the prospects for deep carbon export.

Both MHWs imprinted a warm and stratified water column effect (Fig. 1). Indeed, stratification at the base of the mixed layer was strongest during MHW years (2014-16 and 2019, Fig.1a). Surface temperature anomalies >1˚C during the first MHW persisted between 2014-15, agreeing with previous reports^‎‎5,10. Elevated temperatures propagated to subsurface waters over subsequent years with temperature anomalies up to +0.2˚C observed down to 300 m. During the second MHW, similar anomalies were observed between 2018-19 and propagated to deeper ocean layers down to 500 m persisting until 2021.

BGC-Argo profiling float data collected in the vicinity of stations P20-P26 along the Line P transect (Fig. 2a) revealed peaks of small particles <100 µm, maxima in chlorophyll (Chl) fluorescence concentrations and net community production (NCP) during spring and late summer over the time series, signaling enhanced primary productivity (Figs. S1-S4). Upper ocean stocks of the small POC fraction <100 µm increased during the productive season and were particularly elevated during the second MHW in the Spring of 2019 and Summer of 2020 (Fig. 2d,e). Baseline POC concentrations (i.e. POC during wintertime minima in January) were remarkably different between the two periods: the years preceding and during the first MHW (2010-15) had lower POC baselines in comparison to post-MHW (2018-22). For the first period, the minimum POC ranged between 248 mg m^-2 in 2011 and 395 mg m^-2 in 2012. For the second period, the POC baseline ranged between 596 mg m^-2in 2020 and 949 mg m^-2 in 2021.Fig. 2: Mean monthly chlorophyll to carbon (Chl:C) ratios and POC stocks (small particles) over the float’s time series. Middle panels (b-f) show years before and during the first MHW, and the right panels (c-g) show years after, including the second MHW. (a) Map of the study area along the Line P transect and surrounding NESAP waters indicating BGC-Argo floats tracking (colored dots), the approximate float position at the time of each MHW in 2015 and 2019 (black asterisks), and the Line P stations in which samples were taken for HPLC and genetics (magenta triangles). (b, c) Mean monthly Chl:C from the integrated upper 100 m depth. Red dashed lines denote MHW years. (d, e) Integrated mean monthly POC for each year within the upper 100 m depth. (f, g) Integrated mean monthly POC for each year between 100-300 m depth, denoting the upper mesopelagic.

POC accumulation during each year’s productive season was estimated in relation to the baseline as ∆POC = POC_max- POC_min. For the first half of the time series, upper-ocean ∆POC had similar values regardless of the year corresponded to a MHW or non-MHW year (Fig. 2d,e): ∆POC was 254 mg m^-2 in 2015 (MHW) versus 252 mg m^-2in 2011 (non-MHW). During the second half of the time series, POC accumulation was anomalously high in the spring of 2019 and different from any other year on record. ∆POC increased to 2041 mg m^-2 by April, a three-fold increase from the baseline and 8x higher than in 2015. This accumulation was not proportionally tracked by chlorophyll, resulting in a significant drop in Chl:C ratios (Fig. 2c,e). The late 2019 summer bloom resulted in a small peak in ∆POC, also with a drop in Chl:C ratios. In 2020, ∆POC was the highest for a late summer bloom at 1016 mg m^-2 by October. Similar to 2019, ∆POC was not tracked by chlorophyll, and Chl:C dropped significantly.

The productive seasons were followed by a subsurface increase in POC stocks within the upper mesopelagic layer between 100-300 m (Fig. 2f,g). This was a seasonal trend for the entire time series, resulting in higher POC concentrations at 200 m in relation to 100 m independent of the magnitude of the mixed layer shoaling during spring months (Fig. S5). However, MHW years experienced much higher accumulation during the spring of 2015 and 2019 (Fig. 2f,g). In 2015, there was a 3.2-fold change in mesopelagic POC stocks as ∆POC = 174 mg m^-2. The mesopelagic accumulation did not follow upper-ocean accumulation. In 2019, mesopelagic ∆POC was 1155 mg m^-2, marking a 3.7-fold increase from baseline values and the highest difference between deeper and shallower POC concentrations (Fig. S5b). This mesopelagic enhancement tracked a high upper ocean ∆POC (Fig. 2). During 2015 and 2019, there was a clear particle flux imbalance in the upper mesopelagic in comparison to non-MHW years during spring, when the system experienced lower ∆POC changes between 0.4- and 2.1-fold. The high ∆POC feature was also observed during the summer of 2020 (∆POC = 495 mg m^-2), although smaller in comparison to the spring peak. The accumulated POC during those MHW years slowly decreased over the following months, and concentrations dropped down to near baseline values by winter. This accumulation pattern at a preferential depth does not correspond to particle export by gravitational sinking, when an exponential decrease of particle concentrations over depth follows a Martin curve‎⁶ due to active particle export coupled with microbial respiration.

The timing between POC accumulation in the mesopelagic and its subsequent removal during MHWs is imprinted as concentration anomalies in relation to non-MHW years: overall, stronger anomalies were observed in 2015 and 2019 (Fig. 3) for the upper ocean and the mesopelagic, especially around 200 m. The main source of particles to the deeper layer can be tracked to spring months, with a small contribution of summer months in 2015. In 2014 and 2020, positive anomalies were smaller, shorter in duration and driven by summer particles.

To link the anomalous POC accumulation with phytoplankton production as a potential source of carbon, we computed NCP using the mass balance model based on BGC-Argo profiling float nitrate data (Fig. S4)¹⁹. In 2015, the largest fraction of NCP happened during the spring bloom at about 2.8 mol C m^-2 yr^-1, the highest on record for the season. Elevated NCP did not result in elevated upper ocean ∆POC, but was closely followed by mesopelagic accumulation (Fig. 2d,f). In 2019, spring blooms resulted in NCP around 1.2 µmol C m^-2 yr^-1, modest and comparable in magnitude to non-MHW years. There was a mismatch between modest spring NCP and high ∆POC (Fig. 2c,e), thus, POC accumulation with low Chl:C could not be attributed to increasing phytoplankton productivity and biomass. Summer NCP elevated annual production to about 3.8 C m^-2 yr^-1, but did not contribute significantly to the elevated ∆POC observed in the mesopelagic (Figs. 2g, 3d). Finally, the 2014 and 2020 summer blooms drove most of the NCP in those years at 2.4 C m^-2 yr^-1 and 0.4 C m^-2 yr^-1, respectively (Fig. S4). The 2020 low NCP did not track the elevated ∆POC that year (Fig. 2e,g).

Phytoplankton pigment data (chlorophylls and carotenoids) from stations P20-P26 along the Line P transect were used to estimate seasonal phytoplankton composition (Fig. 4). Total chlorophyll a concentrations were generally low (mean = 0.395±0.187 mg m^-3), but were elevated in the spring of 2019 as compared to other years (0.829 mg m^-3 vs. averaged 0.371 mg m^-3, respectively). Phytoplankton communities were dominated by haptophytes, chlorophytes and pelagophytes. Spring communities had additional contributions from cyanobacteria and diatoms, while summer had contributions from diatoms and dinoflagellates (Fig. 4a,b). Increased chlorophyte abundances were observed during both MHWs (Figs. 4c, S6), while higher diatom abundances were observed only in 2019, tracked by enhanced NCP (Fig. S4) and anomalous ∆POC (Fig. 2e,g).

Prokaryotic and eukaryotic plankton community composition was determined using small subunit ribosomal RNA (16S and 18S rRNA) metabarcode sequencing on samples collected in the surface mixed layer (10-100 m) and upper mesopelagic waters (100-300 m) at stations P20 and P26, flanking the BGC-Argo float paths (Fig. 2) between 2010 to 2019 (2015 to 2019 only for eukaryotes). As expected, community structure was correlated with depth (Figs. S7-S9), which explains 70% of the variation in prokaryotic communities and 24% of the variation in eukaryotic communities (p = 0.001). Season only explained <10% of the community variation (p = 0.001), perhaps not surprising given the weak seasonal harmonic in biochemical parameters at P20 and P26²⁵.

Across the time series, prokaryotic surface communities were dominated by cosmopolitan taxa like SAR11 within the Alphaproteobacteria, SAR86 and OM43 within the Gammaproteobacteria, Formosa within the Bacteroidia, and the cyanobacteria genusSynechococcus. In mesopelagic waters chemoautotrophic taxa were more abundant including Nitrosopumilus and Nitrosopelagicus, SAR324, Marinimicrobia, SUP05, and Nitrospina involved in coupled carbon, nitrogen, and sulfur cycling processes (Fig. S8). Dominant eukaryotic taxa included cyclopoid copepods (especially Oithona) and dinoflagellates distributed throughout the water column while other eukaryotes exhibited depth-specific trends (Fig. S9). For example, Ciliates affiliated with Spirotrichea and phototrophic taxa including the haptophytes Phaeocystis and Chrysochromulina, chlorophyte Bathycoccus, diatoms Bacillariophyceae and Mediophyceae, and Pelagophytes were more abundant in surface waters while radiolarians including Polycystinea, Acantharea, and RAD-B, and siphonophores were more abundant in upper mesopelagic waters. Dinoflagellate abundance shifted from Gymnodiniales in euphotic waters, to parasitic Syndinians in the upper mesopelagic (Fig. S9) - such heterogeneity was similar to observed in the Pacific^25,26,27. Calanoid copepod sequences, especially the genus Paracalanus in the upper ocean, made up a larger fraction of copepod sequences in both heatwave years (2015 and 2019, Fig. S9) compared to non-heatwave years.

The planktonic community inhabiting surface waters was relatively consistent between MHWs (Table S1) with only 8 prokaryotic and 6 eukaryotic taxa identified as significantly differentially abundant (Fig. 5). These included increased abundance of prokaryotic taxa Marinimicrobia, Pseudoalteromonas, Pirellulaceae and Cellvibrionaceae along with the eukaryotic diatoms Thalassiosiraceae and Fragilariopsis, MAST-3 stramenopile, and Syndiniales DG-I-Clade-5 in 2015 as well as various Gammaproteobacteria affiliated with Pseudomonas and Thiotricales and eukaryotic haptophytes in 2019. The magnitude of changes in mesopelagic community structure were more intense (Table S1), with numerous taxa coinciding with POC accumulation (Fig. 2f,g) identified as significantly differentially abundant (Fig. 5). For example, Synechococcus, SAR11 Clade I, SAR116, OCS116, SAR86, and Pseudohongiella were all more abundant in 2015 than in 2019 while HOC36 and UBA10353 were more abundant in 2019 (Fig. 5). Decreased SAR11 abundance in 2019 corresponded to a proportional increase in Nitrosopumilus, Nitrosopelagicus (though these were not identified as significant) and SAR324. Within the eukaryotic community, parasitic Syndiniales and heterotrophic dinoflagellates including Prorocentrum, Gyrodinium, and Gymnodinium taxa were identified as significantly differentially abundant in 2019. Increased abundance of these taxa corresponded to a decrease in Cyclopoid copepod abundance, typically the most abundant copepods in mesopelagic waters, and an increase in the Calanoid copepod Microcalanus (Fig. S9). Polycystinea and Acantharea radiolarians, especially Cladococcus and Acantharea F3, were significantly more abundant in 2015, while radiolarian groups RAD-A and RAD-B were elevated in 2019.

Integrated physical, chemical and biological parameter information collected over more than a decade of time-resolved observations shed light on marine food web structure, POC composition, and carbon export processes in NESAP waters during non-MHW years. A tightly coupled food web was observed in euphotic waters, with microzooplankton as the main consumers of heterotrophic bacteria and phytoplankton^23,24,27,28. Elevated POC during the productive season was composed of detritus and fecal pellets of relatively small sizes ranging between 1–6 µm^7,24,28 with the main drivers of POC export to mesopelagic waters linked to consumption of microzooplankton and picoeukaryotes^8,23,30. Processed at different trophic levels in the euphotic zone, highly-worked particles < 100 µm are typically the main carbon source for mesopelagic zooplankton dominated by heterotrophic and parasitic protists^23,24,33 (Fig. S9).

In response to thermal stress during successive and persistent MHWs, the NESAP experienced changes in plankton community structure and POC production and export, as summarized in Table 1. At the start of the first MHW in 2014, production was nearly balanced by export with no anomalous ∆POC accumulation observed (Figs. 2,4), similar to historical observations in the region^24,28. As thermal stress persisted into 2015, POC export became imbalanced and resulted in mesopelagic particle accumulation (Fig. 3b) likely caused by changes in zooplankton community structure and vertical migration patterns between sunlit and dark ocean waters (Figs. S4, 4, 5). For instance, the Radiolarians Acantharea and Cladococcus were elevated in 2015. Acantharea have been shown to export copious amounts of carbon from the euphotic zone as they migrate through the water column^29,30. Moreover, increased abundance of parasitic Syndiniales in mesopelagic waters may have further altered carbon processing³¹. Syndiniales have been implicated in particulate carbon flux attenuation³² and are associated with RAD-B radiolarians^33,34. Proliferation of Cladococcus^35,36 and the scavenging Copepod Oithonain^37,38 in 2015 (Fig. 1) could have further contributed to the production of small non-sinking fecal pellets in mesopelagic waters, decreasing the prospects for deep particle export

The period following the 2015 MHW was marked by elevated background POC (Fig. 2e) in an increasingly stratified water column (Fig. 1a) promoting further particle accumulation. In 2019 extreme POC accumulation was not proportionally tracked by NCP (Fig. S4) and chlorophyll (Figs. S2,4b), leading to low Chl:C ratios. This anomaly indicates an intensification of food web carbon processing resulting in excess detritus, which is already known to comprise most background POC in the NESAP²⁸. Interestingly, elevated NCP (Fig. S4) and high-chlorophyll levels in 2020 (Fig. 4a) promoted by wildfires and volcaniceruptions^{14,‎17,‎18} did not result in mesopelagic POC accumulation (Fig. 3c), suggesting that phytoplankton production alone cannot account for anomalous ∆POC accumulation. Rather, thermal stress induced by MHWs modulated food web structure causing changes in carbon processing between trophic levels.

Table 1

Cascading effects of MHWs in NESAP waters during spring and summer: changes in POC accumulation, NCP and plankton community structure. Parameters with slash symbols were comparable in magnitude with non-MHW years and not included here. The years of 2015 and 2019, when we saw the biggest changes in the ecosystem, are highlighted in darker red. *Community structure information was only available for the Summer of 2020.
MHW year	Season	∆POC surf	∆POC mesopel	∆NCP	∆Phytoplankton	∆Zoo/protozooplankton
2014	Spring	////	///	///	///	///
2014	Summer	///	///	Moderate	///	///
2015	Spring	///	High	High	○ Chlorophytes ○ Haptophytes ○ Diatom (Thalassiosiraceae, Fragilariopsis)	○ Radiolarians Acantharea and Cladococcus ○ Copepod Oithona
2015	Summer	///	High	///
2019	Spring	Extreme	Extreme	///	○ Diatoms (Chaetoceros) ○ Haptophytes	○ Radiolarians RAD-groups ○ Parasitic alveolates (Syndiales) ○ Pico-protozooplankton ○ Copepod Microcalanus
2019	Summer	///	///	High	○ Diatoms (Chaetoceros) ○ Haptophytes
2020	Spring	///	///	///	///*	No data
2020	Summer	Elevated	Elevated	///	///*	No data

The accumulated particles may have expanded the niche space for chemoautotrophic bacteria and archaea driving coupled carbon, nitrogen, and sulfur cycling processes across microscale oxyclines, with implications for nitrogen loss and climate-active trace gas production³⁹. The decline in cosmopolitan SAR11 clades in the 2019 MHW is noteworthy and could be linked to the heavy recycling and accumulation of organic matter. In NESAP, SAR11 (and SAR86 in 2015) was more abundant when DOM was less degraded⁴⁰. SAR11 and SAR116 (both more abundant in 2015) were more abundant during Phaeocystis (haptophyte) blooms in the Southern Ocean as compared to diatom blooms⁴¹. Thus, the 2019 MHW anomalous diatom bloom may have limited the success of SAR11.

Looking forward, more granular insight into the network topology of marine food webs is needed to resolve functional and activity-dependent changes in carbon processing pathways and the BCP in response to thermal stress and their associated impacts on deep particle export and carbon sequestration. This information can inform future climate models and provide key metrics for ocean health. It is therefore imperative to expand the time-series monitoring programs that provide essential parameter information and samples. At the same time, new methods and blue technologies are needed that link sample collection more closely with BGC-Argo profiling float measurements to more accurately constrain spatiotemporal changes in food web structure and carbon transport processes.

Temperature anomalies

The depth-resolved physical parameters shown in Fig. 1 were obtained from an Argo-based gridded product⁴². We selected data from Ocean Station Papa (50°N, 145°W) to analyze temporal variability over the past decade. Temperature anomalies at each depth were calculated as deviations from the climatological mean during the period 2004–2022, corresponding to the period matching BGC-Argo float and shipboard data. The mixed layer depth (MLD) was defined as the depth where the temperature is 0.2°C lower than that at 10 meters⁴³, which generally corresponds well with the depth at which the maximum Brunt-Väisälä frequency occurs. The gridded product enabled Brunt-Väisälä estimates between 2010 and 2021.

BGC-Argo float data analysis: POC and chlorophyll estimates

Data from the three floats used in this study were freely acquired from the Argo Global Data Assembly Center (USGODAE; usgodae.org/ftp/outgoing/argo/). Only quality-controlled files (Sprof) and data flagged as “good” were considered. Variables used in this study include temperature, salinity, nitrate and bio-optics (i.e., fluorescence chlorophyll and backscatter). We focused the analysis in the upper ocean between the surface and 400 m, with data every 5 to 10 days and acquired between June 2010 and October 2022, with a gap between February 2016 and August 2018. The study region was constrained between 47˚-52.5˚N and 137.5–146˚W (Fig. 2a). Details on quality control and sensor performances have been published^44,45. As previously described‎¹⁵, float-estimated chlorophyll a in our study region is biased when compared to nearby discrete chlorophyll a measurement from repeat Line P cruises, requiring additional post-correction. We applied a correction factor to the adjusted chlorophyll a data for each float, estimated by taking the median of the ratio of float-adjusted chlorophyll a data to the Line P corrected float data¹⁴ in the upper 15 m. The following correction factors were divided by our adjusted chlorophyll data to obtain corrected values extending to 2022.

Backscatter (b_bp) data was separated into small and large particles‎^46,47 (smaller or bigger than 100 µm, respectively), but only small particles were used in this study since large particles were a minor fraction of total b_bp and showed no temporal variability (Figs. S1-S3). To estimate POC stocks (Fig. 2) at the different ocean layers, we first estimated the euphotic zone depth (Ez) assuming it comprises the primary production zone^‎6,47 as:

Ez (m) = Z of 0.1* Chl_Max below Z _{chl_Max} (1)

Where Chl_Max is the depth of maxima chlorophyll for each float profile. Mean euphotic zone depth was 98 ± 20 m for the entire time series, so we assumed Ez = 100 m for the subsequent calculations. We then estimated mean monthly particle stocks between surface-Ez and Ez-300 m for each year and converted b_bp into total POC‎⁴⁴.

POC (mg/m ³ ) = 31,200 × bbp + 3.04 (2)

Chl:C ratios were estimated using mean monthly integrated POC and chlorophyll concentrations within the upper 100 m of depth.

POC anomalies (Fig. 3) derived from BGC-Argo floats were estimated by subtracting POC time series from each MHW year (2014, 2015, 2019 and 2020) from mean monthly anomalies during non-MHW years in a multi-step process. Since POC background signals were significatively different after the first MHW, we compiled non-MHW years into two datasets: the first was composed of float profiles from complete calendar years prior to 2013, and the second was composed of profiles from the years of 2018, 2021 and 2022. For each non-MHW dataset, we computed mean monthly [POC] and then linearly smoothed the data every 20m (floats profiled approximately every 5m for the top 100m of depth, and every 10-20m between 100–400 m). MHW time series were then created for each year by linearly interpolating the data every 20m. Anomalies for each year were computed by subtracting their POC time series from non-MHW datasets (bsxfun function, Matlab®) and daily-interpolated for display. The first non-MHW dataset was used to estimate 2014 and 2015 anomalies, while the second was used to estimate the 2019 and 2020 anomalies.

To estimate the mixed layer rate of change, we first calculated mixed layer depths (MLDs)⁴³, and then ∆MLD (m month^− 1) for each year as:

∆MLD (m month ^− 1 ) = MLD_month – MLD_month-1 (3)

Where MLD_month is the mean monthly MLD, and MLD_month-1 is the mean monthly MLD of the previous month.

Net community production (carbon export) computed from BGC-Argo floats:

The carbon export at the base of the euphotic zone was computed using float-measured nitrate minus a set of abiotic processes acting on the upper-layer nitrate budget, following the method outlined in Huang et al., (2022)¹⁹:

$$\:{\frac{{\partial\:T}_{\left({NO}_{3}^{-}\right)}}{\partial\:t}|}_{Bio}=$$

$$\:\frac{{dT}_{\left({NO}_{3}^{-}\right)}}{dt}{|}_{Obs}-{\frac{{\partial\:T}_{\left({NO}_{3}^{-}\right)}}{\partial\:t}|}_{DIF}-\:{\frac{{\partial\:T}_{\left({NO}_{3}^{-}\right)}}{\partial\:t}|}_{Adv}-\:{\frac{{\partial\:T}_{\left({NO}_{3}^{-}\right)}}{\partial\:t}|}_{EP}-{\frac{{\partial\:T}_{\left(DIC\right)}}{\partial\:t}|}_{BG}$$

where subscripts of Bio, Obs, DIF,Adv, EP, and BG represent the biological activity, float measurements, diapycnal diffusion, wind-induced Ekman pumping, evaporation and precipitation, and background nitrate gradient change alongside the float trajectory, respectively. The biological term solved from the tracer budget represents the amount of photosynthetically-produced carbon exceeding the respiration consumption, and can serve as a metric of carbon export under the assumption of steady-state when integrating over seasonal or annual scales⁴⁸. The parameterizations of abiotic terms and associated uncertainty estimates are detailed in Huang et al., (2022)¹⁹. The nitrate-based biological term estimate was further converted to carbon unit via the Redfield ratio of 6.6. Due to limited knowledge of the seasonality of diapycnal diffusion coefficients, which may exponentially increase by several orders of magnitude in the fall and winter⁴⁹, and considering that most biological production occurs in the spring and summer¹⁹, our analysis is limited to the relatively stratified spring and summer seasons to minimize error in the biological term estimate.

Phytoplankton composition from pigment concentrations

Phytoplankton pigment concentrations (chlorophylls and carotenoids) were measured by high-performance liquid chromatography (HPLC) following the method detailed in Nemcek and Peña (2014)⁵⁰ and analyzed to estimate the phytoplankton composition using CHEMTAX⁵¹, following the procedures detailed in Peña et al. (2019)¹⁶. Pigment samples have been routinely collected since 2011 on the Fisheries and Oceans Canada (DFO) Line P Hydrographic Surveys (waterproperties.ca/linep) which operate three times a year (generally February, June and August) since 1981 from the coast of British Columbia to station P26 (50°N, 145W) 913 km offshore in the Northeast Subarctic Pacific. For this study, we focus on stations P20 to P26 (Ocean Station Papa), as they span the locations of the BGC-Argo floats (Fig. 2a). Pigment concentration data are from 5 m or mixed layer average when there was more than one depth sampled within the mixed layer. Post-CHEMTAX statistical analyses and data visualization were performed using R Statistical Software (v4.3.2)⁵². The observed patterns in phytoplankton pigment groups were generally consistent across these outermost stations of the Line P transect, with no significant difference by station (PerMANOVA p-value = 0.981, although phytoplankton community composition at P26 was somewhat distinct (Fig. S6). Therefore, average concentrations of major phytoplankton functional groups from hydrographic stations P20-P26 were used to examine trends in phytoplankton over the time series.

DNA-derived Community Composition

2-L seawater samples were collected and filtered onto 0.22-µm Sterivex filters for biomolecular analysis from Stations P20 and P26 during the Fisheries and Oceans Canada Line P seasonal occupation program from 2010 through 2019. While samples were collected throughout the entire water column at each station¹⁵, only samples from 10, 25, 50, 100, 150, 200, and 300 m were used in this study. Briefly, DNA was obtained using enzymatic lysis, followed by phenol:chloroform extraction, and purification and concentration using Amicon Ultra-4 Centrifugal Filter Units (described in greater detail Traving et al., 2021). The V4-V5 region of the 16S rRNA gene was amplified from samples collected for the entirety of the biomolecular time series (2010–2019) using the primers 515F-Y(5′-GTGYCAGCMGCCGCGGTAA-3′) and 926R (5′-CCGYCAATTYMTTTRAGTTT-3′), while the V4 region of the 18S rRNA gene was amplified from 2015–2019 samples using the primers V4F (5’-CCAGCASCYGCGGTAATTCC-3’)⁵³ and V4RB (5’-ACTTTCGTTCTTGATYRR-3’)⁵⁴. For 2010-Winter 2016 samples, 16S V4-V5 amplicon libraries were prepared by the Joint Genome Institute (California, USA) following the JGI iTag Library Preparation SOP (jgi.doe.gov/user-programs/pmo-overview/protocols-sample-preparation-information/), and sequenced on an Illumina MiSeq using a 600 cycle V3 kit. 16S V4-V5 libraries prepared for Spring 2016–2019 samples were prepared at the Hakai Institute (BC, Canada). Samples were amplified in duplicate using fusion primers (containing Nextera indices and Illumina adapters in addition to the primers; described in detail here dx.doi.org/10.17504/protocols.io.3byl4bq1jvo5/v1), PCR replicates pooled, cleaned individually using SPRI beads at ratio of 0.8x, quantified using the Quant-iT PicoGreen dsDNA Assay (Invitrogen), and then all samples successfully amplified were pooled in equal amounts before sequencing on an Illumina MiSeq using a 600 cycle V3 kit (dx.doi.org/10.17504/protocols.io.3byl49842go5/v1). The 16S-V4V5 library generated at the Hakai Institute was sequenced twice and data pooled to obtain a read depth similar to those libraries sequenced at the Joint Genome Institute. The 18S-V4 amplicon library was also prepared and sequenced at the Hakai Institute, using methods described in detail here dx.doi.org/10.17504/protocols.io.rm7vzjjjxlx1/v1 and sequenced on an Illumina MiSeq using a 600 cycle V3 kit.

Illumina sequence data were processed in QIIME2 v. 2022.8⁵⁵. Briefly, primer sequences were removed from amplicon libraries using cutadapt (cutadapt trim-paired) and reads were denoised, amplicon sequence variants (ASVs) determined and chimeras removed using dada2 (dada2 denoise-paired, using the following settings for the 16S-V4V5 rRNA libraries run the Hakai Institute trunc-len-f 240, trunc-len-r 190, max-ee-f 4, max-ee-r 6; trunc-len-f 267, trunc-len-r 220, max-ee-f 3, max-ee-r 5 for the 16S-V4V5 libraries run at the Joint Genome Institute; trunc-len-f 220, trunc-len-r 200, max-ee-f 4, max-ee-r 6 for the 18S-V4 library). As multiple 16S-V4V5 libraries were sequenced, each MiSeq run was denoised separately, ASVs within a fixed length window were retained (rescript filter-seqs-length; global-min 300 and global-max 417), rare ASVs removed (feature-table filter-features; present a frequency of at least 0.001% of the average read depth for that run and in at least 2 samples), and then runs were merged at the ASV level (feature-table merge). Rare ASVs were similarly removed from the 18S-V4 library. ASVs were then classified using the QIIME2 Naive Bayes classifier trained to the SILVA database v.138.1⁵⁶ for the 16S rRNA genes. 18S rRNA ASVs were annotated with both the Protist Ribosomal References database (PR2 v5) and the MetaZooGene database (v2023-m07-15)^57,58. MetaZooGene annotations were used for metazoan sequences while PR2 annotations were used for all other eukaryotes. A total of 3058 ASVs were retained and an average read depth of 119,853 sequences per sample were obtained across the five 16S-V4V5 MiSeq libraries. There was a large difference in the library richness between the sample libraries run at the Joint Genome Institute compared to that run at the Hakai Institute. Rarefying the dataset did not improve this richness difference. Therefore, for the 16S-V4V5 dataset we chose to filter the dataset to retain the top 80% of the ASVs using the function filterfun_sample(topf(0.8)) in the R package phyloseq, leaving 858 16S-V4V5 ASVs. For the 18S-V4 dataset, a total of 3068 ASVs were retained with an average read depth of 56,518 sequences per sample. Raw sequence data have been deposited in the NCBI SRA under the BioProjects PRJNA639229 and PRJNA640752.

All data visualization and statistical analysis of the resulting 16S and 18S rRNA ASV tables were performed using R Statistical Software (v4.3.2). Briefly, depth-based and interannual variation in prokaryotic, Opisthokonta (animal), and non-Opisthokonta (predominantly protistan) community composition was quantified using euclidean distance of CLR-transformed ASV counts agglomerated at the genus level. Principal components analysis (microViz R package)⁵⁹ was used to display the variation in this data. PERMANOVA (adonis2 in the vegan R package)⁶⁰ and pairwise PERMANOVAs (pairwiseAdonis R package, using the Benjamini-Hochberg correction)⁶¹ were used to determine the significance of any observed compositional differences between heatwave years within in the surface and upper mesopelagic waters. A Wilcoxon Rank Sum Test (correcting for multiple comparisons using the Benjamini-Hochberg test) as implemented in the metarr R package⁶² was used to determine significant differences in the Spring and Summer prokaryotic and eukaryotic genera (and all higher taxonomic levels) between MHW years in both the surface mixed layer and upper mesopelagic waters. A genus was required to be present in at least 2% of the samples and have a mean relative abundance of 0.0005 to be included in the analysis. Results were visualized using a heat tree, where the log2 of ratio of the median of proportions for each significantly different sample group is displayed by the color of the nodes and edges while the node size is scaled to the average relative abundance of that taxon in the associated depth range.

Competing Interests

SJH is a co-founder of Koonkie Inc., a bioinformatics consulting company that designs and provides scalable algorithmic and data analytics solutions in the cloud.

Acknowledgments

MBB and KSJ thank support from the David and Lucile Packard Foundation and the GO-BGC project under the NSF Award 1946578 with operational support from NSF Award 2110258. YH was supported by China NSF (grant number: 42406099) and Fundamental Research Funds for the Central Universities (grant number: 20720240105). Data was collected and made freely available by the Global Ocean Biogeochemistry Array (GO-BGC) Project funded by the National Science Foundation, Division of Ocean Sciences, and the International Argo Program and the national programs that contribute to it. The Argo Program is part of the Global Ocean Observing System. A heartfelt thanks to the scientists and crews of CCGS John P. Tully, especially Marie Robert for their world-class assistance during Line P cruises, and Fisheries and Oceans Canada for the financial support of Line P program. Additional support for sample collection and processing was provided by the Tula Foundation in support of CTEK and the G. Unger Vetlesen and Ambrose Monell Foundations, the Natural Sciences and Engineering Research Council of Canada, and the Canada Foundation for Innovation through grants awarded to SJH. We’d like to thank Jacki Long for their assistance with chlorophyll fluorescence corrections. CTEK would also like to thank Rose-Lynne Savage for preparing the 2016–2019 metabarcode libraries for sequencing, and Rute Carvalho for running these libraries. SJH would also like to thank the many undergraduate helpers in the Hallam lab and ocean-going technical support staff including Jade Shiller and Chris Payne for their support in sample collection and processing over the years.

Holbrook NJ et al (2020) Keeping pace with marine heatwaves. Nat Reviews Earth Environ 1(9):482–493
Oliver EC et al (2018) Longer and more frequent marine heatwaves over the past century. Nature communications, 9(1), pp.1–12
Smale DA et al (2019) Marine heatwaves threaten global biodiversity and the provision of ecosystem services. Nat Clim Change 9(4):306–312
Frölicher TL, Fischer EM, Gruber N (2018) Marine heatwaves under global warming. Nature 560(7718):360–364
Bif MB, Siqueira L, Hansell DA (2019) Warm events induce loss of resilience in organic carbon production in the northeast Pacific Ocean. Glob Biogeochem Cycles 33(9):1174–1186
Buesseler KO, Boyd PW, Black EE, Siegel DA (2020) Metrics that matter for assessing the ocean biological carbon pump. Proceedings of the National Academy of Sciences, 117(18), pp.9679–9687
Durkin CA et al (2021) visual tour of carbon export by sinking particles. Glob Biogeochem Cycles, 35(10), p.e2021GB006985
McDonnell AMP, Boyd PW, Buesseler KO (2015) Effects of sinking velocities and microbial respiration rates on the attenuation of particulate carbon fluxes through the mesopelagic zone. Glob Biogeochem Cycles 29(2):175–193
Henson SA et al (2022) Uncertain response of ocean biological carbon export in a changing world. Nat Geosci 15(4):248–254
Bond NA, Cronin MF, Freeland H, Mantua N (2015) Causes and impacts of the 2014 warm anomaly in the NE Pacific. Geophys Res Lett 42(9):3414–3420
Scannell HA, Johnson GC, Thompson L, Lyman JM, Riser SC (2020) Subsurface evolution and persistence of marine heatwaves in the Northeast Pacific. Geophys Res Lett, 47(23), p.e2020GL090548
Amaya DJ, Miller AJ, Xie SP, Kosaka Y (2020) Physical drivers of the summer 2019 North Pacific marine heatwave. Nat Commun 11:1903
Taves RC et al (2022) Relationship between surface dissolved iron inventories and net community production during a marine heatwave in the subarctic northeast Pacific. Environ Science: Processes Impacts 24(9):1460–1473
Long JS, Fassbender AJ, Estapa ML (2021) Depth-resolved net primary production in the Northeast Pacific Ocean: A comparison of satellite and profiling float estimates in the context of two marine heatwaves. Geophys Res Lett, 48(19), p.e2021GL093462
Traving SJ et al (2021) Prokaryotic responses to a warm temperature anomaly in northeast subarctic Pacific waters. Commun Biology 4(1):1217
Peña MA, Nemcek N, Robert M (2019) Phytoplankton responses to the 2014–2016 warming anomaly in the northeast subarctic Pacific Ocean. Limnol Oceanogr 64(2):515–525
Bhatt US et al (2021) Emerging anthropogenic influences on the Southcentral Alaska temperature and precipitation extremes and related fires in 2019. Land 10(1):82
Theys N et al (2019) Global monitoring of volcanic SO2 degassing with unprecedented resolution from TROPOMI onboard Sentinel-5 Precursor. Sci Rep 9(1):2643
Huang Y, Fassbender AJ, Long JS, Johannessen S, Bif MB (2022) Partitioning the export of distinct biogenic carbon pools in the Northeast Pacific Ocean using a biogeochemical profiling float. Glob Biogeochem Cycles, 36(2), p.e2021GB007178
Szuwaski CS, Aydin K, Fedewa EJ, Garber-Yonts B, Litzow MA (2023) The collapse of eastern Bering Sea snow crab. Science 382(6668):306–310
Cartwright R et al (2019) Fluctuating reproductive rates in Hawaii's humpback whales, Megaptera novaeangliae, reflect recent climate anomalies in the North Pacific. Royal Soc Open Sci 6(3):181463
Siegel DA et al (2021) An operational overview of the EXport Processes in the Ocean from RemoTe Sensing (EXPORTS) Northeast Pacific field deployment. Elem Sci Anth, 9(1), p.00107
McNair HM, Morison F, Graff JR, Rynearson TA, Menden-Deuer S (2021) Microzooplankton grazing constrains pathways of carbon export in the subarctic North Pacific. Limnol Oceanogr 66(7):2697–2711
Shea CH et al (2023) Small particles and heterotrophic protists support the mesopelagic zooplankton food web in the subarctic northeast Pacific Ocean. Limnol Oceanogr 68(8):1949–1963
Kwong LE, Ross T, Lüskow F, Florko KR, Pakhomov EA (2022) Spatial, seasonal, and climatic variability in mesozooplankton size spectra along a coastal-to-open ocean transect in the subarctic Northeast Pacific. Progress in Oceanography, 201, p.102728
Cohen NR et al (2021) Dinoflagellates alter their carbon and nutrient metabolic strategies across environmental gradients in the central Pacific Ocean. Nat Microbiol 6(2):173–186
Rivkin RB, Putland JN, Anderson MR, Deibel D (1999) Microzooplankton bacterivory and herbivory in the NE subarctic Pacific. Deep Sea Res Part II 46(11–12):2579–2618
McNair HM et al (2023) Quantitative analysis of food web dynamics in a low export ecosystem. bioRxiv, pp.2023–2003
Decelle J et al (2013) Diversity, ecology and biogeochemistry of cyst-forming Acantharia (Radiolaria) in the oceans. PLoS ONE 8(1):e53598
Lampitt RS et al (2023) Deep ocean particle flux in the Northeast Atlantic over the past 30 years: carbon sequestration is controlled by ecosystem structure in the upper ocean. Frontiers in Earth Science, 11, p.1176196
Coats DW (1999) Parasitic life styles of marine dinoflagellates 1. J Eukaryot Microbiol 46(4):402–404
Giering SL et al (2014) Reconciliation of the carbon budget in the ocean’s twilight zone. Nature 507(7493):480–483
Anderson SR, Blanco-Bercial L, Carlson CA, Harvey EL (2024) Role of Syndiniales parasites in depth-specific networks and carbon flux in the oligotrophic ocean. ISME Commun 4(1):014. .ycae
Steinberg DK, Landry MR (2017) Zooplankton and the ocean carbon cycle. Annual Rev Mar Sci 9(1):413–444
Li L, Endo K (2020) Phylogenetic positions of pico-sized radiolarians from middle layer waters of the tropical Pacific. Progress Earth Planet Sci 7:1–10
Ishitani Y, Ujiié Y, de Vargas C, Not F, Takahashi K (2012) Two distinct lineages in the radiolarian Order Spumellaria having different ecological preferences. Deep Sea Res Part II 61:172–178
Kwong LE, Ross T, Lüskow F, Florko KR, Pakhomov EA (2022) Spatial, seasonal, and climatic variability in mesozooplankton size spectra along a coastal-to-open ocean transect in the subarctic Northeast Pacific. Progress in Oceanography, 201, p.102728
Pond DW, Ward P (2011) Importance of diatoms for Oithona in Antarctic waters. J Plankton Res 33(1):105–118
Wright JJ, Konwar KM, Hallam SJ (2012) Microbial ecology of expanding oxygen minimum zones. Nat Rev Microbiol, 10(6), pp.381 – 94.
Stephens BM, Fox J, Liu S, Halsey KH, Nicholson DP, Traylor S, Carlson CA (2023) Influence of amino acids on bacterioplankton production, biomass and community composition at Ocean Station Papa in the subarctic Pacific. Elem Sci Anth, 11(1), p.00095
West NJ, Landa M, Obernosterer I (2024) Differential association of key bacterial groups with diatoms and Phaeocystis spp. during spring blooms in the Southern Ocean. MicrobiologyOpen 13(4):e1428
Cheng L, Trenberth KE, Fasullo J, Boyer T, Abraham J, Zhu J (2017) Improved estimates of ocean heat content from 1960 to 2015. Sci Adv 3(3):e1601545
de Boyer Montégut C, Madec G, Fischer AS, Lazar L, Iudicone D (2004) Mixed layer depth over the global ocean: an examination of profile data and a profile-based climatology. J Geophys Research: Oceans 109:C12003
Johnson KS, Plant JN, Coletti LJ, Jannasch HW, Sakamoto CM et al (2017) Biogeochemical sensor performance in the SOCCOM profiling float array. J Geophys Research: Oceans 122(8):6416–6436
Roesler C, Uitz J, Claustre H, Boss E, Xing X et al (2017) Recommendations for obtaining unbiased chlorophyll estimates from in situ chlorophyll fluorometers: a global analysis of WET labs ECO sensors. Limnol Oceanography: Methods 15(6):572–585
Briggs N, Dall’Olmo G, Claustre H (2020) Major role of particle fragmentation in regulating biological sequestration of CO2 by the oceans. Science 367(6479):791–793
Bif MB, Long JS, Johnson KS (2024) Seasonality modulates particulate organic carbon dynamics in mid-latitudes of South Pacific and South Atlantic Oceans. Journal of Marine Systems, 241, p.103916
Cronin MF, Pelland NA, Emerson SR, Crawford WR (2015) Estimating diffusivity from the mixed layer heat and salt balances in the North Pacific. J Geophys Research: Oceans 120(11):7346–7362
Emerson S (2014) Annual net community production and the biological carbon flux in the ocean. Glob Biogeochem Cycles 28(1):14–28
Nemcek N, Peña MA (2014) Institute of Ocean Sciences protocols for phytoplankton pigment analysis by HPLC. Can. Tech. Rep. Fish. Aquat. Sci. 3117: 80 p. Fisheries & Oceans Canada
Mackey MD, Mackey DJ, Higgins HW, Wright SW (1996) CHEMTAX - a program for estimating class abundances from chemical markers: Application to HPLC measurements of phytoplankton. Mar Ecol Prog Ser 144:265–283
R Core Team (2023) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
Stoeck T et al (2010) Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Mol Ecol 19:21–31
Balzano S, Abs E, Leterme SC (2015) Protist diversity along a salinity gradient in a coastal lagoon. Aquat Microb Ecol 74(3):263–277
Bolyen E et al (2019) Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol 37(8):852–857
Quast C et al (2012) The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41(D1):D590–D596
Guillou L et al (2012) The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucleic Acids Res 41(D1):D597–D604
Bucklin A et al (2021) Toward a global reference database of COI barcodes for marine zooplankton. Marine Biology, 168(6), p.78
Barnett microViz: an R package for microbiome data visualization and statistics. Journal of Open Source Software, 6(63), 3201, Oksanen J et al (2021) 2024. vegan: Community Ecology Package. R package version 2.6-8, https://CRAN.R-project.org/package=vegan. Martinez Arbizu P. 2017. pairwiseAdonis: Pairwise Multilevel Comparison using Adonis. R package version 0.4.1
Foster ZSL, Sharpton TJ, Grünwald NJ (2017) Metacoder: An R package for visualization and manipulation of community taxonomic diversity data. PLoS Comput Biol 13(2):e1005404. https://doi.org/10.1371/journal.pcbi.1005404

Yes there is potential Competing Interest. SJH is a co-founder of Koonkie Inc., a bioinformatics consulting company that designs and provides scalable algorithmic and data analytics solutions in the cloud.

Extendeddata.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Marine Heatwaves Modulate Food Webs and Carbon Transport Processes

Status:

Version 1

Abstract

Figures

Introduction

Physical and biogeochemical impacts

Food web impacts

Discussion

Online Methods

Temperature anomalies

BGC-Argo float data analysis: POC and chlorophyll estimates

Net community production (carbon export) computed from BGC-Argo floats:

Phytoplankton composition from pigment concentrations

DNA-derived Community Composition

Declarations

Competing Interests

Acknowledgments

References

Additional Declarations

Supplementary Files

Status:

Version 1