Small sinking particles control anammox rates in the Peruvian oxygen minimum zone

Anaerobic oxidation of ammonium (anammox) in oxygen minimum zones (OMZs) is a major pathway of oceanic nitrogen loss. Ammonium released from sinking particles has been suggested to fuel this process. During cruises to the Peruvian OMZ in April–June 2017 we found that anammox rates are strongly correlated with the volume of small particles (128–512 µm), even though anammox bacteria were not directly associated with particles. This suggests that the relationship between anammox rates and particles is related to the ammonium released from particles by remineralization. To investigate this, ammonium release from particles was modelled and theoretical encounters of free-living anammox bacteria with ammonium in the particle boundary layer were calculated. These results indicated that small sinking particles could be responsible for ~75% of ammonium release in anoxic waters and that free-living anammox bacteria frequently encounter ammonium in the vicinity of smaller particles. This indicates a so far underestimated role of abundant, slow-sinking small particles in controlling oceanic nutrient budgets, and furthermore implies that observations of the volume of small particles could be used to estimate N-loss across large areas.


Supplementary
Supplementary Table 4: Correlations between anammox rates in the size fractionated water samples (Figure 3).

Compared fractions Spearman ρ (rho) Spearman p* Paired t-test T-test p df
Bulk and 1. 6 Table 5: Statistics for correlations between anammox rates and in situ particle abundances during April 2017 ( Figure 4). For each size class, a Spearman's correlation test was carried out using all anammox rates (if p > 0.05, the rate was assumed as zero) from anoxic water depths and corresponding particle abundances measured using a UVP5 (which was attached to the same CTD that was used to sample water for the rate determinations). A linear regression was carried out between particle abundances or biovolumes in each size class and anammox rates. The sample size was n = 24. In all cases, one anammox rate and associated particle abundances from the euphotic zone (24 m depth, station 341) was excluded as it was > 6 SD from the mean. All correlation tests were carried out using the linear model (lm) function of R.    Table 9: Correlations between organic nitrogen export out of the euphotic zone and integrated anammox rates -see Supplementary Figure 6. A Spearman's correlation test was carried out for nitrogen export estimated from UVP measurements and integrated anammox rates (if p > 0.05, the rate was assumed as zero) from anoxic water depths (where O 2 < 1.5 µM). A linear regression was carried out between particle abundances or biovolumes in each size class and anammox rates. The sample size was n= 24. Supplementary Figure 2: Total particle volume within the size classes and export production in the Peruvian upwelling system during June 2017 (M138). Average particle abundance profiles in a offshore and b onshore stations determined by UVP. Particles were sorted into four size classes with equivalent spherical diameters (ESD) of 128 -256 µm, 256 -512 µm, 512 -1024 µm, and 1 -2 mm. The shaded envelopes correspond to standard deviation in abundances from 20 offshore and 15 onshore stations. They were cropped at 2 mm 3 L -1 . The base of the euphotic zone (PAR < 1%, based on satellite data) and the oxic-anoxic interface (where O 2 dropped to < 1 µM) are indicated. c Average chlorophyll a concentrations from satellite data. d Export production estimated from satellite products (pseudocolor map) and from particle abundances at the base of the euphotic zone (circles, see text). The white line is the 600 m isobath on which basis the onshore and offshore stations were separated. Stations where anammox incubations were performed are circled in red. Supplementary Figure 4: Ratio between anammox and denitrification rates determined on both M136 and M138. a Volumetric rates, b integrated rates at each station. The dashed line represents a 1:1 relationship, therefore any value below this line is representative of a depth or station where anammox rates were higher than denitrification rates. Figure 5: Comparison of organic nitrogen export from the photic zone estimated from satellite imaging and UVP profiles. Organic nitrogen export was estimated from the UVP based on particle concentration at the base of the euphotic zone, measured size specific carbon content of the particles and calculated particle settling velocity (black). Organic nitrogen export estimates from satellite data were calculated using two different methods, firstly based on calculating a primary production export-ratio from temperature and chlorophyll concentrations 4 (red) and secondly based on a regional regional estimate of export to primary production ratio (blue). Correlations for all individual data points are given in a and distribution of values is shown in b, histogram showing the distribution of organic N export based on the three different methods, organic nitrogen export estimates were binned into classes of 0.17 mmol N m -2 d -1 and counts refers to the number of depths where that value was estimated. Red bars are export based on satellite data (regional estimate), blue bars are export based on satellite data (global estimate) and black bars are export based on UVP data. Onshore and offshore are combined in panel b.

