Dinoflagellate vertical migration fuels an intense red tide

Significance Extremely dense harmful algal blooms (HABs) are an increasing problem globally. How microscopic, single-celled organisms can reach such high abundances is still poorly understood. Over 50 y ago, it was postulated that motile dinoflagellates could form dense blooms through their vertical swimming, which gave them a competitive advantage over other phytoplankton. We tested this hypothesis with innovative in situ observations during a dense bloom. The motile dinoflagellates swam downward at night into the deep nutrient pool, and upward during the day to photosynthesize. Decreases of nutrients in deep waters were proportional to increases in organism abundance, directly linking organism behavioral and metabolic activities. These observations demonstrate that vertical migration is central to the formation of exceptionally dense dinoflagellate HABs.

• Step 1: The SUNA was configured to have a duty cycle of 20 minutes per hour with a sampling frequency of ca. 1 Hz.Raw nitrate data (unit: μM; 1 μM = 1 mmol m -3 ) was recorded as a function of time.The 95% confidence interval of raw nitrate data (determined by calculating the probability density function of nitrate data within a 0.5 o C range) was 1 mmol m -3 (about 2 standard deviations).This value was used in subsequent error calculations.
• Step 2: The raw nitrate data contained outliers, such as unexpected zeros or extreme values (about 5% of the total amount of the data).These points were removed.
• Step 3: The SUNA V2 pump stream was oriented with the intake facing in an upwards direction.For this reason, only measurements collected during the upward, free-ascent travel of the Wirewalker were used in this analysis. • Step 4: The nitrate data exhibited a temporal drift -increasing by ca. 3 mmol m -3 over 8 days consistently across all depths (Fig. S3a).This consistency indicates this temporal drift is an instrument error and thus needs to be accounted for.The minimum nitrate value within each duty cycle was identified.Then a low-pass filter with a cut-off frequency of 0.5 cpd was applied to the time series formed by these points to obtain a smoothed baseline.
Assuming that the actual minimum nitrate concentration was consistently 0 over the course of the deployment, this baseline was subtracted from the upcast data to bring the minimum values back to 0, as shown in Fig. S3b. • Step 5: The nitrate-temperature relationship using the temporal-drift-corrected nitrate data showed a positive relationship between nitrate and temperature above 12 o C (Fig. S3d, also appearing as an unrealistic non-monotonic vertical nitrate gradient near the ocean surface in Fig. S3b).Previous literature (e.g., 3) has also shown a temperature-dependent response of the nitrate sensor.This temperature effect was corrected by first fitting a linear function to the temperature and nitrate data in the temperature range where nitrate is expected to be 0 (here above 12 o C), and then subtracting this linear function from the data.
The result is shown in Fig. S3c-d, where nitrate is consistently 0 above a particular temperature threshold, as expected from previous data (e.g., 4).
• Step 6: Given the ca. 1 Hz sampling rate of SUNA, and an average vertical speed of 0.3 m s -1 for WW, nitrate data for each WW profile were gridded into uniform 1-m depth bins from the surface down to 100 m depth.

ChlorophyII fluorescence (ChlF) data analysis
Non-photochemical quenching (NPQ) is a photoprotective mechanism that suppresses fluorescence under high irradiance.Hence, in situ fluorescence data need to be corrected as NPQinduced fluctuations do not represent real changes in biomass.Adapted from other NPQ correction methods applied to a diversity of profiling platforms, including the Wirewalker (e.g., 5-7), we estimated NPQ by first calculating the correlation coefficient between ChlF and PAR using all the data between the surface and 4 m in a two-day window.If the correlation coefficient was less than -0.1, the NPQ value for the full water column would be calculated as a linear function of PAR with the transfer function coefficients estimated from the linear fit of ChlF and PAR in the upper 4 m.The same process was applied incrementally to the whole data set using a 0.1-day offset.The final NPQ value was the average across all the increments that contained that time point.The NPQcorrected ChlF data was obtained by adding the NPQ correction to the raw ChlF data.
ChlF was initially recorded in units of relative fluorescence units (RFU).To convert ChlF into the more common chlorophyII-a (Chl-a) concentrations (units: mg m -3 ), a linear fit was calculated between the WW 31-33 m depth-averaged ChlF time series (RFU) and the nearby Del Mar Mooring 32 m calibrated Chl-a time series (mg m -3 ) (correlation coefficient of 0.91, plot not shown here).This regression was used to estimate a conversion function for Chl-a from RFU to mg m -3 : where  represents Chl-a in units of mg m -3 and  represents ChlF in RFU.Eq.1 was applied to all the NPQ-corrected ChlF measurements.

