Marine heatwaves modulate food webs and carbon transport processes

Domain: earth_science
Published: 2025-10-06
Editorial outcome: major_revision

Abstract
Marine heatwave (MHW) impacts on ecosystem functions and services remain poorly constrained due to limited time-resolved datasets integrating physical, chemical, and biological parameters at relevant scales. Here we show that combining over a decade of autonomous Biogeochemical (BGC)-Argo float measurements with water-column plankton community profiles reveal the impacts of MHWs on particulate organic carbon (POC) production, transformation, and transport in the northeastern subarctic Pacific Ocean. POC concentrations are exceptionally high during the 2015 and 2019 MHWs, linked to detritus enrichment and shifts in plankton community structure. Instead of being rapidly exported to depth, particles <100 µm accumulate in mesopelagic waters, where slow remineralization over the year reduces deep particle flux and carbon sequestration potential. This enhancement is absent in the 2014 and 2020 MHWs, underscoring variability in ecosystem responses to extreme events. These findings highlight the need for sustained, multi-platform observations to assess and predict carbon-cycle responses to thermal extremes.

Introduction
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 cycling1,2. Such impacts can overshoot ecological resilience thresholds resulting in state changes that negatively impact ecosystem functions and services3. Ocean observations and models suggest that MHWs have been expanding in size and intensifying over the past few decades4. Even a reduction in anthropogenic carbon emissions is unlikely to decelerate spatial expansion and duration of these events for many years to come4.

From the standpoint of carbon cycling dynamics, elevated ocean temperature reduces CO2 uptake5 and promotes water column stratification; the latter decreases vertical mixing with a resulting impact on surface nutrient renewal and primary production6. 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 sequestration7–9. 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 limiting10.

Results
Physical and biogeochemical impacts
Both MHWs imprinted a warm and stratified water column effect (Fig. 1a). Indeed, stratification and stability at the base of the mixed layer was strongest during MHW years of 2014, 2015, 2019 and the non-MHW year of 2016, as observed in the Brunt-Väisälä frequency (Fig. 1a). Surface temperature anomalies >1 °C during the first MHW persisted between 2014-15 (Fig. 1b), agreeing with previous reports6,11. Comparisons between three temperature products (BGC-Argo floats, the Argo-based climatology, and the gridded Armor3D reanalysis) confirmed strong agreement among datasets (R=0.88–0.96), with minimal biases and low RMSE values (Fig. 1c, Table S1). The agreement across products suggests that the BGC-Argo data reliably capture regional temperature variability, supporting their use in tracking changes associated with marine heatwaves.

Elevated temperatures propagated to subsurface waters over subsequent years with temperature anomalies up to +0.2 °C observed down to 300 m (Fig. 1d–f). 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 a (Chl), and net community production (NCP) during springs and late summers over the time series, signaling enhanced primary productivity during those seasons (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. 2e). 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-2 in 2020 and 949 mg m-2 in 2021.

POC accumulation during each year’s productive season was estimated in relation to the baseline as ΔPOC=POC -POC . For the first period (Fig. 2d), upper-ocean ΔPOC had similar values regardless of the year corresponded to a MHW or non-MHW year: for example, ΔPOC was 254 mg m-2 in 2015 (MHW) versus 252 mg m-2 in 2011 (non-MHW). During the second period (Fig. 2e), upper ocean POC was anomalously high in the spring of 2019 (MHW) and different from any other year on record. ΔPOC increased to 2041 mg m-2 by April, a 3-fold increase from the baseline and 8x higher than in the MHW year of 2015. The late summer blooms promoted a peak in ΔPOC at 760 mg m-2 in August of 2019, and 1016 mg m-2 in October of 2020 (Fig. 2c). The upper ocean ΔPOC accumulations in 2019 and 2020 were not proportionally tracked by Chl, resulting in a significant drop in Chl:C ratios (Fig. 2c and e).

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). Since this particle accumulation at a deeper depth horizon was unrelated to mixed-layer shoaling and given that the shoaling is minor in the region—only a few meters per month (Fig. S5a)—a mixed-layer pump could not have driven the preferential transfer of particles to deeper waters. Rather, the persistent accumulation at ~200 m likely reflects biologically mediated processes such as vertical migration or mesopelagic retention, which would explain the formation of a stable particle layer at that depth. Additionally, 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 curve6 due to active particle export coupled with microbial respiration.

