Magnitude of nitrate turbulent diffusion in contrasting marine environments

Difficulties to quantify ocean turbulence have limited our knowledge about the magnitude and variability of nitrate turbulent diffusion, which constitutes one of the main processes responsible for the supply of nitrogen to phytoplankton inhabiting the euphotic zone. We use an extensive dataset of microturbulence observations collected in contrasting oceanic regions, to build a model for nitrate diffusion into the euphotic zone, and obtain the first global map for the distribution of this process. A model including two predictors (surface temperature and nitrate vertical gradient) explained 50% of the variance in the nitrate diffusive flux. This model was applied to climatological data to predict nitrate diffusion in oligotrophic mid and low latitude regions. Mean nitrate diffusion (~ 20 Tmol N y−1) was comparable to nitrate entrainment due to seasonal mixed-layer deepening between 40°N–40ºS, and to the sum of global estimates of nitrogen fixation, fluvial fluxes and atmospheric deposition. These results indicate that nitrate diffusion represents one of the major sources of new nitrogen into the surface ocean in these regions.

"Follow the nitrogen", the mantra proposed by Capone et al. 1 for searching for life on Mars, summarizes the relevance of this element for life. By limiting both land and marine primary productivity, this element controls the flow of energy through the ecosystem, and the atmospheric CO 2 uptake by photosynthetic organisms. Marine photosynthesis can be supported by nitrogen inputs supplied from outside the photic layer (named new production, mainly in the form of nitrate), or by nitrogen remineralized within the photic layer (regenerated production) 2 . In steady state, when large temporal and spatial scales are considered, and for fixed stoichiometry of organic matter, new production is equivalent to the amount of organic carbon that can be exported to the deep ocean, where it remains isolated over the time scale of deep ocean circulation (100 s-1000 s years) 3 . The export efficiency of organic carbon production in the ocean can be characterized in terms of the f-ratio, defined as the ratio of new production to the sum of new plus regenerated production (total production) 4 .
Nitrate turbulent diffusion has been traditionally considered one of the main mechanisms by which nitrogen is supplied to phytoplankton, particularly in permanently stratified ocean environments such as tropical and subtropical regions 5 , and in temperate areas during summer stratification 6,7 . Quantifying nitrate diffusion requires estimates of vertical turbulent diffusivity (Kz), which can be derived from observations of microstructure, tracer release experiments and acoustic measurements of flow velocities 8 . These experimental approximations present methodological difficulties that have historically limited our knowledge about the magnitude and variability of turbulent mixing in the ocean. Alternatively, constant values of Kz have been used to estimate nitrate turbulent diffusion 9 . However, the increased availability of commercial microturbulence profilers has enlarged the number of observations, which have revealed a large temporal and spatial variability of Kz in the upper layer of the ocean [10][11][12][13][14][15][16][17] .
The limited number of studies that have simultaneously quantified the relevance of different new nitrogen supply mechanisms, by using observations at different spatial scales, have reported contradictory results about the role of nitrate diffusion. With measurements across the tropical and subtropical Atlantic, Pacific and Indian oceans, Fernández-Castro et al. 18 showed that nitrate diffusion (171 ± 190 mmolm −2 d −1 ) clearly dominated over N 2 fixation (9.0 ± 9.4 mmolm −2 d −1 ) at the time of sampling. When comparing with atmospheric nitrogen deposition and biological nitrogen fixation across a longitudinal section in the Mediterranean Sea, the contribution of nitrate diffusion to new production ranged from 0.1 to 25% 19 . In the northeast subtropical Atlantic, Painter et al. 20 reported nitrate diffusive fluxes comparable to N 2 fixation rates. Finally, according to Caffin et al. 21 N 2 fixation was the major source of new nitrogen (> 90%), compared to nitrate diffusion and atmospheric deposition, in the western tropical South Pacific. The bulk estimates of seawater nutrient concentration have been also frequently used as a proxy for nutrient availability in the euphotic zone 22,23 . However, nitrate concentrations and actual nitrate supply into the euphotic zone can be disconnected in near-steady-state systems, such as the subtropical gyres, where nitrate diffusion is slow, and nitrate concentrations are kept close to the detection limit due to phytoplankton uptake 24,25 .
We used a large dataset of microturbulence observations and a multivariable fractional polynomial method to investigate the relationship between nitrate turbulent diffusion and environmental variables that are routinely measured during oceanographic cruises. The observed relationship was used to build the first large-scale map of the distribution of the supply of nitrate into the euphotic zone through turbulent diffusion.

