Between February 2014 and December 2015, within the framework of the NICANOR (Nitrogen fixation and diffusive fluxes in the NW Iberian Peninsula) project, 10 one-day samplings were carried out on board B/O Lura at station 2 (43.42°N, 8.44°W; depth = 80 m) of the RADIALES time-series project (http://www.seriestemporales-ieo.com) (Figure 1). This station is located in the adjacent shelf off Ría de A Coruña (Golfo Ártabro, NW Iberian Peninsula). The samplings were designed to cover the main hydrographic stages in this system over a seasonal cycle, and spanned over a period of 20 months. To simplify the interpretation of the results, we classified the samplings into three conditions: five cruises were characterized by downwelling (February 2014, December 2014, May 2015, November 2015, and December 2015), three by upwelling (May 2014, April 2015, and June 2015), and another two by relaxation, i.e., intermediate conditions between upwelling and downwelling pulses (July 2015, and September 2015). This classification was based on the wind-driven upwelling index and the hydrographic conditions of the water column (see results). The upwelling index (I_W) was estimated by the Ekman transport (m³ s⁻¹ km⁻¹) and computed from wind data for the sector Rías Altas (44°N 9°W), near Ría de A Coruña, using the NAVGEM model of Fleet Numerical Meteorology and Oceanography Center (FNMOC) (http://www.indicedeafloramiento.ieo.es), averaged over the 3 days period before each cruise. On each sampling day, profiles of temperature, salinity and fluorescence were obtained with a CTD (Conductivity-Temperature-Depth) probe SBE25plus (SeaBird Electronics) attached to a rosette of Niskin bottles. Water samples were collected for the determination of inorganic nutrients, chlorophyll a, primary production, picoplankton composition, biological N₂ fixation, as well as for DNA extraction. In 7 out of the 10 surveys, measurements of dissipation rates of turbulent kinetic energy were conducted using a microstructure turbulence profiler (MSS, Prandke and Stips, 1998).

In order to facilitate the comparison between hydrographic conditions, we provide mean and standard deviation for the variables quantified under the relaxation condition, despite the low number of samplings (n = 2). We are aware that a higher number of samplings would be desirable to better characterize averaged physical, chemical and biological fields in this system, where variability happens at short temporal scales of a few days (Casas et al., 1997; Gilcoto et al., 2017).

Unfiltered samples for the determination of dissolved inorganic nutrients (NO3-, PO43-) were collected from 3 to 7 depths in rinsed polyethylene 15-mL tubes and stored frozen at −20°C, until further analysis by standard colorimetric methods on a segmented flow AutoAnalyzer 3 Bran Luebbe (Aminot and Kerouel, 2007). Nutrient analysis were performed at the facilities of the University of Oviedo (Spain). Data for nutrient concentrations were not available for December 2014 and April 2015. In these cases, nitrate concentration was computed from a nitrate-density (σ_t) relationship built by using all samples (n = 52) collected during the sampling period. The relationship showed a linear behavior (NO3- = 9.948σ_t − 260.850; R² = 0.930; p < 0.05) for σ_t ranging between 26.1 and 27.1 Kg m⁻³.

During all samplings except February 2014, when only the total fraction was available, size-fractionated (< and > 10 μm) chlorophyll a and primary production were determined at 0, 20, and 40 m. 150-mL samples were filtered onto Whatman GF/F filters for the total fraction, and onto Whatman polycarbonate filters for the >10 μm fraction. Filters for chlorophyll a quantification were extracted in 90% acetone at 4°C overnight and measured using the spectrofluorometric method, as described in Bode et al. (2011). Fluorometrically determined chlorophyll a concentrations, ranging from 0.1 to 5.8 mg m⁻³, were used to calibrate the fluorometer attached to the CTD-rosette (Chl a = 0.533 × fluorescence − 0.212; R² = 0.830, n = 121), which was used to plot the chlorophyll a distribution in February 2014, December 2014, and December 2015. In May 2014, and from April to November 2015, we used the fluorometer included in the MSS profiler (see below), also calibrated with chlorophyll a concentrations ranging from 0.1 to 6.5 mg m⁻³ (Chl a = 3.520 × fluorescence–5.125; R² = 0.770, n = 36), obtained in the same station during the same sampling period. Size-fractionated primary production was measured by ¹⁴C-uptake. Seawater samples were transferred to triplicate 300–mL polycarbonate bottles (Nalgene) (2 light and 1 dark bottles), which were spiked with 2–10 μCi of NaH¹⁴CO₃ and incubated during 24 h in refrigerated incubators equipped with a system of recirculating water, and covered with neutral density screens to simulate the corresponding in situ irradiance. Sample filtration, filter processing and primary production calculations were conducted as described in detail in Bode et al. (2011).

Samples for the determination of picoplankton abundance and cell properties were taken at 0, 40, and 70 m depth. Picoplankton samples (1.8 ml) were preserved with 1% paraformaldehyde + 0.05% glutaraldehyde (final concentration), flash-frozen in liquid N₂ for 10 min and stored at −80°C until analysis with a FACSCalibur flow cytometer (Becton-Dickinson) equipped with a laser emitting at 488 nm. To estimate the abundance of the different groups (cell mL⁻¹), calibration of the cytometer flow rate was performed before each use. A suspension (approximately 1 × 10⁵ mL⁻¹) of fluorescent latex beads 1 μm in diameter (Molecular Probes, Life Technologies) was added to all the samples as internal standard. Two aliquots from the same sample were used for the study of picophytoplankton (0.6 ml) and heterotrophic bacteria (0.4 ml), analyzed at high (~mean 58 μL min⁻¹) and low (~mean 22 μL min⁻¹) flow rate, respectively. Before the analysis, the DNA of heterotrophic bacteria was stained with the nucleic acid fluorochrome SYTO-13 (2.5 μM; Molecular Probes, Life Technologies) in the dark for 10 min.

Autotrophic cells were separated into 2 groups of cyanobacteria (Synechococcus and Prochlorococcus) and 2 groups of picoeukaryotes (small and large), based on their orange (FL2, 585 nm) and red (FL3, >650 nm) fluorescence and side-scattered light (SSC) (Calvo-Díaz and Morán, 2006). Two groups of heterotrophic bacteria, with high and low nucleic acid content (HNA and LNA, respectively), were distinguished based on their relative green fluorescence (FL1, 530 nm), which was used as a proxy for nucleic acid content.

To estimate biovolume, we used an empirical calibration between SSC and cell diameter, assuming spherical shape for all groups. Finally, picoplankton biomass was computed using the following conversion factors of biovolume to carbon: 350 fg C μm⁻³ for heterotrophic bacteria (Bjørnsen, 1986), 230 fg C μm⁻³ for Synechococcus, 240 fg C μm⁻³ for Prochlorococcus and 237 fg C μm⁻³ for picoeukaryotes (Worden et al., 2004).

Water samples from 0, 20, 40, and 70 m depth were collected with Niskin bottles to determine size-fractionated (< and >10 μm) N₂ fixation following the ¹⁵N₂-uptake technique (Montoya et al., 1996), with the modifications described in Rees et al. (2009). In February 2014, only surface water (0 m) was sampled. Triplicate 2-L acid-cleaned clear polycarbonate bottles (Nalgene) were filled directly from the Niskin bottle using acid-washed silicone tubing. After carefully removing all air bubbles, bottles were closed with caps provided with silicone septa, through which 3 mL of ¹⁵N₂ (98 atom%, Isotec Stable Isotopes lot TV0533 in February 2014, and Cambridge Isotope Laboratories lots I-16727 and I-19168A in the remaining samplings) were injected with a gas-tight syringe. The pressure across the septum was equilibrated by allowing the excess water to escape through a syringe needle piercing the septum. Incubation bottles were gently shaken by manually inverting them fifty times, and then incubated for 24 h in the same incubators used for primary production measurements. The incubation of seawater from 70 m (collected in the aphotic zone) was performed under complete darkness, by covering the incubation bottles with black duct tape and placing them inside a black opaque plastic bag.

Incubations ended by filtration through 25 mm Whatman GF/F filters to obtain total N₂ fixation rates. N₂ fixation rates of the <10 μm size-fraction were determined by pre-filtering the water sample through 47 mm Whatman polycarbonate filters of 10 μm pore-size, and subsequently collecting onto GF/F filters. Total and pre-filtered 2-L seawater samples from each depth were also filtered at time zero for the determination of natural abundance of N stable isotopes (δ¹⁵N). After filtration, filters were dried at 40°C during 24 h and stored at room temperature until pelletization in tin capsules. Measurement of particulate organic nitrogen and carbon (PON and POC) content and ¹⁵N atom% were carried out with an elemental analyser combined with a continuous-flow isotope ratio mass-spectrometer (FlashEA112 + Deltaplus, ThermoFinnigan), and using an acetanilide standard as reference. The precision of the analysis, expressed as the standard deviation of the ¹⁵N values determined in a series of 10 standards, was 0.15‰. Isotopic analyses were made at the stable isotope facility (SAI) of the University of A Coruña (Spain). Assuming that the minimum acceptable change of δ¹⁵ N between the initial and the final PON sample is 4‰ (Montoya et al., 1996), and given that the detection limit of the used elemental analyzer is 0.15–0.20 μg N, the detection limit of our method is approximately 0.0005 nmol N L⁻¹ d⁻¹. The equations of Weiss (1970) and Montoya et al. (1996) were used to calculate the initial N₂ concentration (assuming equilibrium with atmosphere) and N₂ fixation rates, respectively. The fact that we obtained in several samplings at different depths BNF rates equal to zero, or below the detection limit, allowed us to discard any effect of potential contamination in the commercial ¹⁵N₂ gas stocks (Dabundo et al., 2014).

Average cell-specific N₂ fixation rates by UCYN-A that would be required to explain the observed BNF rates, were calculated based on their qPCR-estimated abundances and the BNF measured, following the assumptions and calculations indicated in Turk-Kubo et al. (2014). Representative cell-specific rates for UCYN-A ranging from 0.02 to 2.6 fmol N cell⁻¹ h⁻¹ have been directly estimated or modeled by Goebel et al. (2010). Considering these known cell-specific rates, we roughly estimated whether the UCYN-A abundances in our samples could account for the measured bulk BNF rates.

During all samplings except February 2014, when only the sample at 0 m was collected, seawater samples for DNA were collected at 0, 20, 40, and 70 m depth, transferred to acid-cleaned carboys using acid-washed silicone tubes and kept in darkness until further processing in the laboratory. Microbial biomass was collected by filtering 7.5-10 L of seawater through a 0.22 μm Sterivex™ filter unit (Millipore) using a peristaltic pump. The filters were preserved with 1.8 mL of lysis buffer (50 mM Tris-HCl pH 8.3, 40 mM EDTA pH 8.0, 0.75 M sucrose) and stored at −80°C until extraction. DNA was extracted using the PowerWater® DNA Isolation Kit (MO BIO, Laboratories, Inc.), quantified and quality-checked (according to the A₂₆₀/A₂₈₀ ratio) using a spectrophotometer NanoDrop 2000™ (Thermo Fisher Scientific).

The abundance of diazotrophs was determined by quantitative polymerase chain reaction (qPCR) using TaqMan primers-probe sets (PrimeTime® qPCR Assays, Integrated DNA Technologies) targeting Gammaproteobacteria-affiliated phylotype γ-24774A11 (Moisander et al., 2008), UCYN-A1 (Church et al., 2005), and UCYN-A2 (Thompson et al., 2014). However, we are aware that the primers-probe set designed for UCYN-A2 does not contain enough mismatches to avoid cross-hybridization with the UCYN-A3 and UCYN-A4 sublineages (Farnelid et al., 2016). The probes were 5′FAM labeled and double-quenched with ZEN™/3′-IBFQ (Integrated DNA Technologies). The reactions were run on a MyiQ2™ Real-Time PCR Detection System (Bio-Rad Laboratories) at the Instituto Español de Oceanografía (IEO). Final reaction volume of 20 μL contained 10 μL PrimeTime® Gene Expression Master Mix (IDT), 2 μL DNA template, 2 μL of primers-probe set (0.5 and 0.25 μM final concentrations, respectively), and 6 μL PCR grade water (Sigma-Aldrich). Thermal cycling conditions were 95°C for 3 min, followed by 45 cycles of 95°C for 15 s and 60°C for 1 min. Standards, consisting of nine 10-fold serial dilutions containing the targeted nifH fragments (gBlocks® Gene Fragments, IDT), were included with each qPCR run. All samples and standards were run in duplicate. Standard curves were obtained by linear regression of the threshold cycle (Ct) and log₁₀ gene copies per reaction using standards ranging from 10⁸ to 10⁰ gene copies per reaction. Amplification efficiencies were >90% for all reactions. No-template control sample was run in duplicate in all plates and it did not amplify in any runs. We performed inhibition tests by 2- and 5-folds dilutions of all samples and we concluded that our samples were not inhibited. The limit of detection (LOD) and detected but not quantified (DNQ) limits were 4 nifH copies L⁻¹ and 41 nifH copies L⁻¹, respectively. Abundances below LOD were considered 0 copies L⁻¹, whereas nifH copies higher than LOD but less than DNQ were assigned a conservative value of 1 nifH copy per liter.

Measurements of dissipation rates of turbulent kinetic energy were conducted on each cruise during 6–10 vertical casts down to a maximum depth of 70 m. The microstructure turbulence profiler was equipped with 2 velocity microstructure shear sensors (type PNS98), a microstructure temperature sensor (FPO7), a sensor to measure horizontal acceleration of the profiler and a high-precision CTD probe. Fluorometrically determined chlorophyll a concentrations were used to calibrate the fluorometer sensor included in the MSS profiler, as previously indicated. The profiler was balanced to have negative buoyancy and a sinking velocity of ~0.4 to 0.7 m s⁻¹. The frequency of data sampling was 1024 Hz. The sensitivity of the shear sensors was checked after each use. Due to significant turbulence generation close to the ship, only the data below 5 m were considered reliable. Data processing and calculation of dissipation rates of turbulent kinetic energy (ε) was carried out with the commercial software MSSpro as described in detail in Fernández-Castro et al. (2014). The squared Brunt Väisälä frequency (N²), a proxy for water column stratification, was computed from the CTD profiles according to the equation: (1)N2=-(gρw)(∂ρ∂z) (s-2) where g is the acceleration due to gravity (9.8 m s⁻²), ρ_w is seawater density (1,025 kg m⁻³), and ∂_ρ/∂_z is the vertical potential density gradient. After averaging ε and N² over depth intervals of 1 m length, vertical diffusivity or mixing (Kz) was estimated as: (2)Kz=eεN2 (m2s-1) where e is the mixing efficiency, here considered as 0.2 (Osborn, 1980).

Vertical diffusive fluxes of nitrate into the euphotic layer were calculated following the Fick's law, from the product of the nitrate gradient, obtained by linearly fitting nitrate concentrations in the 10–40 m depth layer, and the averaged Kz for the same depth interval. For those samplings in which the microstructure profiler was not deployed (February 2014, December 2014, and December 2015), estimates of Kz under similar hydrographic conditions were considered: October 2014 (for February 2014), and November 2015 (for December 2014 and December 2015). In order to quantify the biogeochemical relevance of BNF, we compared it with the supply of nitrate into the euphotic zone through vertical diffusion.

A simplified calculation of nitrate supply by vertical advection due to upwelling was computed considering the Golfo Ártabro as a single box divided into two layers (Álvarez-Salgado et al., 2000; Villamaña et al., 2017), the deeper one influenced by upwelled water 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³ s⁻¹) would be equal to the incoming bottom flux (Q_B, m³ s⁻¹), computed as the product of the upwelling index (I_W, m³ s⁻¹ km⁻¹) and the length of the mouth of the Golfo Ártabro (ca. 11.5 km). Finally, the nitrate transport into the euphotic layer by vertical advection was computed as: (3)NO3- advective flux = QZAbasin [NO3-]70 m where A_basin is the surface area of the Golfo Ártabro (ca. 400 km²), Q_Z is the vertical advective flux, and [NO3-]_{70 m} is the nitrate concentration determined at 70 m depth for each cruise in the sampling station.