Supplementary
Supplementary Figure 6: Comparison between integrated anammox rates and N-loss rates and the export of organic nitrogen out of the euphotic zone ("export"). Three methods were used to determine export, one based on particle abundances at the base of the euphotic zone from UVP data, measured size specific carbon content of the particles and estimated excess densities 5 (a, c), and two based on calculating a primary production export-ratio from temperature and chlorophyll concentrations as determined by satellite imagery, either based on a global estimate of export to primary production ratio 4 (b, e), or a regional estimate (c, e). Integrated anammox rates (in N 2 ) and shown in the top row, and integrated N loss (i.e. anammox and denitrification) in the second row. Lines indicate the linear regression, for statistics see Supplementary Table 9. Figure 7: Ratio of anammox rates in bulk water with rates in 1.6 µm filtered (panel a) and 10 µm filtered water (panel b) from the same sampling bottle. As there is no pattern towards increased rates in the bulk water, we assume that anammox cells are free living (see also Figure 3 in the main document), The standard error of the slope for each data point is shown in Supplementary Table 4. 23 samples did not have a significant anammox rate in any of the fractions and are therefore on top of each other at the zero point in this figure.

Supplementary
Supplementary Figure 8: Relationship between volumetric anammox rates and particle biovolumes from different size classes in the Peruvian upwelling system during April 2017. Particles throughout the water column were quantified with a UVP and binned into four size classes (each depicted in a different color). Anammox rates were determined from the slope of 29 N 2 production over time in anoxic incubations after addition of 15 NO 2 -(taking into consideration any contribution to 29 N production from denitrification if a denitrification rate with p < 0.05 could be detected). For all particle size classes, there was a significant positive correlation (Spearman's rank correlation; p = < 0.05 between particle biovolume and anammox rates from the anoxic part of the OMZ (O 2 < 1.5 µM in situ). The line is the linear regression from which the slope and R 2 was calculated (see Supplementary Table 5 for all relevant statistics). One outlier was removed from the figures and correlation as it was more than 6 SD from the mean in the smallest size fraction and it was the only sample from the euphotic zone (24.3 m depth). Samples from depths where in situ oxygen concentrations were > 1.5 µM were excluded. Supplementary Figure 10: Sensitivity of the trend of Organic N Export, N Release and Anammox Particle Encounters along increasing particle sizes based on variations of the particle settling velocity. For each UVP depth, the slope between the four particle sizes was calculated and plotted against the total particle abundance. Slopes below zero indicate that small particles dominate the respective process and -vice versa -slopes above zero indicate that large particles dominate. a -c Reference cases based on Stoke's Law with a variable porosity as described in equations (2)(3)(4)(5). d -f Power law approach for the settling velocity after Kriest 6 : U= 132r 0.62 , which is based on a compilation of various datasets. The results indicate some variability based on the applied relationship, however, overall the results are in the same range and the dominance of small particles for the respective processes persists. Figure 11: Organic N Export and N release for variable C:N ratios. a+e Reference case for the C:N ratio applied in this study. b+f An increased C:N ratio does not change the overall trend but reduces the total organic N export and N release. c+g Size-class dependent C:N ratio, where small particles have a C:N ratio of 12 and larger particles a C:N ratio of 6.6. d+h Hypothetical scaling relationship for the C:N ratio required to compensate for the effect of small particles. The size dependent C:N ratios are 29, 18, 12, 8 and 5 for 128 µm -256 µm, 256 µm -512 µm, 512 µm -1024 µm and 1024 µm -2048 µm, respectively. Figure 12: Schematic illustration indicating the applied boundary layer theory. The diffusive boundary layer thickness of a settling particle is determined by, both, the ratio of inertial to viscous forces, represented by the Reynolds-number, and the ratio of viscous to diffusive forces, represented by the Schmidt number (see main text). The Reynolds number determines the thickness of the flow boundary layer, while the Schmidt number, determines the amount of solutes diffusing from the particle interstitials into the surrounding seawater. For that reason, the diffusive boundary layer thickness relative to the particle size δ/r depends on both the Reynolds number and Schmidt number, this relationship is represented by the Sherwood (Sh) number (see main text). a Applying previously derived scaling laws reveals that for small particles viscous forces are more dominant resulting in a wider DBL relative to the particle size. b In case the diffusive boundary layer is hypothesized to be constant relative to the particle size (Sh = const), larger particles will have a strongly increased diffusive boundary layer thickness compared to small particles. Please notice that horizontal lines of the same color indicate the same length. Figure 13: Encounter rates of anammox bacteria with the diffusive boundary layer of a sinking particle in response to changing diffusion coefficients. a Assuming a realistic parametrization for the diffusion coefficient, anammox bacteria predominantly encounter small particles, resulting from the large amounts of particles and a larger boundary layer thickness of the particle relative to their size (see also Supplementary Figure 12). b The importance of small particles is less prominent with smaller diffusion coefficients or increased viscosities which under in situ conditions could result from colder water temperatures. c Assuming a hypothetical, constant ratio of boundary layer thickness and particle size (Sh = const), anammox bacteria would roughly encounter as many small particles as they would encounter larger particles. Figure 14: Encounter rates of sinking particles with anammox bacteria, assuming motility or non-motility. Within the encounter model, we assume that anammox bacteria are non-motile to estimate a conservative range of encounter rates. However, recently it was shown that encounter rates increase if bacteria are motile 7 . As motility values for anammox bacteria are not described, we adapted a motility model based on an intermediate motility of 1x10 -9 m 2 s -1 (ref. 7 , eq. 3, green violin plots). Inclusion of motility increases encounter rates by 406 times for smaller particles and 179 times for larger particles. The ratio of encounter with smaller particles compared to encounters with larger particles increased from 2.6 to 3.8. The outlines of the violin plots depict the kernel density estimation of the data points shown. The black line and red line indicate mean and median respectively. We estimated the size specific remineralization length-scale based on the exponential decay of the carbon fluxes (eq. 6). We used three reference depths, namely euphotic depth (z eup ), upper boundary of the OMZ (z omz ) and the core of the OMZ (z omz-50 ) which is assumed to start 50 m below the upper boundary of the OMZ. All profiles for the offshore stations (> 600 m water-depth) and onshore stations (< 600 m water-depth) were aligned to the reference depths and subsequently averaged. The information is summarized in Supplementary Table 6. Figure 16: Carbon content determined for 66 aggregates sampled using the Marine Snow Catcher. This information was used to rescale the volume-carbon relationship described by Alldredge 8 to match the particles we found in the Peruvian OMZ. Subsequently, the volume-carbon relationship was used to estimate the carbon export for different particle size classes based on their volume (see methods). Boxplots depict the 25-75 % quantile range, with the center line depicting the median (50 % quantile); whiskers encompass data points within 1.5 times the interquartile range. Figure 17: Encounter rates of anammox bacteria with particles correlated with anammox rates. Encounter rates were estimated from particle abundance in four different size classes (UVP data) combined with cell numbers determined by Hamersley et al. 9 , statistics of the linear regressions can be found in Supplementary Table 10.

Comparison of UVP-based and satellite-based estimates of organic nitrogen export from the euphotic zone
Organic nitrogen export was estimated from UVP measurements based on particle concentrations at the base of the euphotic zone or based on satellite imagery (see methods and Supplementary Figure 5).
UVP-based estimates were highly variable, with a median of 16.6 ± 10.9 mmol N m -2 day -1 (± is median absolute deviation, see Supplementary Figure 5) , which is line with the large variations in particle numbers ( Figure 2). The export estimates based on satellite imagery were less variable with a median of 3.7 ± 2.1 mmol N m -2 day -1 . Satellite-based export production and UVP based organic nitrogen export were significantly correlated (Spearman's rank test, p-value < 2.2 ・10 -16 , slope = 1.1), however median values were significantly different from each other (Wilcoxen signed rank test, p-value < 2.2・10 -16 ). This discrepancy likely results from errors associated with the primary production to export ratio (pe-ratio) utilized in the satellite-based estimate. The pe-ratio is derived from a previously described relationship between export, chlorophyll concentrations and temperature 4 and is based on an average from a global dataset that does not take the regional and temporal variability of the phytoplankton community and heterotrophic activity into account (see ref. 4 ). Therefore, we tested the impact of using a constant pe-ratio of 0.4, which is based on measurements in the Peruvian upwelling system 4 .
This doubled the estimated median export production of 7.3 ± 3.0 mmol N m -2 day -1 , highlighting the sensitivity of satellite-based estimates to the poorly constrained pe-parameter. The satellite data based approach was also limited by data resolution (spatially: 4 km and temporally: 8 days) and cloud-coverage. These make it difficult to capture the dynamic conditions resulting from mesoscale eddies and lateral intrusions that are present in the Peruvian upwelling system 10-12 . In contrast, the UVP profiles represent a snapshot of the system at the time of sampling, have a higher spatial and temporal resolution and provide detailed information about particle abundances and sizes throughout the water column. Therefore we used the UVP estimates of organic N-export to compare to anammox rates, and to investigate the mechanisms that lead to the relationship between particle volume and anammox rates. It should be noted that for both satellite and UVP estimates organic N-export was high enough to support the measured anammox rates (Supplementary Figure 6).

Evidence that anammox bacteria are free-living in the bulk water incubations
Bulk water anammox rates were determined from a time series of measurements conducted in 12 mL glass vials that were incubated with 15 NO 2 -(see methods). Based on the in situ particle abundance determined from the UVP, it is likely that there were between zero and three particles in each of the 12 mL glass vials. If anammox bacteria were directly attached to the particles, then linear rates would not have been expected in the time series. Instead we would expect N 2 production from anammox in some vials and no N 2 production in other vials, rather than a steady increase over time. As we consistently found linear rates of anammox using this approach, despite the likely variation in particle numbers between individual vials, this suggests that anammox rates were not directly affected by the presence or absence of particles. The size fractionation experiments, which showed the same rates in all three treatments (all of which also had linear rates), confirm this result.

Carbon specific respiration rates in oxic waters and anoxic waters
Carbon-specific degradation rates were estimated using two different models both of which are based on the carbon flux attenuation calculated from the measured particle profiles (Supplementary Figure 15). The first model represents a simple exponential decay based on a first-order kinetic (equation (6) in methods) and the second model introduces an additional parameter that takes declining particle degradation with depth into account 3 , equation (7)  While sinking, particles are remineralized, but also undergo aggregation and fragmentation, which can impact carbon specific degradation rates when they are calculated from particle profiles. Fragmentation and disaggregation is driven by zooplankton feeding, microbial degradation and shear in the turbulent water column [13][14][15] . In effect this leads to the "disappearance" of larger particles in profiles (potentially leading to an overestimation of degradation rates in the larger size classes) and results in the formation of smaller particles masking remineralization (potentially leading to an underestimation of degradation in the smaller size classes). The difference in degradation rate estimated for the offshore stations by the two models might indicate that fragmentation was playing a role in the upper water column which declined with depth. This is more likely to be captured in the Extended degradation model as it estimates the overall shape of the flux attenuation profiles better (see methods, eq. 7).

Organic nitrogen export into the OMZ and remineralization
We tested how varying the parameters in the model impacted the estimates of organic nitrogen export into the OMZ and remineralization of organic matter within the OMZ (equations (1)(2)(3)(4) in the methods). The major assumptions included in these models used to make these estimates are the particle settling velocity, carbon specific degradation rate (see above), and the conversion factor that relates carbon to nitrogen export (i.e. the C:N ratio).
We compared two widely applied methods to determine settling velocity based on particle size. The first equation (equations (2)(3)(4)(5)  classes, respectively. The comparably higher export of organic N associated with smaller particles persisted using the power law approach but the ratio associated with smaller vs larger particles changed from to 1.2. The remineralization (or N release) ratio between the smaller to larger particles changed from 4.6 to 3.7.
We assumed a C:N ratio of 7.2 based on previous measurements from the same region 16 . However, preferential N-remineralization seems to occur in sinking organic matter 17,18 .
To test the effect that this would have on organic N-export estimates, we first changed the C:N ratio in the model to 9.1 (Supplementary Figure 11 panels b and f) based on measurements from deep sediment traps near Peru 18 . The increased C:N ratio decreased the organic N-export by 21 % for each of the four different size classes but did not affect the ratio between organic N-export by smaller and larger particles.
There is also ample evidence in the literature that the preferential remineralization and uptake of N that occurs during heterotrophic organic matter degradation leads to increasing C:N ratios as material ages [19][20][21] . Therefore, theoretically, preferential N-remineralization over time could lead to the smaller, slower sinking and older particles having a higher C:N ratio than that of larger, faster sinking, younger particles. It should be noted however, that most of these studies focus on the transformation of particulate to dissolved organic matter 20,22 , or ageing within the dissolved organic matter pool 19,21 , and that differences in C:N ratio between particle size classes are poorly constrained. However, to investigate the impact of such an effect we altered the C:N ratios in our model for the individual size classes (to 6.6 for the larger particles and 12 for the smaller particles; Supplementary Figure 11 panels c and g). This reduces the organic N-export by the smaller particles by 68 % and increases it by 7 % for the larger size classes. Based on this extreme case, the ratio of organic N-export between the smaller to larger particles changed from 3.4 to 1.9. In fact, in the model a hypothetical, size-dependent C:N scaling relationship of 396 r -0.61 would be needed to compensate for the trend between the particle size fractions. In this case, the smallest particle size class would have a 5.8 times higher C:N ratio than the largest size class.
Concerning the N-release estimates -if preferential N-remineralization of the organic matter is taking place while particles sink, then this implies that our N-release rates are underestimated. Preferential N-remineralization and the associated enhanced release of N into the water column would lead to a C:N ratio that gradually increases with depth. Therefore, the C and N mass fluxes would be decoupled leading to element specific degradation rates.
Implementation of C and N specific degradation rates within the models would require accurate measurements of C:N ratios through-out the water column and for different particle types.

Encounters between anammox bacteria and particles
Encounter rates were calculated based on particle size, abundance, settling velocity and boundary layer theory. We also assume that anammox bacteria are passive tracers in the water column (i.e non-motile) and that the bacteria only need to be near the boundary layer of a particle to utilize the released N. This differs from previous literature which has mainly focused on the colonization of particles by motile bacteria (e.g. refs. 7,23 ). The boundary layer associated with particles of different size classes and particle settling velocities are the key parameters determining the results from the model.
The changes in the flow-regime lead to large variability in the boundary layer thicknesses of small and large particles. Smaller and slower sinking particles have a wide boundary layer thickness relative to their size compared to larger and faster settling particles (Supplementary Figure 12). Therefore, anammox bacteria encountered the boundary layer of small particles 2.6 times more often. When we applied a different size-sinking velocity relationship (see above), the trend remains similar but anammox bacteria encounter small particles 1.5 times more often (Supplementary Figure 10). In case the boundary layer thickness would, hypothetically, be assumed to be constant relative to their size, anammox bacteria would encounter as many small particles as they encounter large particles (Supplementary Figure 13 c).
We also tested the impact of assuming that anammox bacteria are non-motile, using a particle encounter model developed for motile bacteria 7 . In comparison to the non-motile model, encounter rates within each size class increased by 2 to almost 3 orders of magnitude (Supplementary Figure 14). Encounter rates between anammox bacteria and smaller particles showed a stronger increase compared to those between anammox bacteria and larger particles, increasing the ratio from 2.6 (as found in the first iteration of the non-motile model; see above) to 3.8.

The influence of suspended particles
Suspended particles, which we assume are not imaged by the UVP and are therefore not included in the models discussed above, might have an additional influence on N-loss in the OMZ that is not addressed by our approach. Suspended particles can be present in the water column through 1) spontaneous assembly of free polymer chains into microgels (see ref. 24 ), however such transparent gels are formed of material which likely has low nitrogen content 24 ) and as such would contribute little to organic nitrogen release 2) through fragmentation or degradation of sinking particles to the extent that they stop sinking. In the second case, assuming steady state conditions, the presence of suspended particles is unlikely to have had a large impact on our organic N-export rates. This is supported by the similar export estimates we found using a UVP independent approach (see above) and the similarity of our fluxes to those in Cavan et al. 2 . The influence of suspended particles on N-release and encounter rates is harder to predict, as it requires a better understanding of the distribution of suspended particles, their nitrogen contents and their remineralization rates. However, it seems likely that remineralization would continue in very small particles that have become suspended. As the current approach does not take degradation and release of nitrogen from suspended particles into account it is likely that our modelled estimates are conservative.