Results and discussion
Characteristics of the sampled regions. This Table 1  Most stations carried out in the Galician coastal upwelling were conducted inside three different Rías (Ría de Vigo, Ría de Pontevedra and Ría de A Coruña). The Rías are coastal bays influenced by seasonal wind-driven coastal upwelling which bring cold, nutrient-rich North Atlantic Central water to the surface 26 . In this region all stations considered were deeper than 40 m. One cruise sampled 10 stations located in the Antarctic Peninsula during the austral summer (COUPLING Jan-2010). Finally, the PERFILC cruise sampled 27 stations in the Bay of Biscay in July 2008. hHydrographic and microstructure turbulence profiles collected during both cruises were averaged arithmetically to obtain mean profiles which were assigned to the mean geographical location of all the stations conducted during each cruise (see Methods). Additional information about the sampling design of these cruises is described in Aranguren-Gassis et al. 27 , Mouriño-Carballido et al. 28 , Fernández-Castro et al. 18 , Mouriño-Carballido et al. 24 , Cermeño et al. 29 , Teira et al. 30 , Villamaña et al. 15 , Díaz et al. 31 , and Moreira-Coello et al. 32 .
This database covered a wide environmental gradient from oligotrophic to eutrophic conditions. Temporal variability due to seasonal changes in environmental forcing was partially captured in some regions www.nature.com/scientificreports/ (Northwestern Mediterranean Sea and Galician coastal upwelling) but mainly missed in others (tropical and subtropical Atlantic and Pacific oceans and the Antarctic Peninsula. The 70 stations sampled in the tropical and subtropical Atlantic and Pacific oceans were, on average, characterized by relatively warm surface temperature (26 ± 3 °C, mean ± SD), weak vertical diffusivity (0.03 ± 0.05 × 10 -3 m 2 s −1 ) and low nitrate diffusive supply (0.4 ± 1 mmol N m −2 d −1 ) ( Fig. 2 and Table 2). Accordingly, surface chlorophyll-a and photic layer depth-integrated net primary production were low (0.2 ± 0.1 mg m −3 and 202 ± 98 mg C m −2 d −1 , respectively). The Mediterranean Sea, was characterized by cooler surface waters (16 ± 4 °C), and relatively large diffusivity (2 ± 6 × 10 -3 m 2 s −1 ) and nitrate diffusive supply (39 ± 110 mmol N m −2 d −1 ). Surface chlorophyll-a and net primary production took intermediate values (0.9 ± 0.9 mg m −3 and 581 ± 454 mg C m −2 d −1 , respectively). The stations sampled in the Galician coastal were characterized, on average, by relatively cold surface water (16 ± 2 °C), and intermediate diffusivity (0.4 ± 0.5 × 10 -3 m 2 s −1 ) and nitrate diffusive supply (5.3 ± 7.1 mmol N m −2 d −1 ). However, due to the influence of coastal upwelling the total nitrate supply (including both diffusive and advective processes) was significantly higher (31 ± 42 mmol N m −2 d −1 ). As the result of the input of new nitrogen by both processes this system was characterized by high values of surface chlorophyll-a and net primary production (2.5 ± 2.9 mg m −3 and 2950 ± 2330 mg C m −2 d −1 , respectively). Finally, the 10 stations sampled in the Antarctic Peninsula were characterized by cold surface water (0.1 ± 1.0 °C), high diffusivity (66 ± 126 × 10 -3 m 2 s −1 ) and high nitrate diffusive supply (102 ± 165 mmol N m −2 d −1 ). However, due to the limitation of control factors other than nitrate 30 , these regions exhibited intermediate values of surface chlorophyll-a and relatively low net primary production (1.1 ± 0.8 mg m −3 and 153 ± 136 mg C m −2 d −1 , respectively).
Multivariable fractional polynomial models. At each station a set of twelve physical, chemical and biological variables regularly obtained during oceanographic cruises was selected as potential predictors of the nitrate supply by turbulent diffusion: sea surface temperature (SST), sea surface salinity (SSS), vertical stratification in the pycnocline (maxN 2 , defined as the maximum value of the squared buoyancy frequency), depth of the pycnocline (dmaxN 2 , depth of the maximum value of the squared buoyancy frequency), average stratification in the nitracline (avrN 2 , mean squared buoyancy frequency in this layer), depth of the mixed layer (MLD, calculated using a density difference criterion of 0.1 kg/m 3 with respect to the surface value), depth of the nitracline (nitraD, the shallowest depth were nitrate concentration was equal to 1 mmol m -3 ), nitrate gradient across the nitracline (grNO 3 ), surface nitrate concentration (sNO 3 ), depth (DCM) and value (maxCHL-a) of the chlorophyll-a maximum, and surface chlorophyll-a (sCHL-a).
In order to exclude cross-correlation between predictors and consider the possibility of non-linear relationships, a Multivariable Fractional Polynomial (MFP) method 33 was used to select the predictors in the four www.nature.com/scientificreports/ contrasting environments and the complete dataset excluding those stations located in the Antarctic Peninsula (GM) ( Table 3), where surface temperature was below 15 °C (Fig. 3). In a first exploratory step we allowed the models to introduce any of the twelve potential predictors (exploratory models). In tropical and subtropical regions a model including four predictors (grNO 3 , SSS, avrN 2 and sNO 3 ) explained 66% of the variance (Adj-R 2 ) in nitrate diffusive fluxes, whereas in the Mediterranean and the Galician coastal upwelling two predictors (SST and sCHL-a, and grNO 3 and MLD; respectively) explained more than 55% of the variance. In the Antarctic Peninsula, a model including only SST explained 75% of the variance. Finally, six predictors were selected for the GM global model (grNO 3 , sNO 3 , avrN 2 , SST, dmaxN 2 , maxN 2 ; Adj-R 2 = 0.63).
In order to build more simple models by limiting the number of potential predictors, we investigated the correlation between predictors by using a correlogram (Fig. 4). The cluster analysis created using a dendrogram based on the correlation matrix between the predictors showed that variables could be clustered into three groups. www.nature.com/scientificreports/ Among these clusters we selected three variables (SST, grNO 3 , and sCHL-a) to build the 3-variable models. The three selected predictors were chosen because they are routinely measured during oceanographic cruises, and also derived from satellite observations or available on global oceanographic databases. This was not necessary for the Mediterranean and the Antarctic, where the exploratory models already selected a combination of these predictors. For tropical and subtropical regions the algorithm retained one variable (grNO 3 ), which explained about half of the variance (27%) compared to the model with all the potential predictors included. In the Galician coastal upwelling only grNO 3 was retained, but the decrease in explained variance was smaller (Adj-R 2 = 0.49). The global model maintained its accuracy very close (Adj-R 2 = 0.50) after decreasing the number of predictors to two variables (SST and grNO 3 ). Several factors could explain the variability in the variance explained by the different regional 3-variable models. First, the dataset collected in the Antarctic and, specially, in the Mediterranean included contrasting hydrographic conditions and a wider range of nitrate diffusive fluxes than the other regions ( Fig. 2). Second, in the Galician upwelling nitrate turbulent diffusion represents, on average, a minor source of new nitrogen into the euphotic zone compared to the nitrate flux driven by vertical advection of deeper waters through upwelling 32 . Moreover, it is noteworthy that the explained variance in tropical and subtropical regions decreased significantly from the exploratory (66%) to the 3-variable models (27%), whereas the reduction in the number of predictors had a limited effect in the other regions. This could be related to the fact that the observations included in our database barely cover seasonal changes in hydrographic conditions in this oceanic domain. The dependence of nitrate diffusion on the different explanatory variables could be interpreted in terms of the underlying physical and ecological processes. All models in which SST was included exhibited negative linear relationships between nitrate diffusive fluxes and this variable (Fig. 3). The vertical nitrate gradient exhibited a positive linear correlation with nitrate diffusive flux in the model built for the Galician coastal upwelling. For tropical and subtropical regions and the global model, the relationship with the vertical nitrate gradient was also positive but logarithmic. Finally, surface chlorophyll-a showed a linear or nearly linear positive relationship with nitrate diffusive fluxes in the Mediterranean. Sea surface temperature is related to vertical stratification, as stratification usually increases with higher surface temperature. As expected, a statistically significant positive correlation was found between the two variables in our dataset (Fig. 4). In turn, assuming no systematic compensatory changes by increased dissipation rates of turbulent kinetic energy or enhanced vertical nitrate gradient, vertical diffusivity, and therefore nitrate diffusive fluxes, decreases when stratification increases (see Methods). This could explain the negative relationship observed between surface temperature and nitrate diffusive flux for all models where this predictor was included. The relationship between the vertical nitrate gradient and nitrate diffusive flux was positive, as this term is directly involved in the calculation of diffusive fluxes (see Methods). Finally, when included, the relationship between surface chlorophyll-a and nitrate diffusive fluxes was positive, as this predictor is an estimator of phytoplankton biomass, and therefore reflects the response of these organisms to nitrate availability. Table 2. Mean ± standard deviation of selected variables. Sea surface temperature (SST), sea surface salinity (SSS), vertical stratification in the pycnocline (maxN 2 ), depth of the pycnocline (dmaxN 2 ), average stratification in the nitracline (avrN 2 ), depth of the mixed layer (MLD), vertical turbulent diffusivity (K), depth of the nitracline (nitraD), nitrate gradient across the nitracline (grNO 3 ), surface nitrate concentration (sNO 3 ), diffusive nitrate supply (dNO 3 ), total nitrate supply (including diffusive and advective fluxes, TNO 3 ), depth (DCM) and value (maxCHL-a) of the maximum chlorophyll-a, surface chlorophyll-a (sCHL-a), and photic layer depth-integrated net primary production (PP), computed for the tropical and subtropical regions (T), the Mediterranean Sea (M), the Galician coastal upwelling (G), and the Antarctic Peninsula (A). A nonparametric one-way ANOVA (Kruskal-Wallis, KW) was performed to test the null hypothesis that independent groups come from the same distribution. The Bonferroni multiple comparison test was applied a posteriori to analyze the differences between every pair of groups.  5). As expected based on the regional differences observed in our estimates of nitrate diffusive fluxes derived from observations (Fig. 2 Table 4. Because SST and grNO3 data obtained from WOA09 include observations collected at different seasons, mean nitrate fluxes computed for the biogeographical provinces can be considered as spatio-temporal means. In general, a good correspondence was found between our model estimates and the limited number of studies that so far have quantified the magnitude of nitrate diffusive fluxes derived from observations of microstructure, tracer release experiments and high-resolution horizontal currents (Fig. 5 and Figure S2). However, a meaningful direct comparison between our mean estimates of nitrate diffusive fluxes with those derived from local Table 3. Regression equations obtained by the MFP method in each domain (D) by using exploratory and 3-variable models. T indicates tropical and subtropical regions, M Mediterranean Sea, G Galician coastal upwelling, A Antarctic Peninsula, and GM global model built with the complete dataset except those collected in the Antarctic Peninsula. Exploratory models used all the predictors, and 3-variable models only used sea surface temperature (SST), surface chlorophyll-a (sCHL-a) or vertical nitrate gradient (grNO 3 ). Adj-R 2 is adjusted squared coefficient of determination, and Std Error standard error. The number of stations used to build each model is indicated in Table 1. www.nature.com/scientificreports/ observations, whose representativeness is limited in space and time, is very difficult to achieve. Turbulent mixing is highly intermittent in time and patchy in space due to its chaotic nature, and the multiple scales of variability of the forcing mechanisms. In the upper ocean, energy dissipation and mixing are linked to the exchanges of momentum and buoyancy with the atmosphere, which vary between seasons, day and night cycles and also following the synoptic weather patterns 35 . In the pycnocline, mixing also depends on highly unpredictable processes, such as the generation of shear by internal waves 36,37 . In coastal systems, tidal forcing can be a predictable source of mixing 38 , but it can also lead to complex dissipation patterns as a result of the production of non-linear internal waves 11,15 , and the interaction with wind-driven circulation 39 . Furthermore, the magnitude of nitrate fluxes is very sensitive to the different methodological choices made for the calculation, in particular to the depth-interval. Different depth criteria used for the estimates reported in Fig. 5 and Table 5 include the base of the euphotic zone, the base of the mixed layer, and also the nitracline. In fact, significant differences in nitrate diffusive fluxes computed in part of our dataset (the TRYNITROP cruise carried out in the subtropical Atlantic) were noticed depending on whether the base of the euphotic zone 28 or the nitracline 18 were chosen as depth-range for the calculation. A few studies have presented similar attempts to model nitrate availability in the surface waters by using different oceanographic predictors. In an early effort, Kamykowski and Zentara 40 investigated the capability to predict nutrient concentrations from temperature or density using a global dataset. With satellite SST observations, Switzer et al. 41 generated an index of nitrate availability in the surface waters of the global ocean. Steinhoff et al. 42 proposed a proxy for nitrate in the North Atlantic Ocean, using multiple linear regression and observed nitrate and SST, and model-based MLD. Arteaga et al. 43 developed a method for estimating global monthly mean surface nitrate by using local multiple linear regressions and including satellite SST and sCHL-a, and modelled MLD. Finally, Liang et al. 44 used multivariate empirical orthogonal functions constructed from modeled nitrate, salinity, and potential temperature fields to build nitrate maps in the Southern Ocean. All these studies have assumed implicitly that the nitrate concentration alone is sufficient to characterize the variability in the supply  Table 3 Contribution of nitrate diffusion to primary production and export. We next investigated how our estimates fit in the nitrogen budget for tropical and subtropical regions, where most of our observations were collected. We first computed the spatial mean of nitrate diffusive fluxes between 40° N and 40° S using the global model. In order to limit our calculation to open ocean oligotrophic waters, we selected those regions where   Table 3. The colored dots represent nitrate diffusive fluxes previously reported in the literature and described in Table 5. The white color indicates those regions where data were outside of the range covered by the dataset used to build the prediction model.  Figure S3). These results indicate that, together with mixed-layer entrainment, nitrate diffusive flux represents a major source of new nitrogen into the euphotic zone in these regions. Assuming Redfield stoichiometry (C:N = 6.6), this flux may support at least 13-32% of the global export production of organic carbon, which is estimated to be 5-12 PgC y −149 . The sum of our estimates of nitrate diffusion and nitrate entrainment for the region 40° N-40° S, and the global estimates of nitrogen fixation, fluvial nitrogen fluxes and atmospheric nitrogen deposition (~ 46-50 Tmol N y −1 ) represented ~ 37-40% (f-ratio = 0.37-0.40) of the regional mean net primary production derived from satellite data (183 ± 108 Tmol N y −1 , see Methods). This is higher than previous estimates of the f-ratio in these regions, which according to several studies using sediment traps and modelling approaches fall within the range 0.1-0.2 2,4,50-52 . However, it is consistent with higher estimates of f-ratio derived for the Sargasso Sea from considering most identified physical and biological nitrogen supply mechanisms (see Table 2 in Lipschultz et al. 5 ).