All the statistical analysis were performed by using SPSS v. 22 for Windows (IBM SPSS Statistics).

During our survey we encountered the main hydrographic features of a temperate coastal region subjected to seasonal variability and influenced by upwelling pulses, which in this system induce large variability at short temporal scales (Bode et al., 1996; Casas et al., 1997). During downwelling conditions, we sampled winter (February 2014 and December 2014) and autumn (November 2015 and December 2015) mixing, and also a spring transitional period under prevailing predominance of downwelling (May 2015) (Figure 2). Autumn mixing was characterized by thermohaline mixing, whereas during winter mixing the slight haline stratification caused mixed layer depth being shallower (Figure 3, Table S1). In May 2015, the higher vertical gradient in temperature and salinity at the surface produced that the mixed layer depth was even shallower. However, during downwelling averaged vertical stratification (determined as Brunt Väisälä frequency) between 10 and 40 m was relatively low (4.4 ± 7.1 × 10⁻⁵ s⁻²), and as a result the mixed layer depth was significantly deeper (50 ± 22 m), than during upwelling and relaxation conditions (Table 1, Figure S1). Microstructure turbulence observations, which were only available for May and November 2015, showed relatively low values of averaged vertical (10–40 m) dissipation rates of turbulent kinetic energy (ε) in both cruises (1.3–2.2 × 10⁻⁸ m² s⁻³).