The light field, with optical backscatter
The photosynthetically available radiation (PAR) was estimated by integrating the observed irradiance at wavelengths of 380 nm, 412 nm, 490 nm, and 532 nm and then interpolating/extrapolating to the range between 400 nm to 700 nm (i.e., the full PAR range).The equation is: where  is the measured irradiance in units of   '  ' ,  represents each wavelength,  is the wavelength interval size, and 0.046 is for conversion from   ' to   '  ' .Fig. S4a shows that even though the surface PAR had some spikes, the shortwave radiation (measured independently at the Scripps Pier) was about two to ten times larger than surface PAR -suggesting that about 50%-90% of the input solar energy was absorbed in the surface 3 m.Below the ocean surface, PAR was modeled as decaying exponentially in depth: where   is the PAR value at the ocean surface,  is the diffuse attenuation coefficient of irradiance, and  is depth.For each WW profile,  can be estimated from linear fits of [()] vs. depth, giving values between 0.08 m -1 and 0.6 m -1 .The euphotic depth (  ), defined as the depth where PAR is 1% of its surface value, was estimated from: A larger absolute   means that the water is clearer, with deeper penetration of light.A negative correlation between   and upper ocean (0-30 m) phytoplankton biomass (or particle load) (Fig. S4b) demonstrates that the subsurface light field was strongly modulated by upper ocean biological properties during the red tide: that is, the irradiance was attenuated (and the euphotic zone shoaled) by the dense phytoplankton community in addition to the expected decay of light in clear water.
Turbidity (suspended particle load) was estimated from the optical backscatter data using the instrument's factory calibration: where  is in NTU (Nephelometric Turbidity Units), and  is the measured optical backscatter at 532 nm.Fig. S4d showed similar spatial-temporal patterns as Chl-a (Fig. 1e), indicating that the suspended particles in the upper ocean were dominated by phytoplankton during the observational period.Note that we use turbidity as a relative measure, and the sensor is assumed to be stable in time.

Depth-to-isopycnal coordinate transformation
Raw data were recorded as a function of depth (pressure) and time, as seen in Fig. 1c-f.Within a depth range, raw measurements have variations both along and cross the isopycnals as a function of time.An example of nitrate concentrations in T/S space (Fig. S5) shows both a general trend of decreasing nitrate with decreasing density (i.e., cross-isopycnal variability), and variations of nitrate over time at the same density value (i.e., along-isopycnal variability).When internal waves were present, cross-isopycnal variability appears as temporal changes of nitrate (and other variables like temperature and salinity) in the depth coordinate system due to the vertical heaving of isopycnals.This complicates calculations of temporal changes at any depth, especially for dynamically passive but biologically active tracers such as nitrate.Transforming from depth coordinates to isopycnal coordinates removes the vertical displacements created by internal waves, decoupling the fluctuations induced by vertical displacements of isopycnals from those caused by biological uptake and horizontal (along-isopycnal) gradients.
The depth-to-isopycnal coordinate transformation was executed by interpolating raw physical and biogeochemical variables from each individual WW profile onto uniformly gridded density bins (here, the density interval was set to be 0.02 kg m -3 ).Similarly, raw pressure measurements were also interpolated to give the instantaneous depths of each density bin for each profile.The mean isopycnal depth was calculated as the time-averaged interpolated depth for each gridded density bin over the course of the study period.As the mean isopycnal depth accounts for the different vertical spacings between different density layers, it provides a more realistic representation of the vertical distribution of density layers than using density as the y axis.
This different spacing of density bins in depth also leads to different numbers of raw measurement points in each gridded density bin, which results in different errors associated with the transformed values.To estimate the errors, the standard error of the mean (SEM; the standard deviation of the mean of multiple samples means) was calculated for each density bin: where  is the standard deviation of the raw measurements, and  is the number of points within each density bin.An example of the coordinate transformation result is shown in Fig. S6, where vertical displacements of the isopycnals disappear after applying the depth-to-isopycnal transformation.While temperature and salinity are relatively constant along isopycnals (Fig. S6b,d), the along-isopycnal nitrate displays considerable variability (Fig. S6f).

Along-isopycnal velocities
As Eulerian measurements (fixed point mooring), the along-isopycnal WW time series for each variable contains temporal gradients driven by horizontal advection of existing spatial gradients.As seen in Fig. S7c, from noon 05/06 to 05/08, the along-isopycnal nitrate and horizontal velocities (southeast-ward current corresponding to decreasing of nitrate) were correlated, suggesting that changes in the along-isopycnal nitrate were partially driven by horizontal advection.Crucially, the advective signal will confound estimates of biologically driven nitrate rates of change, and thus must be accounted for.The detailed procedures are shown below.