Conclusions
Because nitrate availability is not necessarily correlated with nitrate concentration, it is crucial to quantify turbulent diffusion, traditionally considered one of the main mechanisms supplying new nitrogen to the phytoplankton cells inhabiting the surface layers of the ocean. By using a large dataset of microturbulence observations collected in contrasting marine environments and a multivariable fractional polynomial method, we built regional and global models that included a maximum of only three variables as predictors (surface temperature, SST; nitrate vertical gradient, grNO 3 , and surface chlorophyll-a, sCHL-a). The global model, which included the first two predictors, explained 50% of the variance in nitrate turbulent diffusion. Regional models for the Northwestern Mediterranean (including SST and sCHL-a) and for the Antarctic Peninsula (including SST) explained, respectively, 63% and 75% of the variance. In contrast, regional models for tropical and subtropical regions and for the Galician coastal upwelling (both including grNO 3 ) explained, respectively, 27% and 49% of the variance. The global model applied to temperature and nitrate collected from databases provided the first large-scale map of nitrate supply into the euphotic zone through turbulent diffusion.
Given that chlorophyll is an estimator of phytoplankton biomass, which in turn depends on nutrient availability, including this variable as a predictor limits the application of these models for the study of the role that nitrate availability plays in controlling phytoplankton biomass and productivity. However, this approach can be very useful for constraining biogeochemical budgets. Mean nitrate diffusive flux (~ 20 Tmol N y -1 ) in oligotrophic tropical and subtropical regions was comparable to nitrate entrainment due to seasonal mixed-layer deepening between 40° N-40° S, and to the sum of global estimates of nitrogen fixation, fluvial fluxes and atmospheric deposition. This result confirms this process as one of the major mechanisms of new nitrogen supply into the surface www.nature.com/scientificreports/ oligotrophic ocean. The model presented here allows to quantify nitrate availability in the euphotic zone, globally and in selected regions, which is critical for understanding the export efficiency of organic carbon production, and its spatio-temporal variability, in response to natural and anthropogenic ocean change.