Cruises which sampled spring-summer upwelling (May 2014, April 2015, and June 2015) were characterized by shallow mixed layer depths (9 ± 3 m) (Table 1, Figure S1). Due to the enhanced values observed in the upper 20 m in May 2014, averaged dissipation rates were slightly higher during this cruise (4.3 ± 7.3 × 10⁻⁸ m² s⁻³) (Figure 3, Table S1).

Relaxation conditions included two cruises (July 2015 and September 2015), which sampled the transition between intense upwelling and downwelling events during the summer thermal stratification period (Figure 2). Both cruises were characterized by intense vertical stratification (19.7 ± 21.9 × 10⁻⁵ s⁻²), and shallow mixed layer depths (12 ± 1 m) (Table 1, Figure 3). As a result of the relatively low vertical stratification and enhanced dissipation rates observed between 25 and 40 m depth layer, averaged vertical diffusivity (12.0 ± 33.0 × 10⁻⁴ m² s⁻¹) was higher, although not statistically significant, than during downwelling (4.5 ± 7.6 × 10⁻⁴ m² s⁻¹) and upwelling conditions (1.5 ± 4.2 × 10⁻⁴ m² s⁻¹).

In downwelling conditions, due to the low stratification and deeper mixed layers, nitrate and phosphate vertical profiles were relatively homogeneous, whereas during upwelling, the input of cold and high-nitrate waters at depth caused significant enhanced nutrient vertical gradients (Figure 4). Vertical gradients of nitrate were also higher during relaxation, and the concentrations at depth of both nitrate and phosphate were higher during September 2015, compared to July 2015, due to the intensifying upwelling observed one week before the September 2015 sampling (Figure 2). On average, the vertical nitrate gradient between 10 and 40 m depth during upwelling (148 ± 41 μmol m⁻⁴) and relaxation (139 ± 2 μmol m⁻⁴) was significantly higher than during downwelling (15 ± 13 μmol m⁻⁴) (Table 1).