MHW years experienced much higher accumulation during the spring of 2015 and 2019 than all other years (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 (MHW year, Δ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.

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)20. 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 alone. Summer NCP elevated annual production to about 3.8 mol C m-2 yr-1, but did not contribute significantly to the elevated ΔPOC observed in the mesopelagic (Figs. 2g, 3c). Finally, the 2014 and 2020 summer blooms drove most of the NCP in those years at 2.4 mol C m-2 yr-1 and 0.4 mol C m-2 yr-1, respectively (Fig. S4). The 2020 low NCP did not track the elevated ΔPOC that year (Fig. 2e, g).

Food web impacts
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 the Spring of 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−2019 (2015−2019 only for eukaryotes). Community structure correlated with depth (Figs. S7–S9), which explained 70% of the variation in prokaryotic communities (PERMANOVA, p=0.001) and 24% of the variation in eukaryotic communities (PERMANOVA, p=0.001). Season only explained <10% of the community variation (PERMANOVA, p=0.001), perhaps not surprising given the weak seasonal harmonic in biochemical parameters at P20 and P2626.

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 cyanobacterial genus Synechococcus. 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, signaling a complex food-web26–29. 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, however copepod sequences as a whole were notably lower in the 2019 upper mesopelagic layer compared to 2015 (Fig. S10). Typically found on the continental shelf in the Northeast Pacific30, the presence of sequences closely related to Paracalanus sp. warrants further investigation and validation.

The planktonic community inhabiting surface waters was relatively consistent between MHWs (Table S2) 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 S2), 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.