Methods
Sampling. During all cruises except COUPLING hydrographic properties and turbulent mixing were derived from a MSS microstructure profiler 53 . The MSS is equipped with a high-precision Conductivity-Temperature-Depth (CTD) probe, two microstructure shear sensors (type PNS06), and also a sensor to measure the horizontal acceleration of the profiler. The fluorometer included in the MSS profiler was calibrated with fluorometrically determined chlorophyll a concentrations ranging from 0.03 to 8.60 mg m −3 (Chl a = 2.255 × fluorescence -0.527; R 2 = 0.859, number of samples (ns) = 134), obtained during 12 cruises (see below). Measurements of dissipation rates of turbulent kinetic energy (ε) were conducted to the bottom, or to 243 ± 23 m over deep waters ( Table 1). The MSS profiler was balanced to have negative buoyancy and a sinking velocity of ~ 0.4 to  www.nature.com/scientificreports/ 0.7 m s −1 . The frequency of data sampling was 1024 Hz. The sensitivity of the shear sensors was checked after each use. Data processing and calculation of dissipation rates of ε was carried out using the commercial software MSSpro 53 . During COUPLING turbulent mixing was derived from a TurboMAP microstructure profiler 54 . Data processing and calculation of dissipation rates of turbulent kinetic energy ε were carried out with the commercial software of the TurboMAP, as described in Sangrá et al. 55 . Microstructure turbulence profiles used for computing nitrate fluxes at each station were always deployed successively. Sets include 2-11 in the tropical and subtropical regions, 6-94 in the Mediterranean, 2-402 in the Galician coastal upwelling, and 2-3 in the Antarctic Peninsula. Due to significant turbulence generation close to the ship, only the data below 5 m (HERCULES2, HERCULES3, DISTRAL, REIMAGE, STRAMIX, ASIMUTH, CHAOS, and NICANOR) and 10 m (CARPOS, TRYNITROP, COUPLING, MALASPINA, PERFILC, PERFILM, FAMOSO1, FAMOSO2, FAMOSO3) were considered reliable. Data were depth-averaged by calculating mean values within 1 m bins.
The squared buoyancy frequency (N 2 ) was computed from the CTD profiles according to the equation: where g is the acceleration due to gravity (9.8 m s −2 ), ρ w is a reference seawater density (1025 kg m −3 ), and ∂ρ/∂z is the vertical potential density gradient. Vertical diffusivity (K z ) was estimated as: where Ŵ is the mixing efficiency, here considered as 0.2 56 for all cruises except MALASPINA, TRYNITROP and CARPOS, which included stations with favorable conditions for double diffusion through salt-fingers. During these cruises vertical diffusivity including mechanical turbulence and the effect of salt-fingers mixing was calculated according to St. Laurent and Schmitt 57 (see details in Fernández-Castro et al. 18 ).
Nitrate supply. Samples for the determination of nitrate (NO 3 ) + nitrite (NO 2 ) were collected from 5 ± 2 (Galician coastal upwelling), 7 ± 1 (Mediterranean), 11 ± 2 (Antarctic Peninsula) and 11 ± 2 (tropical and subtropical regions) different depths in rinsed polyethylene tubes and stored frozen at − 20 °C until analysis, according to standard methods using the automated colorimetric technique 58 . Analyses were performed on land except during the MALASPINA expedition where the samples were analyzed on board. When nitrate concentrations were not available, during STRAMIX and one sampling during the NICANOR cruises, concentration values were obtained from a nitrate-density (σ t ) relationship built by using all samples collected during CHAOS (ns = 624) and NICANOR (ns = 52), respectively. In CHAOS the nitrate-density relationship (NO 3 = 11.93 × σ t − 310.60; Adj-R 2 = 0.92; p < 0.001) was valid for the density range between 25.9 and 27.2 kg m −3 . During NICANOR the relationship showed a linear behavior (NO 3 = 9.78 × σt − 256.38; Adj-R 2 = 0.87; p < 0.001) for density ranging between 26.1 and 27.1 kg m −332 . Nitrate concentration profiles for PERFILM and PERFILC cruises were obtained from the World Ocean Atlas 2009 (WOA09,. http:// www. nodc. noaa. gov/) as the closest available bin to the mean geographical location of each cruise (40.770° N-2.675° E and 43.865° N-2.173° W, respectively), and considering the mean values for the month when each cruise was conducted. The WOA09 was also used to obtain nitrate profiles at four stations during MALASPINA where samples for the determination of nutrients were not collected 18 .
Vertical diffusive fluxes of nitrate into the euphotic zone were calculated following the Fick's law as: where NO 3 is the nitrate vertical gradient obtained by linear fitting of nitrate concentrations in the nitracline, determined as a region of approximately maximum and constant gradient, and K z is the mean turbulent mixing in the same depth interval. In the Galician coastal upwelling, nitrate diffusive fluxes were estimated between 10 and 40 m depth using the same procedure. The total nitrate supply in the Galician Rías was computed as the sum of nitrate vertical diffusion plus nitrate vertical advection due to coastal upwelling. A simplified estimate of nitrate supply through vertical advection due to upwelling was computed considering the Galician Rías as single boxes divided into two layers, the deeper one influenced by upwelled inflowing waters and the surface layer dominated by the outgoing flow. Assuming that the bottom layer volume is conservative and stationary, the vertical advective flux (Q Z , m 3 s −1 ), would be equivalent to the incoming bottom flux (Q B , m 3 s −1 ), computed as the product of the upwelling index (I W , m 3 s −1 km −1 ) and the lengths of the mouth of the Rías (ca. 10-11.5 km). I W was averaged by calculating the mean value over the three-day period before each cruise from wind data recorded by meteorological buoys located in Cabo Vilano (HERCULES, NICANOR) and Cabo Silleiro (DISTRAL-REIMAGE, ASIMUTH, CHAOS, STRA-MIX, REIMAGE), or modeled by the Fleet Numerical Meteorology and Oceanography Center (FNMOC) model when buoy data were not available (http:// www. indic edeafl oram iento. ieo. es). Finally, the transport of nitrate into the euphotic zone through vertical advection was computed as: where A basin is the surface area of the Galician Rías, Q Z is the vertical advective flux, and [NO 3 ] D is the mean nitrate concentration at the base of the euphotic layer. A basin is 141 km 2 for Ría de Pontevedra (ASIMUTH), 174 www.nature.com/scientificreports/ km 2 for Ría de Vigo (CHAOS, ASIMUTH, DISTRAL-REIMAGE, STRAMIX), and 145 km 2 for Ría de A Coruña (HERCULES, NICANOR). Total nitrate supply included diffusive and advective vertical processes in the Galician costal upwelling, and only nitrate turbulent diffusion in the other regions. These results were only used to describe the magnitude of total nitrate supply in the four regions, whereas only nitrate diffusion was used as the response variable in the prediction models.
Chlorophyll-a and primary production. Samples for the determination of chlorophyll-a were collected at four to eight depths during CARPOS (number of stations (n) = 8), TRYNITROP (n = 18), MALASPINA (n = 47), FAMOSO (n = 18), NICANOR (n = 13), HERCULESIII (n = 1), DISTRAL (n = 9), CHAOS (n = 2), and COUPLING (n = 9) cruises. 50-250 mL of seawater were filtered through 0.2 μm pore-size polycarbonate or Whatman GF/F (FAMOSO) filters which were later frozen at −20 °C until analysis. The fluorescence emitted by the chlorophyll-a was measured from pigments extracted in 90% acetone overnight. The fluorescence due to chlorophyll-a was measured using a Turner TD-700 fluorometer previously calibrated with pure chlorophyll-a. For those cruises where samples for the determination of chlorophyll-a were not collected, chlorophyll-a was derived from the calibrated fluorescence sensor included in the MSS profiler.
Samples for the determination of net primary production were collected at selected stations and several depths during TRINITROP (n = 18), MALAPINA (n = 37), FAMOSO (n = 19), NICANOR (n = 13), HERCU-LESIII (n = 1), DISTRAL (n = 8) and COUPLING (n = 3) cruises. During TRYNITROP, MALAPINA, NICANOR, HERCULESIII, and DISTRAL water samples from two to seven depths were collected for the determination of primary production with the 14 C-uptake technique during on-deck incubations. For each depth, 75-mL (250ml MALAPINA, NICANOR & HERCULESIII) polystyrene bottles (3 light and 1 dark bottles, or 2 light and 1 dark bottles for FAMOSO, NICANOR & HERCULESIII) were filled with seawater just before dawn and spiked with 2-15 μCi (50-100 μCi MALASPINA) of NaH 14 CO 3 . Samples were incubated during 2 h (DISTRAL), 24 h (MALASPINA, FAMOSO, NICANOR & HERCULESIII) or from dawn to dusk (TRYNITROP) in flow-through incubators provided with neutral density and blue (Lee Mist Blue) filters that simulated the PAR levels experienced by the phytoplankton in their natural location within the water column. The incubators were cooled with running seawater pumped from the surface or with water circulating through a refrigerator. At the end of the incubations, samples were filtered through 0.2 μm polycarbonate filters under low-vacuum pressure (< 100 mm Hg). Inorganic carbon on the filters was removed by exposing them to concentrated HCl fumes overnight. After removal of inorganic 14 C, filters were placed into scintillation vials to which 4 mL of scintillation cocktail were added. The radioactivity on each sample was determined on scintillation counters which used an internal standard for quenching correction. Detailed methods for the experiments conducted during these cruises is included in Mouriño-Carballido et al. 28 , Marañón et al. 59 , Estrada et al. 60 , Moreira-Coello et al. 32 , and Cermeño et al. 29 .
During COUPLING primary production was measured using the 13 C method 61 at two sampling depths (surface and at the depth of the deep chlorophyll maximum). Water samples were transferred to 2 L polycarbonate bottles, and after addition of NaH 13 CO 3 at about 10% of total inorganic carbon in the ambient water, the samples were incubated for about 12 h in a tank on deck. Initial and final particulate organic carbon, and particulate material used for isotope analysis were filtered through GF/F filters. They were frozen and stored at − 20 °C until analysis on land. Detailed methods are included in Teira et al. 30 .
For those experiments where incubations lasted less than 24 h, net primary production was computed assuming the ratio of phytoplankton respiration to gross photosynthesis (~ 20%; Geider 62 ), and the number of light hours for each sampling date. Depth-integrated values were computed by trapezoidal integration considering the rates measured at the different depths.

Model regression variables and statistical analysis. A Multivariable Fractional Polynomial (MFP)
method was used to select the predictors of nitrate turbulent diffusion in the four contrasting environments (tropical and subtropical Atlantic Ocean, Northwestern Mediterranean Sea, Galician coastal upwelling ecosystem and Antarctic Peninsula), and for the complete dataset. Regression analyses commonly assume that relationships are linear between the predictors and the response variable, or assume a log-transform without assessing these assumptions. When addressed, nonlinear relationships are often based on generalized additive models or use standard polynomials. However, these models neither can be expressed as an explicit equation nor handle general nonlinearity, respectively. The MFP method solves these problems and determines whether a predictor is important for the model, and its functional form. The algorithm use fractional polynomials based on box-cox transformation 63 of predictors. It uses the eight powers − 2, − 1, − 0.5, 0, 0.5, 1, 2 and 3 (with the 0 case corresponding to the natural log transform). A back-forward procedure is applied to include/exclude the terms of the model following the next steps: (I) backward elimination of predictors that are statistically insignificant, and (II) iterative examination of the scale of all continuous predictors. Therefore, we need two significance levels: α 1 , for the exclusion and inclusion of a predictor, and α 2 for the determination of significance of fractional transformation of continuous predictors. Successive cycles improve the accuracy of the model until two cycles converge. Goodness of fit was assessed via quantile-quantile plots ( Figure S1). All calculations were done using mfp package 33 in R 64 .
Global spatial distribution of nitrate turbulent diffusion. The large-scale spatial distribution of nitrate turbulent diffusion was computed by using the MFP model obtained by using the complete dataset except those stations collected in the Antarctic Peninsula, and including as predictors sea surface temperate and the vertical nitrate gradient across the nitracline derived from the WOA09. At each bin from WOA09 the vertical nitrate gradient was calculated as the maximal vertical gradient between 0 and 250 m. Regions where predictors www.nature.com/scientificreports/ data were outside the range covered by the dataset used to build the model were removed. Mean nitrate fluxes for most of the biogeographical provinces defined by Longhurst 34 were calculated using the 3-variable model built with the complete dataset except those stations collected in the Antarctic Peninsula (GM). Mean nitrate diffusive flux between 40° N and 40° S was calculated using the 3-variable model GM, and only considering those regions where sCHL provided by the GlobColour Project from the European Space Agency (ESA, hermes.acri.fr) and nitrate diffusive fluxes were lower than 0.5 mg m −3 and 1 mmol m −2 d −1 , respectively.
Contribution of nitrate turbulent diffusion to net primary production in the 40° N-40° S region. Spatial mean total net primary production for the region between 40º N and 40º S was derived from satellite monthly means (period 1998-2012) by using the model proposed by Uitz et al. 65 , which was depthintegrated down to the base of the photic layer, computed from the algorithm proposed by Morel et al. 66 . In order to convert net carbon production to nitrogen units we assumed mean stoichiometry relationship for carbon and nitrogen proposed by Galbraith and Martiny 67 . Only those regions where sCHL and nitrate diffusive fluxes were lower than 0.5 mg m −3 and 1 mmol m −2 d −1 , respectively, were considered.
Nitrate entrainment due to seasonal mixed-layer deepening. The entrainment flux of nitrate into the mixed-layer was calculated from monthly climatological nitrate from the World Ocean Atlas 2013 as: where dMLD/dt is the rate of deepening of the mixed layer between two consecutive months (i and i + 1, with i = 1,…,12); NO ml 3 and NO th 3 are the mean nitrate concentrations on month i in the first mixed layer (i.e. between the surface and MLD i ) and in the thermocline (between MLD i and MLD i+1 ), respectively, and is the Heaviside function, being = dMLD/dt when dMLD/dt > 0 and = 0 otherwise. The entrainment fluxes were calculated for each gridpoint of the WOA13 and for every pair of months and then time-averaged to obtain an annual mean estimate. Mixed layer depths where calculated from the monthly profiles of potential density as the depth were potential density exceeded the near-surface value by 0.1 kg m -3 .