Vertical profiles of total chlorophyll a and primary production were also homogeneous during downwelling (Figure 5). During these conditions depth-integrated chlorophyll a and primary production were, on average, relatively low (31 ± 16 mg m⁻² and 1,359 ± 1,032 mg C m⁻² d⁻¹, respectively), and dominated (>80 %) by small (<10 μm) phytoplankton cells (Table 1, Figure S1). During upwelling the fertilization effect of deep waters stimulates phytoplankton growth, and increases the contribution of larger cells (>10 μm). Averaged depth-integrated chlorophyll a (65 ± 35 mg m⁻²) and primary production (5,192 ± 538 mg C m⁻² d⁻¹) were slightly and significantly, respectively, higher than during downwelling and relaxation, and the contribution of larger cells to both chlorophyll a and primary production was higher than 60%. During relaxation depth-integrated chlorophyll a (20 ± 2 mg m⁻²) and primary production (1,611 ± 295 mg C m⁻² d⁻¹) were relatively low. In September 2015, small phytoplankton cells contributed more than 80% to both chlorophyll a and primary production, whereas in July 2015 the contribution was lower than 50%.

Flow cytometry data showed that picoplankton biomass was largely dominated by heterotrophic bacteria (HNA and LNA), followed by picoeukaryotes (large and small), and finally cyanobacteria (Synechococcus and Prochlorococcus) (Table 1, Table S1, Figure S2). The averaged contribution of HNA bacteria to total biomass peaked (ca. 60%) during upwelling and relaxation conditions. Consistently with the information provided by size-fractionated chlorophyll a, which showed the predominance of small (<10 μm) cells during downwelling (Figure 5), the averaged depth-integrated biomass of total picophytoplankton (i.e., large and small picoeukaryotes plus Synechococcus and Prochlorococcus cyanobacteria) (88 ± 48 mg C m⁻²), was during this condition slightly higher than during upwelling (54 ± 23 mg C m⁻²), and relaxation (31 ± 26 mg C m⁻²).