Temperature/Salinity (T/S) binning
To separate the effects of horizontal advection from biological activities in the along-isopycnal nitrate time series, we applied a Temperature/Salinity binning method, which utilized spiciness as a tracer to follow the same water parcel (defined by the conservative tracers temperature and salinity) in time.Spiciness is a quantity related to Temperature/Salinity co-variability on a density layer: more positive means warmer and saltier water (spicier), while less positive means cooler and fresher water (mintier) (8).Thus, temporal changes in nitrate concentration driven only by biological processes can be estimated by calculating changes in nitrate concentration within each water mass.An example of how this method was performed is shown in Fig. S8.Along each isopycnal, data points were rearranged into evenly spaced spiciness bins (interval of 0.005 kg m -3 ).Points in the same spiciness bin were assumed to represent the same water mass.Note that the number of data points in each spiciness bin will vary.

Estimation of Nitrate rate of change
Rates of change of nitrate,   (units: mmol N m -3 day -1 ), quantify how fast nitrate changes over time, and are estimated from the specific nitrate rate of change,  (unit: day -1 ), following the procedures described below.
After the T/S binning, within each water parcel, nitrate can be modeled as changing exponentially over time (similar to (9)).The equation for the temporal evolution of nitrate can be expressed as: where   is the initial nitrate concentration,  is the nitrate-specific rate of change of nitrate (days - 1 ), and  is time (days).Using consecutive time points in each spiciness bin (i.e., each water mass),  was estimated as: where  is the measured nitrate concentration, and  represents the time index of the data.The error for  was calculated as: An example is shown in Fig. S9a.Note that this exponential fitting was only applied when the absolute difference between two points was larger than 2 mmol m -3 (2 times raw nitrate data's confidence limit) -ensuring a statistically significant estimate of  that can be further interpreted as biological processes.Due to this criterion, some estimates of  were undefined.Overall, the probability density functions (PDFs) of the estimated  (Fig. S9d) show a clear distribution of the majority of the data being negative, consistent with the observations of nitrate loss over time.Furthermore, estimates of  give a distinct mode value of -1 day -1 (Fig. S9d), indicating that the biologically driven loss rate of nitrate is relatively constant among different water parcels, across different isopycnals, and over time.Although experiments were done under different conditions (lab vs. nature), and with different species (Heterocapsa niei vs. Lingulodinium polyedra), our estimated  is consistent with the values reported in (9).
A temporal integration of the estimated (, ) was performed to estimate the changes in the nitrate field driven only by biological activities -that is, without the effects of advection.Due to the fact that estimates of  had uneven time intervals along each isopycnal (primarily due to the inconsistent presence of different water masses at the WW), a rebinning in time was executed by setting a uniform time bin of 4 hours, and choosing the median  value within each 4-hour window to be the rebinned , denoted as   .Its error is: where  is the number of estimates of  in each time bin.Even though the time resolution of the nitrate estimates decreased after this rebinning process, there were two distinct advantages: 1) the process generated a uniform time interval for further nitrate calculations; 2) the process compensated for non-significant values previously excluded in the estimates of .
The predicted biologically mediated nitrate concentration,   was calculated using the estimates of  along each isopycnal as: where  is the time interval between estimates of  (4 hours),  represents the index of the data, and the initial value is taken from the first profile shown in Fig. 2c from the observations.The error associated with the predicted nitrate is estimated as: is shown in Fig. S9b, with signatures of clear deepening of nitrate contours in time, and also decreasing of nitrate along isopycnals in time between 20 m and 33 m.These patterns are consistent with the observations (Fig. 2c), but have had advection-associated fluctuations removed.
To estimate the rate of change of nitrate,   , a linear differentiation in time was applied to the predicted nitrate concentrations (Fig. S9b) along each isopycnal: . [13] The associated error was estimated by: A modal value of -6 mmol m -3 day -1 was obtained from the PDF of   (Fig. S9e), which is consistent with previous literature (1 -4 mmol m -3 day -1 , with L. polyedra concentrations of 10 5 cells per liter, 10, 11).The estimated nitrate loss (  ) was simply obtained by integrating   in time with a zero initial value, and had the same error bars as  , .Generally, the error bars for depth integrals,   , can be expressed as: where  is the number of density layers,  is the depth difference between each layer, and  is the error bar of the integrated variable.
In all, these approaches allow us to accurately decouple vertical and horizontal advection signals from biologically driven changes of the nitrate field -even though we had only a single mooring and so could not directly estimate horizontal spatial gradients.The results produce features in the data that are consistent with other independent data (PAR, Chl-a, turbidity, etc.), further supporting the feasibility of this T/S binning method; this method could potentially be applied to other research problems in similar situations.
Finally, to emphasize how important it is to remove horizontal advection when estimating nitrate uptake rate, we compare the depth-integrated nitrate loss calculated from the raw observational nitrate data to the "de-advected" nitrate data (  ) (Fig. S10).This comparison shows clearly that there are large excursions in nitrate concentration in the raw nitrate data (the blue curve), which are not seen in the "de-advected" data (red line).Our calculations show that these fluctuations were driven by horizontal advection of horizontal, along-isopycnal gradients of nitrate.