Methods
Temperature anomalies
The depth-resolved physical parameters shown in Fig. 1 were obtained from an Argo-based gridded product46- (RG09) Argo Climatology (https://sio-argo.ucsd.edu/RG_Climatology.html, last access: 06 January 2024). We selected data from Ocean Station Papa (50°N, 145°W), a core region used to define MHWs in NESAP11–13, to assess decadal-scale variability. Temperature anomalies at each depth were calculated as deviations from the climatological mean during the period 2004−2023, 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 meters47, which generally corresponds well with the depth at which the maximum Brunt-Väisälä frequency (N2) occurs (calculated using sw_bfrq Matlab® function).

To verify the consistency of temperature estimates used in the anomaly calculations, and to assess potential spatial biases arising from the semi-lagrangian nature of the BGC-Argo floats, we compared monthly temperature records at 50 m from three sources: BGC-Argo floats, the gridded Roemmich and Gilson (2009) Argo climatology product, and the Armor3D Copernicus reanalysis48,49. The BGC-Argo data were averaged monthly within the study domain. Armor3D temperature fields were averaged over the spatial bounds of 47°N-53°N and 147°W-137°W, corresponding to the spatial extent sampled by the BGC-Argo floats, using monthly reprocessed data.

BGC-Argo float data analysis: POC and Chl 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–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 published50,51. As previously described15, float-estimated chlorophyll a (Chl) in our study region is biased when compared to nearby discrete Chl measurement from repeat Line P cruises, requiring additional post-correction. We applied a correction factor to the adjusted Chl data for each float, estimated by taking the median of the ratio of float-adjusted Chl data to the Line P corrected float data15 in the upper 15 m. The following correction factors were divided by our adjusted Chl data to obtain corrected values extending to 2022.

Backscatter (bbp) data was separated into small and large particles52,53 (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 bbp 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 zone7,53 as:

Ez(m)=Z of 0.1*Chl Max below Z chl Max (1)

Where Chl_Max is the depth of maxima Chl 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 bbp into total POC50.

POC(mg=m3)=31,200×bbp+3.04 (2)

Chl:C ratios were estimated using mean monthly integrated POC and Chl 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 20 m (floats profiled approximately every 5 m for the top 100 m of depth, and every 10−20 m between 100−400 m). MHW time series were then created for each year by linearly interpolating the data every 20 m. 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)47, 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.20:

∂T(NO3)j / ∂t = ∂T(NO3)j / dt Obs - ∂T(NO3)j / ∂t DIF - ∂T(NO3)j / ∂t Adv - ∂T(NO3)j / ∂t EP - ∂T(DIC)j / ∂t BG (4)

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 scales54. The parameterizations of abiotic terms and associated uncertainty estimates are detailed in Huang et al.20. 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 winter55, and considering that most biological production occurs in the spring and summer20, our analysis is limited to the relatively stratified spring and summer seasons to minimize error in the biological term estimate.

Spatiotemporal sensitivity analyses
We conducted a series of sensitivity analyses to assess whether methodological assumptions or float-specific variability introduced bias into the observed patterns of POC production and accumulation. First, to test the robustness of using a fixed 100 m integration depth for surface POC stocks, we compared estimates based on a variable Ez versus the fixed 100 m depth. Small-particle backscatter (bbp) was integrated from the surface to both depths and converted to POC using empirical coefficients (Eq. 2). The resulting regression analysis revealed a strong linear relationship (R²=0.94) with a slope of 0.96 (Fig. S11), supporting the use of a fixed depth for consistent comparisons across the time series.

Next, we tested whether float-specific spatial variability could influence our results. We compared Chl concentrations from two floats deployed during overlapping periods but spatially apart (#5903274 and #5903714; Figs. 2a, S12a). Profiles collected within 3 days of each others showed a strong positive correlation in integrated Chl over the upper 100 m. We also calculated NCP from all floats equipped with nitrate sensors using the same method as in Fig. S4, and included floats lacking bio-optical sensors but profiling the region concurrently (Fig. S12b). A linear mixed-effects model with year and float ID as random intercepts showed that most variance in NCP was explained by year (σ²=0.81), while float ID contributed negligible variance (σ²≈0), indicating that floats operating during the same time periods yielded comparable NCP estimates.

Finally, we examined spatial variability in surface chlorophyll using monthly mean concentrations from the Aqua-MODIS satellite product56 (4km resolution, 2008-2023) across four Line P stations located within the float trajectories. Seasonal chlorophyll patterns were highly consistent among stations (Pearson correlation coefficients=0.85–0.98, Bonferroni-adjusted p-values<0.05; Fig. S12c), and within-month differences across stations were not statistically significant (Repeated Measures ANOVA, F(4,50)=0.213, p=0.93). Discrete chlorophyll samples collected during Line P cruises further confirmed these seasonal trends (Fig. S12d). While individual profiles were strongly correlated regionally (P21 and P22, Spearman’s ρ=0.77, adjusted p-value<0.001; P24 and P25 Spearman’s ρ=0.89, adjusted p-value<0.001), there was weaker correspondence between more distant stations, reflecting the inherently patchy nature of chlorophyll distributions, these results collectively indicate that the observed trends in POC production and accumulation were not driven by float-specific or spatial biases.

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ña57 and analyzed to estimate the phytoplankton composition using CHEMTAX58, following the procedures detailed in Peña et al.17. 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 averages 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)59. 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 station16, 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 (describing in greater detail Traving et al.). 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’)60 and V4RB (5’-ACTTTCGTTCTTGATYRR-3’)61. 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 https://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 (https://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 https://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.862. 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.163 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)64,65. 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 data set 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)66 was used to display the variation in this data. PERMANOVA (adonis2 in the vegan R package)67 and pairwise PERMANOVAs (pairwiseAdonis R package, using the Benjamini-Hochberg correction)68 were used to determine the significance of any observed compositional differences between heatwave years within 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 package69 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript examines the effects of marine heatwaves (MHWs) on carbon cycling dynamics in the Northeastern subarctic Pacific Ocean. The authors used over a decade of Biogeochemical-Argo float observations with data about the plankton community composition, effectively integrating complementary datasets to address this pressing ecological issue. The study reveals that MHWs in 2015 and 2019 resulted in elevated particulate organic carbon (POC) concentrations, driven by increased detritus production and shifts in plankton communities. Interestingly, MHWs in 2014 and 2020 did not exhibit similar impacts, highlighting the variable nature of these events.\nThe manuscript is well-structured and addresses a significant gap in understanding the effects of MHWs on marine ecosystems, making it a valuable contribution to the field. However, given the broad readership of Nature Communications, certain sections could benefit from revisions to improve clarity and accessibility. Additionally, some figures would benefit from refinement (see minor comments). Beyond these points, I have several major concerns regarding specific methods and data processing approaches that the authors should address.\nMajor comments:\nA) My first major concern with the authors' approach is their use of data collected from different platforms and locations within a relatively large area (47˚–52.5˚N and 137.5˚–146˚W). For instance, a recent study by Niebergall et al. (2023), which analyzed data collected in the same region (near Ocean Station Papa), concluded that “the primary driver of differences between the datasets (i.e., they evaluated new production and NCP from thirteen independent datasets) is spatial variability around Station P.” Furthermore, in Figure S6 of this manuscript, it is evident that the diatom spring bloom in 2019 observed at Station P20 is not detected at Station P26 (as noted by the authors, line 431), further highlighting the significant spatio-temporal variability within this region.\nSuch variability could potentially explain the differences in the observed impacts of marine heatwaves in different years (e.g., 2014 and 2020 versus 2015 and 2019).\nI strongly recommend a more thorough analysis and consideration of the spatio-temporal variability present in the area. For example, comparing similar data measured by BGC-Argo floats, cruise data, and remote sensing observations (e.g., chlorophyll-a concentration, temperature, and salinity) could provide valuable insights and strengthen the study's conclusions.\nB) My second major concern relates to the way in which physical data have been used to define MHWs:\n1b)Author decided to use an “Argo-based gridded product” associated with the paper Cheng et al. (2017). However, no link (or doi) is provided to know which version of the data was used. For example, note that the product is available in two different horizontal resolutions (0.5° and 1°) and that the 1° resolution data were updated in January 2024.\n2b)Because this product only provided temperature data, authors could not calculate density, which is typically used to evaluate the Mixed Layer Depth (MLD) and the Brunt-Väisäla frequency (N2). However, there are other gridded products that provide temperature and salinity:\n•EN4 (https://www.metoffice.gov.uk/hadobs/en4/index.html)\n•ARMOR3D (https://data.marine.copernicus.eu/product/MULTIOBS_GLO_PHY_TSUV_3D_MYNRT_015_012/description)\n3b)Authors decided to extract from one of the gridded products from Cheng et al. (2017), the temperature data from one location: at the Ocean Station Papa (50°N, 145°W). But their manuscript use data from BGC-Argo floats that cover a larger region between 47˚-52.5˚N and 137.5-146˚W, which probably present a spatio-temporal variability that should be considered. Moreover, no comparison with the temperature measured by the float or during the Line P cruises are performed to ensure that the floats and the cruises effectively collected data during marine heatwaves and in similar water masses.\n4b)The quality of Figure 1 is relatively poor. The x-axis in panel (a) is not aligned with the other panels of the figure. Interpolation should have been carried out in panels a-b, instead of representing the data in a scatter plot where the circles are enlarged and overlap. Panels c-e are not really useful and are not mentioned in the manuscript.\nC) My third major concern is the processing of data from the BGC-Argo floats:\n1c)Line 361: “only small particles were used in this study since large particles were a minor fraction of total bbp and showed no temporal variability (Figs. S1-S3)”. I disagree with the authors' interpretation because the range of values (colorbar) is consistent between the graphs for small and large particles (middle panels in Figures S2 and S3). Additionally, there appears to be interannual variability in the concentration of large particles. Could this variability be significant and/or linked to the marine heatwaves?\n2c)Approximating the productive surface layer as being between 0 and 100 m is crude. This fixed definition of the surface productive layer also impacts the definition of the mesopelagic layer (line 199: “…100-300 m depth, denoting the upper mesopelagic”).\n- Has the definition of the euphotic zone (line 366) been published before? Some BGC-Argo floats are measuring Photosynthetic Available Radiation, which could be used to calculate the depth of the euphotic zone, and to compare it with this definition. It could also be compared with estimates based on remote sensing data.\n- Line 367: “Mean euphotic zone depth was 98 ± 20 m for the entire time series, so we assumed Ez = 100 m for the subsequent calculations”. Although that there is an assumed variation of the depth of the euphotic zone between 80-120m, the authors decided to define the surface productive zone between 0-100m. However, with BGC-Argo float data, it is possible, for a given profile, to define the depth of the productive layer more precisely as the greater of the MLD or the depth of the euphotic zone.\n3c)In a recent published study, Niebergall et al (2023), that used data collected in the same area (near Ocean Station Papa), the NCP was calculated with NO3, O2 and DIC data estimated by BGC-Argo floats. But here the authors only used nitrate data. Why? Using another estimate could provide an uncertainty and give more confidence in their conclusions.\nMinor comments:\n1) Abstract: After reading the manuscript, I find it surprising that the years 2014 and 2020 are not mentioned in the abstract. It would be valuable to at least acknowledge that the cascading effects of marine heatwaves in this oceanic region may be more complex and variable than currently described in the abstract.\n2) Line 51: The study by Mignot et al (2022", "reviewer_2": "This review is based on the evaluations of two reviewers. As part of the Nature Communications initiative to support Early Career Researchers, an Early Career researcher also participated in this co-review process.\nStatement from the Early Carrier researcher: “I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.”\nGeneral comments\nThe manuscript “Marine Heatwaves Modulate Food Webs and Carbon Transport Processes” by Biff et al. investigates the anomalous POC accumulation that occurred in the mesopelagic zone of the Northeast subarctic Pacific during 2 Marine Heatwave (MHWs) events in 2015 and 2019. The authors use a very interesting dataset comprised from both BGC-Argo deployed in the area and in-situ data of plankton community composition (Pigments and metabarcoding data). I find this kind of multidisciplinary approach very interesting and I believe is one of the most appropriate ways to properly address ecological questions in areas like NESAP and generally the open ocean, where access, and thus in-situ data, is limited.\nHowever, I believe that there are some very important issues that the authors should try to resolve. Particularly the connection between MHW-induced community composition changes and the impacts of zooplankton vertical migration on the abnormal particle accumulation during 2015 and 2019 is not clear and cannot be unambiguously deduced from the results presented. Moreover, alternative mechanisms that could contribute to the abnormal POC accumulation in the mesopelagic zone, such as the Mixed Layer Pump are not considered at all. Therefore, I believe that it is essential to address the outstanding issues to provide a robust and convincing explanation for the anomalous POC accumulation under MHW conditions before the manuscript can be considered for publication.\nSpecific comments\nLines 79 – 81: “During both MHWs…..in the region”.\nI cannot see how the accumulation of small particles in the mesopelagic zone was actually related to an intensification of zooplankton vertical migration. Although according to the references provided by the authors (23, 24) zooplankton vertical migration contributes to POC in the mesopelagic zone, the authors provide no clear evidence showing that in the cases of MHWs examined the accumulated POC was actually due to zooplankton vertical migration. Another potential source of small particle supply to the mesopelagic zone which is the mixed layer pump (Omand et al., 2015; Dall’Olmo et al., 2016; Lacour et al., 2019) that can have a significant contribution to the carbon export in the Northeast Pacific (Nowicki et al., 2022) is not considered or at least discussed.\nLine 86: “…stratification at the base of the mixed layer was strongest during MHW years \"(2014-16 and 2019, Fig.1a).”\nWouldn’t that help isolate the euphotic zone from the mesopelagic zone thus promoting the accumulation of Organic Matter within the latter?\nLines 120-121: “POC accumulation during each year’s productive season was estimated in relation to the baseline as ∆POC = POCmax - POCmin .”\nA table is needed showing ΔPOC for each year so the reader can directly compare the values for each year.\nLines 121-124: “For the first …. in 2011 (non-MHW).”\nDo you mean for the first period? Better use consistent phrasing. Also, the reference to Figure 2 panels is confusing (I suppose you talk about Figure 2d).\nLine 124: “During the second half of the time series…” same issue as in the previous comment. Please use consistent phrasing.\nLine 128: “…a significant drop in Chl:C ratios (Fig.2 c,e).”\nIn spring 2019 Chl:C was lower than all other years during the second period (and maybe the first but is hard to tell from the figures due to different y-axis limits). I am wondering if this is related to the anomalously high diatom concentration observed during spring 2019 (Fig. 4). (Such high diatoms concentration was also observed during spring 2018 but no BGC-Argo data are available then.).\nLines 128-131: “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.”\nThis part is not clear - Please rephrase. It took some time to understand where the sentence was referring to. Maybe a better correspondence between the text and the figure it refers to would help. Also, it has been stated earlier in the text that ΔPOC is max-min for each year, so how is it different on April 2019 than in late summer? Maybe you are referring to POC and not the ΔPOC?\nLine 132 - 135: “The productive seasons were followed ….. during spring months (Fig. S5).”\nThis is hard to see from Figure S5. A alternative better plot for Fig.S5 might be one single panel with bars of ΔMLD and on top to plot the interannual variability of POC 200m/POC 100m ratio.\nLines 145-146: “The high ∆POC feature was … the spring peak.”\nAs I stated before, I don’t understand how you compare seasonal ΔPOC values, since ΔPOC is calculated on an annual basis. Please rephrase or provide further details.\nLines 169-171: “There was a mismatch … to increasing phytoplankton productivity and biomass.”\nThere are indications though of increased contributions of diatoms during spring 2019 (Fig. 4) that have fast sinking rates and would result in increased export fluxes to the mesopelagic. Also, according to Fig. S1 Chl-a was particularly high during spring 2019 which is indicative of enhanced phytoplankton growth despite the moderate NCP estimates.\nLines 202-204: “As expected, community structure …..of the community variation (p = 0.001)…”\nAgain I cannot understand how this conclusion was reached. For example, it is true that the first two principal components in Fig. S7 explain 70% of the variability between samples observed but it is the community structure that was used to build the principal components, therefore this sentence is not correct or at least these results cannot be deduced from Figs S7-S9. Please provide further information. Also, in Fig. S7 the panel labels (A-F) are messed up. Please correct them.\nLine 277-278: “As thermal stress persisted …dark ocean waters (Figs. S4, 4, 5).”\nFig. S4 shows NCP and Fig. 4 shows phytoplankton community structure. How are they relevant to this sentence?\nLine 278-281: “For instance, the ….water column29,30. ”\nIt is true that Cladococcus was particularly elevated during 2015 MHW but according to Fig. S9 it was also elevated during non-MHW years. Also, both Cladococcus and Acantharea were elevated during August, that is at the ending phase of the anomalous POC accumulation, thus it cannot unambiguously be shown that the proliferation of these groups was the cause of the anomalous POC accumulation in the mesopelagic zone during 2015 MHW.\nLine 281-284: “Moreover, increased … RAD-B radiolarians33,34 .”\nI cannot see the increased abundance of Synidiales in Fig. S9. Their contribution appears similar under MHW and non-MHW years. Could you please be more specific?\nLines 284-287: “Proliferation of ……….deep particle export.”\nAgain, here Oithona copepods were also particularly elevated during non-MHW years according to Fig. S9, therefore this result does not provide support the linkage between their presence and the abnormal POC accumulation under MHW contitions. Also, the reference to Fig. 1 is irrelevant here.\nLine 294-296: “Interestingly, elevated NCP (Fig. S4) and high-chlorophyll levels in 2020 (Fig. 4a) promoted by wildfires and volcanic eruptions14,​17,​18 did not result in mesopelagic POC accumulation (Fig. 3c), …”\nThis is not correct. According to Fig. S4 NCP was not at all elevated during 2020. Additionally, according to Fig. 4a, Chl-a does not appear particularly elevated during 2020. It is lower than both MHW years (2025 and 2020), although spring sampling of 2020 is unfortunately missing.\nLines 288–299: While the increased thermal stratification during MHW clearly contributed to the significant retention of POC in the mesopelagic zone, the presented results do not allow for conclusions about changes in food web structure or carbon processing.\nMinor comments\nLine 43: ‘’as a function of climate change.’’: I think as a result is more appropriate wording.\nLine 54: “NESAP”-Please redefine abbreviation in the main text. Currently the definition is only in the Abstract.\nLine 88 -92: You may add a reference to Fig. 1b-e.\nLine 94: A map would be useful in Fig.1 showing the domain of the study area. Also, please add (black line) next to “overlain by mixed layer”.\nLine 99-100: “maxima in chlorophyll (Chl) fluorescence concentrations” I think it is more accurate to use the term Chlorophyll-a (Chl-a). Also please delete “fluorescence”. It is Chl-a concentration.\nLine 101-102: “…signaling enhanced primary productivity (Figs. S1-S4)”. Do you mean during 2015 and 2019, i.e., under MHW conditions? If yes please state the years explicitly.\nLine 137: “…2015 and 2019 (Fig. 2f,g).” Please add: “compared to all other years” for clarity.\nLine 173: Fig. 3d is not relevant here, maybe you mean Fig. 3c.\nLine 179: Please be consistent on how you are referring to Chl-a. Sometimes abbreviated, sometimes not, sometimes italics, sometimes not.\nLine 189: (Fig. 4 legend). I assume you show Chl-a (total and per chemotaxonomic group as estimated by chemtax). Please fix legend and axis labels accordingly.\nLine 210: Add space between “genus” and “Synechococcus”\nLine 249: “Heat trees showing changes…” Please be more specific. I think it should clearly state “…changes between 2015 and 2019 MHWs.”\nLine 315 – 317: “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.” Please add “ relative to 2015 MHW” to be clearer.\nLine 352: “As previously described​15, …” Please check the reference number. I think you probably you mean Ref. 14 (Long et al., 2021).\nLine 504: “…within in…” delete “in”.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We appreciate the reviewer’s comment. We conducted a sensitivity analysis and now show an additional figure (S11) in the supplementary material. We also added a paragraph on the methodology section, talking about this analysis, starting around Line 574. The seasonal patterns and magnitudes of the different variables analyzed are not significantly different among nearby stations. We compared different floats, stations, and even used a satellite product to look at spatial variability. Therefore, we believe that our dataset is sufficient to explain the interannual variability in the constrained region and hope that Fig.S11 suffices the reviewer’s concern.\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "One comment may not have been fully resolved. Specifically:\nComment:\nLine 294–296: “Interestingly, elevated NCP (Fig. S4) and high-chlorophyll levels in 2020 (Fig. 4a) promoted by wildfires and volcanic eruptions14,​17,​18 did not result in mesopelagic POC accumulation (Fig. 3c), …”\nThis is not correct. According to Fig. S4 NCP was not at all elevated during 2020. Additionally, according to Fig. 4a, Chl-a does not appear particularly elevated during 2020. It is lower than both MHW years (2015 and 2023), although spring sampling of 2020 is unfortunately missing.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for the observation. The correct period is the summer of 2019, not 2020. The text is now updated.\nHowever, we could not find this updated information reflected in the main text.\nWe apologize for the mistake. The text is now updated around line 344: “Interestingly, elevated NCP (Fig. S4) and high Chl levels in 2019 (Fig. 4a) promoted by wildfires and volcanic eruptions15,18,19 did not result in mesopelagic POC accumulation (Fig. 3c), suggesting that phytoplankton production alone cannot account for anomalous ∆POC accumulation.”"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed several comments from the reviewers; however, they have not thoroughly responded to some major concerns. As a result, I have significant reservations about whether the current manuscript has reached the level required for acceptance in Nature Communications.\n[...]"}, "author_response": ""}]

Decision letter