The vertical distribution of picoeukaryotic biomass, which potentially would include the UCYN-A diazotroph-prymnesiophyte host (Thompson et al., 2012, 2014), was relatively homogeneous during downwelling, with the exception of May 2015, when the highest biomass (>3 μg C L⁻¹) was observed in the upper 40 m (Figure 6). During upwelling conditions, the vertical distribution of picoeukaryotes biomass showed small vertical variability in May 2014, whereas in April and June 2015 enhanced values (>2 μg C L⁻¹) were observed at the surface. During relaxation, small vertical variability was observed in July 2015, whereas in September 2015 enhanced values (>1.5 μg C L⁻¹) also occurred at surface waters.

No significant differences (Kruskal-Wallis, p > 0.05) were found between the BNF rates corresponding to the total and small (<10 μm) size-fractions (data not shown), suggesting that small diazotrophs, presumably unicellular, were responsible for all the BNF activity detected in the region. For this reason, only BNF rates corresponding to the total size-fraction are described here. Whereas no significant differences were observed in depth-integrated rates, which ranged from 0.1 ± 0.1 to 1.6 ± 0.5 μmol N m⁻² d⁻¹ (Table 1), contrasting vertical patterns were observed during downwelling, upwelling and relaxation conditions (Figure 6).

During downwelling, the vertical distribution of BNF was relatively homogeneous, and no significant differences were found between the rates determined at different depths (Kruskal-Wallis, p > 0.05) (Figure 6). Depth-integrated rates were slightly higher in December 2014 (1.3 ± 0.4 μmol N m⁻² d⁻¹) and May 2015 (1.0 ± 0.3 μmol N m⁻² d⁻¹), compared to November 2015 (0.3 ± 0.1 μmol N m⁻² d⁻¹) and December 2015 (0.1 ± 0.1 μmol N m⁻² d⁻¹). During both upwelling and relaxation conditions, maximum rates were observed in general at the surface, decreasing between 20 and 40 m and slightly increasing again at 70 m. Depth-integrated rates were higher in April 2015 (1.6 ± 0.5 μmol N m⁻² d⁻¹, upwelling) and September 2015 (0.8 ± 0.1 μmol N m⁻² d⁻¹, relaxation) compared to June 2015 (0.6 ± 0.3 μmol N m⁻² d⁻¹, upwelling). As a result of the hydrographic patterns described, when considering all the samplings together, BNF was negatively correlated with mixed layer depth (Pearson's r = −0.364, p < 0.05), whereas a positive significant relationship was found with the primary production corresponding to the small (<10 μm) size-fraction (Pearson's r = 0.573, p < 0.01).

Due to the increased vertical diffusivity and vertical nitrate gradient determined during relaxation conditions, averaged nitrate diffusive fluxes during these samplings were more than six-fold higher (14.3 ± 3.1 mmol m⁻² d⁻¹), compared to upwelling (2.2 ± 3.1 mmol m⁻² d⁻¹), and downwelling (0.8 ± 0.9 mmol m⁻² d⁻¹) conditions (Table 1, Figures 3, 4). These numbers were between 2 and 5 orders of magnitude higher than the estimates of depth-integrated BNF (0.1 ± 0.1–1.6 ± 0.5 μmol m⁻² d⁻¹) carried out during the NICANOR samplings. From a biogeochemical perspective, the comparison between these two processes evidences that the relevance of BNF in this system, which was relatively higher during the upwelling conditions sampled in April 2015, was always <2% (Table S1).

Clone library analysis of the nitrogenase gene (nifH) recovered from DNA samples obtained during the NICANOR samplings revealed that most of sequenced clones belonged to unicellular cyanobacterial and heterotrophic bacterial diazotrophs (Moreira-Coello et al., in preparation). nifH clusters 1G (Gammaproteobacteria), 1B (Cyanobacteria) and 3E (Verrucomicrobia and Deltaproteobacteria) were the most abundant, representing 29, 21, and 22% of total nifH clones, respectively. The other clusters represented less than 8% of relative abundances each. Most of clones belonged to Pseudomonas stutzeri (1G) (21%), UCYN-A (1B) (21%), Opitutaceae bacterium (3E) (9%) and Desulfovibrio aespoeensis (3E) (8%), whereas the rest of diazotrophs represented <7% each. Considering that the two most abundant groups were UCYN-A and Gammaproteobacteria-affiliated phylotypes (1G), we used the qPCR primers-probe sets available for quantifying UCYN-A sublineages (-A1 and -A2) and Gammaproteobacteria γ-24774A11 phylotype, widely distributed and considered one of the most important heterotrophic diazotrophs (e.g., Moisander et al., 2014; Langlois et al., 2015). The quantification of these diazotrophs allowed examining their relative importance through the year. All the abundances were above the limit of detection, except the abundance of UCYN-A1 determined at the surface in February 2014 (downwelling), and the abundances of Gammaproteobacteria determined at the surface and 40 m depth in May 2014 and April 2015 (upwelling) (Figure 6).

The Gammaproteobacteria phylotype was significantly more abundant than UCYN-A1 and UCYN-A2 during downwelling and relaxation (Bonferroni comparisons, p < 0.05) (Table 1, Figure 6). Although not statistically significant, both surface (3.1 ± 5.3 × 10² nifH copies L⁻¹) and depth-integrated (13 ± 23 × 10⁵ nifH copies m⁻²) abundances of Gammaproteobacteria were slightly lower during upwelling, as in May 2014 and April 2015 the abundances between 0 and 40 m depth were below LOD or DNQ limits (Table S1, Figure S1). During downwelling and relaxation conditions, Gammaproteobacteria showed vertical variability, reaching maximum abundances at 70 m depth. In November 2015 (downwelling) and September 2015 (relaxation) we found the highest surface (9.0 × 10³ and 1.9 × 10³ nifH copies L⁻¹, respectively) and depth-integrated abundances (380 × 10⁵ and 220 × 10⁵ nifH copies m⁻², respectively) measured on all cruises.

During relaxation, UCYN-A1 abundance increased at the surface, and both surface (6.2 ± 0.7 × 10² nifH copies L⁻¹) and depth-integrated (13 ± 0 × 10⁵ nifH copies m⁻²) abundance were significantly higher than during downwelling (33 ± 64 nifH copies L⁻¹ and 1.4 ± 1.3 × 10⁵ nifH copies m⁻², respectively) and upwelling conditions (1.9 ± 0.3 × 10² nifH copies L⁻¹ and 1.3 ± 0.2 × 10⁵ nifH copies m⁻², respectively), when the vertical distribution of this group did not exhibit important vertical variability (Table 1, Figure 6).

The abundance of UCYN-A2 was higher than UCYN-A1 (Kruskal-Wallis, p < 0.05) and showed larger vertical variability (Figure 6, Table S1). On average, surface and depth-integrated abundances of UCYN-A2 were slightly higher during relaxation (8.7 ± 8.6 × 10² nifH copies L⁻¹ and 21 ± 16 × 10⁵ nifH copies m⁻², respectively), compared to upwelling (4.3 ± 5.6 × 10² nifH copies L⁻¹ and 10 ± 10 × 10⁵ nifH copies m⁻², respectively) and downwelling (1.4 ± 2.4 × 10² nifH copies L⁻¹ and 7 ± 7 × 10⁵ nifH copies m⁻², respectively) (Table 1). During downwelling, higher values occurred at 70 m in November 2015 (6.9 × 10² nifH copies L⁻¹), and also at the surface in December 2015 (5.0 × 10² nifH copies L⁻¹). During upwelling and relaxation, the maximum abundances were quantified at the surface in June 2015 (1.1 × 10³ nifH copies L⁻¹) and September 2015 (1.5 × 10³ nifH copies L⁻¹), respectively.

Our results demonstrate the existence of detectable BNF activity (0.001 ± 0.002–0.095 ± 0.024 μmol N m⁻³ d⁻¹) under contrasting hydrographic regimes in a nitrogen-rich temperate coastal system subjected to upwelling pulses. During downwelling-influenced cruises the vertical pattern of BNF was rather homogeneous, whereas during upwelling and relaxation the BNF rates peaked in the well-lit surface waters, where we measured the highest rates of all cruises. Comparing with other studies using the same technique (¹⁵N₂-gas tracer method, Montoya et al., 1996) and performed in similar nitrate-rich and temperate areas, our volumetric BNF rates were of similar magnitude to those measured by Benavides et al. (2011) in Cape Silleiro (NW Iberian Peninsula) in summer 2009, but between 1 and 3 orders of magnitude lower than those measured by Mulholland et al. (2012) in the temperate NE American coast during several summer and autumn cruises, and by Rees et al. (2009) in the western English Channel in summer. On a global context, our BNF rates estimated in the outer part of Ría de A Coruña are similar to the lower-end described for the subtropical NE Atlantic (Luo et al., 2012). Although our rates were relatively low, our results contradict the traditional paradigm about marine BNF. This view established that BNF is mainly restricted to the nitrogen-poor, warm (>20°C) (sub)tropical ocean, as high nutrient inputs and the relatively low surface temperatures in temperate coastal systems would not favor N₂ fixation (Howarth et al., 1988; Conley et al., 2009; Stal, 2009). The BNF activity detected in this study is a further evidence that high nitrate concentrations (and high N:P ratios) do not inhibit BNF activity, such as it was also observed in other upwelling regions (e.g., Voss et al., 2004; Sohm et al., 2011; Subramaniam et al., 2013; Benavides et al., 2014). Indeed, recent studies demonstrated that BNF can occur upon exposure to elevated DIN (nitrate and/or ammonium), as long as phosphate is available (Knapp, 2012), and point toward iron as the primary control factor of BNF (Knapp et al., 2016; Bonnet et al., 2017). Thus, the sensitivity of marine BNF to high DIN is not clear. Along with the wide distribution and large diversity of diazotrophs (Zehr et al., 2000; Moisander et al., 2010), our results indicate that restricting BNF measurements exclusively to oligotrophic (sub)tropical regions may result in an underestimation of the global magnitude of this flux, contributing to the apparent imbalances observed in the marine N budget (Brandes and Devol, 2002; Mahaffey et al., 2005; Codispoti, 2007).

The contribution of BNF to new N supply into the euphotic zone was negligible (<2%) in this system, as it was between 2 and 5 orders of magnitude lower than diffusive fluxes. These numbers represent the first estimate of the contribution of both processes in a coastal upwelling system. Few studies have simultaneously quantified BNF and nitrate diffusion (Mouriño-Carballido et al., 2011; Painter et al., 2013; Fernández-Castro et al., 2015), mainly due to the methodological difficulties to estimate vertical mixing (Kz) in the field, which has frequently motivated the use of constant values of Kz (Capone et al., 2005), and empirical parameterizations (Planas et al., 1999; Fernández-Castro et al., 2014). Capone et al. (2005) estimated in the tropical North Atlantic that BNF could equal or even exceed nitrate vertical diffusion into the euphotic zone. In the subtropical NE Atlantic, Mouriño-Carballido et al. (2011) reported a contribution of BNF to new nitrogen input of 2 ± 2%, whereas Painter et al. (2013) indicated that BNF represented a significant source of new N, but not as large as the diffusive fluxes. In their study across the (sub)tropical Atlantic, Pacific and Indian oceans, Fernández-Castro et al. (2015) pointed out that, on average and compared to BNF, nitrate diffusion was the main mechanism supplying new N. The magnitude of averaged nitrate diffusion estimated in our study (3.9 ± 5.8 mmol m⁻² d⁻¹) was slightly higher than the estimations for the subtropical NE Atlantic by Mouriño-Carballido et al. (2011) (0.8 ± 1.3 mmol m⁻² d⁻¹) and Painter et al. (2013) (0.06 ± 0.02 mmol m⁻² d⁻¹), and by Fernández-Castro et al. (2015) (0.2 ± 0.2 mmol m⁻² d⁻¹) in the (sub)tropical Atlantic, Pacific and Indian oceans. The range of diffusive fluxes here reported (0.1–16.5 mmol m⁻² d⁻¹) was comparable to that obtained in a nearby coastal system (Ría de Vigo, NW Iberia) by Cermeño et al. (2016) (0.8–20.1 mmol m⁻² d⁻¹) and Villamaña et al. (2017) (9.4–34.6 mmol m⁻² d⁻¹).

We cannot discard that our N₂ fixation rates may be underestimated (Mohr et al., 2010; Großkopf et al., 2012) due to the use of the ¹⁵N₂-bubble addition technique (Montoya et al., 1996) instead of the modified tracer method with ¹⁵N₂-enriched seawater (Mohr et al., 2010). However, the biogeochemical relevance of BNF in this system would be still very low since other important mechanisms of new N supply, such as advective vertical and horizontal transport, atmospheric deposition or internal waves (Oschlies and Garçon, 1998; Duce et al., 2008; Villamaña et al., 2017), were not taken into account in this study.

In order to get a rough estimation of the fertilization effect of the upwelling, we compared a simplified estimate of nitrate supply through vertical advection (see Materials and Methods) with the input of nitrogen through nitrate vertical diffusion. As expected, nitrate vertical advection was the main mechanism of N supply during upwelling (79 ± 119 mmol m⁻² d⁻¹) and during relaxation in September 2015 (77 mmol m⁻² d⁻¹), which was probably influenced by an upwelling pulse happening a few days before the sampling. This nitrogen supply was more than 180-fold higher than nitrate diffusion during upwelling, and up to 6-fold higher in September 2015. Obviously, non-upwelling situations implied higher contribution to new N supply by diffusion. In any case, the consideration of this process would decrease, even more, the contribution of BNF to the new nitrogen supply in this system. Overall, the proportion of gross primary production (GPP) that could be sustained by diazotrophy, calculated by assuming Redfield stoichiometry, would be irrelevant (<0.1%) when compared to the proportion of GPP supported by nitrate (41.9 ± 31.2%) and ammonium (6.8 ± 5.0%) uptake rates, measured in the same region in an earlier study (Bode et al., 2004).

Small diazotrophs (<10 μm), presumably unicellular, were responsible for all BNF activity detected in this region. This is consistent with analysis of the nifH clone library, which revealed that most of sequenced clones belonged to unicellular cyanobacterial and heterotrophic diazotrophs (Moreira-Coello et al., in preparation). Most of nifH clones belonged to Gammaproteobacteria-affiliated phylotypes of the cluster 1G (29%) and UCYN-A (21%).

The abundance of Gammaproteobacteria γ-24774A11 phylotype (0–2.4 × 10⁴ nifH copies L⁻¹) was significantly higher than UCYN-A1 and UCYN-A2 during downwelling and relaxation (Bonferroni comparisons, p < 0.05). These values were in the medium-high range of global abundances reviewed by Luo et al. (2012), and by Benavides and Voss (2015) specifically from studies performed in the North Atlantic. In fact, heterotrophic diazotrophs abundances higher than 10⁴ nifH copies L⁻¹ have seldom been reported in the marine environment (e.g., Moisander et al., 2008; Halm et al., 2012). The abundance of the heterotrophic Gammaproteobacteria was positively correlated with temperature (Pearson's r = 0.473, p < 0.05) and mixed layer depth (Pearson's r = 0.495, p < 0.05), since slightly higher abundances were measured during downwelling conditions, when the water column was warmer, less stratified and probably depleted in bioessential nutrients. These hydrographic features could cause the dominance of Gammaproteobacteria as the primary productivity was low and dominated by small cells, and likely sustained by regenerated nutrients (mainly nitrogen) through bacterial remineralisation of organic matter, such as it was observed previously in this region (Bode et al., 2004, 2011). The production based on bacterial recycling of nutrients may also have occurred in September 2015, when the abundance of Gammaproteobacteria increased significantly, probably fueled by the organic matter produced during the upwelling pulse happening before the sampling (Bode et al., 1996; Casas et al., 1997).

The symbiosis of UCYN-A with picoeukaryotic prymnesiophyte cells also thrives in these nitrogen-rich and relatively cold (<17°C) waters. Similar findings were also reported by Agawin et al. (2014) and Benavides et al. (2011) in the same upwelling region, and by Mulholland et al. (2012) in the N-rich temperate NE American coast. Furthermore, UCYN-A1 and its small uncultured host were detected in the Galician Bank, an open ocean environment relatively close to our region (Cabello et al., 2015). Foster et al. (2009) also found abundant UCYN-A in the nutrient-rich Equatorial Atlantic upwelling area in summer 2007.

The pattern of vertical distribution of picoeukaryotes was consistent with those of UCYN-A. When the abundances of UCYN-A peaked, the biomass of the picophytoplankton also peaked (but not vice versa), except in July 2015, when the consistency was observed in the profile of abundance (data not shown). Indeed, qPCR analysis of picoeukaryotic populations sorted with flow cytometry unequivocally demonstrated that UCYN-A are associated with picoeukaryotes (Thompson et al., 2012, 2014). Furthermore, we did not find a link between UCYN-A abundance and phytoplankton biomass (Chlorophyll a), which could be expected on the basis of the symbiotic association of UCYN-A with a photosynthetic prymnesiophyte in regions where picoplankton dominates phytoplankton biomass, such as it was previously observed in the oligotrophic (sub)tropical SW Pacific (Moisander et al., 2010). This decoupling was probably due to the fact that our study was carried out in an upwelling system, where usually most biomass corresponds to large phytoplankton (>10 μm), mainly chain-forming diatoms (Casas et al., 1997; Varela et al., 2001).

The abundance of UCYN-A2 sublineage (1–1.5 × 10³ nifH copies L⁻¹) was on average 3-fold higher than that of UCYN-A1 (0–6.7 × 10² nifH copies L⁻¹). The higher UCYN-A2 abundance is consistent with recent findings pointing out this sublineage, although globally distributed (Cabello et al., 2015; Martínez-Pérez et al., 2016), mainly inhabit coastal regions, whereas UCYN-A1 dominates in the oligotrophic open ocean (Thompson et al., 2014; Messer et al., 2015; Turk-Kubo et al., 2017). The UCYN-A abundances detected were in the medium-low range of global abundances reported by Luo et al. (2012), Martínez-Pérez et al. (2016), and Benavides and Voss (2015) in the North Atlantic. In the temperate NE American coast, Mulholland et al. (2012) showed that UCYN-A was the most abundant diazotroph, reaching 3.5 × 10⁷ nifH copies L⁻¹, among the highest abundances ever observed. Several studies performed in the North Atlantic reported that, after Trichodesmium, UCYN-A is the most abundant diazotrophic phylotype (Benavides and Voss, 2015; Ratten et al., 2015; Benavides et al., 2016). Even though it extends along a wide latitudinal and longitudinal range, it dominates mainly the eastern basin, where iron inputs by Saharan dust are higher (Benavides et al., 2013), which favors its growth and N₂-fixing activity (Krupke et al., 2015).

Higher abundances of UCYN-A1 and UCYN-A2 occurred at the well-lit surface waters during upwelling and relaxation conditions, suggesting a possible seasonal variability on the growth of the UCYN-A–prymnesiophyte symbiosis, such as it was previously observed in the (sub)tropical SW Pacific (Moisander et al., 2010; Bonnet et al., 2015). These high values coincided with the highest BNF rates, and in fact UCYN-A abundances were positively correlated with BNF rates (Pearson's r = 0.534, p < 0.05), consistent with the fact that small diazotrophs (<10 μm) were responsible for all BNF activity in this system. This relationship has also been observed in other regions, as the NE American coast (Mulholland et al., 2012), or the tropical SW Pacific (Bonnet et al., 2015). Therefore, UCYN-A–prymnesiophyte associations may be responsible for most BNF activity, or at least a large fraction, in our region, since Gammaproteobacteria abundance and BNF rates were not correlated or consistent. Turk-Kubo et al. (2014) reported that Gammaproteobacteria phylotypes do not contribute significantly to marine BNF, based on their analysis of heterotrophic cell-specific N₂ fixation rates required to explain previously reported BNF rates in different marine environments (e.g., Zehr et al., 2007; Großkopf et al., 2012; Halm et al., 2012; Farnelid et al., 2013). In our study, average cell-specific N₂ fixation rates required to account for BNF rates (see Material and Methods) ranged from 0.8 fmol N cell⁻¹ h⁻¹ (November 2015, 40 m) to 19.3 fmol N cell⁻¹ h⁻¹ (April 2015, 0 m). Therefore, in several samples UCYN-A may have been responsible for all the BNF activity, since cell-specific rates estimated for this diazotrophic group range from 0.02 to 2.6 fmol N cell⁻¹ h⁻¹ (Goebel et al., 2010). However, in some samples the required cell-specific rates are too high to be due solely to the UCYN-A quantified. Thus, BNF rates may result from small contributions from Gammaproteobacteria phylotypes or other diazotrophs that we have not quantified by qPCR.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.