Biological N2 Fixation in the Upwelling Region off NW Iberia: Magnitude, Relevance, and Players

The classical paradigm about marine N2 fixation establishes that this process is mainly constrained to nitrogen-poor tropical and subtropical regions, and sustained by the colonial cyanobacterium Trichodesmium spp. and diatom-diazotroph symbiosis. However, the application of molecular techniques allowed determining a high phylogenic diversity and a wide distribution of marine diazotrophs, which extends the range of ocean environments where biological N2 fixation may be relevant. Between February 2014 and December 2015, we carried out 10 one-day samplings in the upwelling system off NW Iberia in order to: 1) investigate the seasonal variability in the magnitude of N2 fixation, 2) determine its biogeochemical role as a mechanism of new nitrogen supply, and 3) quantify the main diazotrophs in the region under contrasting hydrographic regimes. Our results indicate that the magnitude of N2 fixation in this region was relatively low (0.001±0.002 – 0.095±0.024 µmol N m-3 d-1), comparable to the lower-end of rates described for the subtropical NE Atlantic. Maximum rates were observed at surface during both upwelling and relaxation conditions. The comparison with nitrate diffusive fluxes revealed the minor role of N2 fixation (2 fixation activity detected in the region. Quantitative PCR targeting the nifH gene revealed the highest abundances of two sublineages of Candidatus Atelocyanobacterium thalassa or UCYN-A (UCYN-A1 and UCYN-A2) mainly at surface waters during upwelling and relaxation conditions, and of Gammaproteobacteria γ-24774A11 at deep waters during downwelling. Maximum abundance for the three groups were up to 6.7 × 102, 1.5 × 103 and 2.4 × 104 nifH copies L-1, respectively. Our findings demonstrate measurable N2 fixation activity and presence of diazotrophs throughout the year in a nitrogen-rich temperate region.

The classical paradigm about marine N 2 fixation establishes that this process is mainly constrained to nitrogen-poor tropical and subtropical regions, and sustained by the colonial cyanobacterium Trichodesmium spp. and diatom-diazotroph symbiosis. However, the application of molecular techniques allowed determining a high phylogenic diversity and wide distribution of marine diazotrophs, which extends the range of ocean environments where biological N 2 fixation may be relevant. Between February 2014 and December 2015, we carried out 10 one-day samplings in the upwelling system off NW Iberia in order to: (1) investigate the seasonal variability in the magnitude of N 2 fixation, (2) determine its biogeochemical role as a mechanism of new nitrogen supply, and (3) quantify the main diazotrophs in the region under contrasting hydrographic regimes. Our results indicate that the magnitude of N 2 fixation in this region was relatively low (0.001 ± 0.002 -0.095 ± 0.024 µmol N m −3 d −1 ), comparable to the lower-end of rates described for the subtropical NE Atlantic. Maximum rates were observed at the surface during both upwelling and relaxation conditions. The comparison with nitrate diffusive fluxes revealed the minor role of N 2 fixation (<2%) as a mechanism of new nitrogen supply into the euphotic layer. Small diazotrophs (<10 µm) were responsible for all N 2 fixation activity detected in the region. Quantitative PCR targeting the nifH gene revealed the highest abundances of two sublineages of Candidatus Atelocyanobacterium thalassa or UCYN-A (UCYN-A1 and UCYN-A2), mainly at surface waters during upwelling and relaxation conditions, and of Gammaproteobacteria γ-24774A11 at deep waters during downwelling. Maximum abundance for the three groups were up to 6.7 × 10 2 , 1.5 × 10 3 , and 2.4 × 10 4 nifH copies L −1 , respectively. Our findings demonstrate measurable N 2 fixation activity and presence of diazotrophs throughout the year in a nitrogen-rich temperate region.

INTRODUCTION
Dinitrogen (N 2 ) is the most abundant form of nitrogen (N) in aquatic and terrestrial ecosystems, however only a limited, but diverse, group of bacteria and archaea (termed diazotrophs) can use this large N reservoir by the energy-costly process of biological N 2 fixation (BNF) (Karl et al., 2002). Diazotrophs possess nitrogenase, the enzyme that catalyzes this process, which reduces atmospheric N 2 to bioavailable ammonium (NH + 4 ). With an estimated global flux of ∼140-177 Tg N year −1 , BNF is the main mechanism that supplies new N to the ocean, fueling phytoplanktonic photosynthesis and the subsequent export of organic matter into the deep waters through the "biological carbon pump" (Voss et al., 2013). The limited number of studies simultaneously quantifying BNF vs. the transport of nitrate into the euphotic zone through turbulent mixing indicate that BNF can equal or even exceed nitrate diffusion in some regions of the subtropical gyres Mouriño-Carballido et al., 2011;Painter et al., 2013;Fernández-Castro et al., 2015). Therefore, BNF is a key process controlling oceanic productivity, the carbon cycle, and climate (Gruber and Galloway, 2008).
Traditionally, it has been considered that marine BNF is mainly restricted to nitrogen-poor tropical and subtropical regions, and sustained by the filamentous cyanobacterium Trichodesmium spp. and diatom-diazotroph associations. However, the development and application of molecular techniques allowed to determine an extremely high phylogenic diversity of nitrogenase genes (nifH) (Zehr et al., 2000). In addition to Trichodesmium and diatom-diazotroph associations, marine diazotrophic groups include unicellular cyanobacteria (UCYN-A, -B, -C) and non-cyanobacterial diazotrophs (heterotrophic bacteria and archaea) (Zehr et al., 2003). Diazotrophic heterotrophic bacteria comprise a wide range of phylogenetic groups, including Firmicutes, Alpha-, Beta-, Gamma-, Delta-and Epsilonproteobacteria (Zehr et al., 2003).
Although numerous studies have investigated the biogeography and controlling mechanisms of marine BNF in recent years (e.g., Moore et al., 2009;Monteiro et al., 2011;Luo et al., 2014;Fernández-Castro et al., 2016), the limiting factors of this process are still under debate. Due to the high energetic cost of fixing N 2 in comparison to the assimilation of dissolved inorganic nitrogen (DIN) (mainly ammonia and nitrate), DIN-replete environments have been considered as inhibitory for the growth of diazotrophs and BNF activity (Holl and Montoya, 2005). However, the discovery of a large diversity of organisms able to fix N 2 , and their wide distribution, indicated that BNF may be more widespread than previously thought (Moisander et al., 2010). Recent studies demonstrate the presence and activity of diazotrophs in DIN-enriched environments such as the eastern tropical North Atlantic (Voss et al., 2004), the western English Channel (Rees et al., 2009), the Mekong River plume in the South China Sea (Grosse et al., 2010), the temperate NE American coast (Mulholland et al., 2012), and mesohaline temperate estuaries between the Baltic Sea and the North Sea (Bentzon-Tilia et al., 2015). Evidences of BNF have also been reported in upwelling systems such as the eastern equatorial Atlantic (Voss et al., 2004;Foster et al., 2009;Subramaniam et al., 2013), the Benguela upwelling system (Sohm et al., 2011;Benavides et al., 2014), and the region off central Chile (Raimbault and Garcia, 2008;. The NW Iberian Peninsula is part of the Iberia-Canary upwelling domain, one of the eastern boundary upwelling ecosystems in the global ocean (Arístegui et al., 2009). This area is characterized by the presence of coastal embayments named Rías and the seasonal influence of wind-driven upwelling and downwelling events (Wooster et al., 1976;Fraga, 1981). From April to September, the prevailing northerly winds cause upwelling, and the subsequent rise of cold and nutrient-rich subsurface waters into the coastal euphotic layer. These nutrient inputs, and their amplification due to remineralisation inside the Rías (Bode et al., 1996;Álvarez-Salgado et al., 2002), fuel the high phytoplankton production that supports a rich marine ecosystem and productive fisheries (Fraga, 1981;Varela, 1992). From October to March, the prevailing southerly winds favor downwelling events (Wooster et al., 1976;Fraga, 1981). Knowledge about diazotrophic activity in the NW Iberia region is limited to a single observation carried out in summer 2009, when relatively low BNF rates, mainly attributed to UCYN-A, were reported by Agawin et al. (2014) and Benavides et al. (2011). Here we describe the results of a survey conducted in the outer part of Ría de A Coruña between February 2014 and December 2015, which covered contrasting hydrographic regimes. Our main goals are to: (1) investigate the seasonal variability in the magnitude of BNF, (2) determine its biogeochemical role as a mechanism of new N supply, and (3) quantify the main diazotrophs.

Sampling and Hydrographic Measurements
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, andSeptember 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 3 s −1 km −1 ) 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 2 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).

Inorganic Nutrients, Chlorophyll a, and Primary Production
Unfiltered samples for the determination of dissolved inorganic nutrients (NO − 3 , PO 3− 4 ) 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 (NO − 3 = 9.948σ t − 260.850; R 2 = 0.930; p < 0.05) for σ t ranging between 26.1 and 27.1 Kg m −3 .
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 −3 , were used to calibrate the fluorometer attached to the CTD-rosette (Chl a = 0.533 × fluorescence − 0.212; R 2 = 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 −3 (Chl a = 3.520 × fluorescence-5.125; R 2 = 0.770, n = 36), obtained in the same station during the same sampling period. Size-fractionated primary production was measured by 14 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 14 CO 3 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).

Flow Cytometric Counting
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 2 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 −1 ), calibration of the cytometer flow rate was performed before each use. A suspension (approximately 1 × 10 5 mL −1 ) 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 −1 ) and low (∼mean 22 µL min −1 ) 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.

Biological N 2 Fixation
Water samples from 0, 20, 40, and 70 m depth were collected with Niskin bottles to determine size-fractionated (< and >10 µm) N 2 fixation following the 15 N 2 -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 15 N 2 (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 gastight 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 2 fixation rates. N 2 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 (δ 15 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 15 N atom% were carried out with an elemental analyser combined with a continuous-flow isotope ratio massspectrometer (FlashEA112 + Deltaplus, ThermoFinnigan), and using an acetanilide standard as reference. The precision of the analysis, expressed as the standard deviation of the 15 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 δ 15 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 −1 d −1 . The equations of Weiss (1970) and Montoya et al. (1996) were used to calculate the initial N 2 concentration (assuming equilibrium with atmosphere) and N 2 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 15 N 2 gas stocks (Dabundo et al., 2014).
Average cell-specific N 2 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 −1 h −1 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.

DNA Samples Collection and Extraction
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 TM 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 R DNA Isolation Kit (MO BIO, Laboratories, Inc.), quantified and quality-checked (according to the A 260 /A 280 ratio) using a spectrophotometer NanoDrop 2000 TM (Thermo Fisher Scientific).

Quantification of nifH Gene by qPCR
The abundance of diazotrophs was determined by quantitative polymerase chain reaction (qPCR) using TaqMan primers-probe sets (PrimeTime R 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 TM /3 ′ -IBFQ (Integrated DNA Technologies). The reactions were run on a MyiQ2 TM 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 R 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 R 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 10 gene copies per reaction using standards ranging from 10 8 to 10 0 gene copies per reaction. Amplification efficiencies were >90% for all reactions. Notemplate control sample was run in duplicate in all plates and it did not amplify in any runs. We performed inhibition tests by 2and 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 −1 and 41 nifH copies L −1 , respectively. Abundances below LOD were considered 0 copies L −1 , whereas nifH copies higher than LOD but less than DNQ were assigned a conservative value of 1 nifH copy per liter.

Dissipation Rates of Turbulent Kinetic Energy and Estimates of Nitrate Fluxes through Vertical Diffusion
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 −1 . 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 2 ), a proxy for water column stratification, 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 seawater density (1,025 kg m −3 ), and ∂ ρ /∂ z is the vertical potential density gradient. After averaging ε and N 2 over depth intervals of 1 m length, vertical diffusivity or mixing (Kz) was estimated as: 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 ( In order to quantify the biogeochemical relevance of BNF, we compared it with the supply of nitrate into the euphotic zone through vertical diffusion.

Nitrate Supply through Vertical Advection
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 3 s −1 ) would be equal 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 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: where A basin is the surface area of the Golfo Ártabro (ca. 400 km 2 ), Q Z is the vertical advective flux, and [NO − 3 ] 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).

Hydrography and Turbulent Mixing
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 −5 s −2 ), 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 −8 m 2 s −3 ).

Nutrients, Chlorophyll a, and Primary Production
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 −4 ) and relaxation (139 ± 2 µmol m −4 ) was significantly higher than during downwelling (15 ± 13 µmol m −4 ) ( 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 −2 and 1,359 ± 1,032 mg C m −2 d −1 , 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 depthintegrated chlorophyll a (65 ± 35 mg m −2 ) and primary production (5,192 ± 538 mg C m −2 d −1 ) 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 −2 ) and primary production (1,611 ± 295 mg C m −2 d −1 ) 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%. MLD corresponds to the mixed layer depth, estimated from an increase in water column density of 0.125 Kg m −3 relative to surface values. Squared Brunt Väisälä frequency (N 2 ), dissipation rates of turbulent kinetic energy (ε), and vertical diffusivity (Kz) correspond to averaged values for the 10-40 m depth interval. The same interval was used to compute vertical NO − 3 gradients and diffusive fluxes. Depth-integrated values (Int) for chlorophyll a, primary production (PP), contribution of <10 µm size-fraction, N 2 fixation rates, Gammaproteobacteria γ -24774A11 (Gamma) and UCYN-A nifH abundances, and biomass of picoplankton groups were calculated down to 40 m. BNF contribution represents the contribution of biological N 2 fixation to the supply of new N (considering N 2 fixation plus nitrate vertical diffusion). The contribution of each group to total picoplankton biomass is also included. A nonparametric one-way ANOVA (Kruskal-Wallis) was performed to test the null hypothesis that independent groups come from distributions with equal medians (*p < 0.05; **p < 0.01; ***p < 0.001). The Bonferroni multiple comparison test was applied a posteriori to analyse the differences between every pair of groups (D, U, R). n, Indicates the number of cruises included in each hydrographic condition.

Picoplankton Community Composition
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 −2 ), was during this condition slightly higher than during upwelling (54 ± 23 mg C m −2 ), and relaxation (31 ± 26 mg C m −2 ).
The vertical distribution of picoeukaryotic biomass, which potentially would include the UCYN-A diazotrophprymnesiophyte host (Thompson et al., 2012(Thompson et al., , 2014, was relatively homogeneous during downwelling, with the exception of May 2015, when the highest biomass (>3 µg C L −1 ) 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 −1 ) 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 −1 ) also occurred at surface waters.

Magnitude and Biogeochemical Relevance of Biological N 2 Fixation
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 −2 d −1 (Table 1), contrasting vertical patterns were observed during downwelling, upwelling and relaxation conditions (Figure 6).
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 −2 d −1 ), compared to upwelling (2.2 ± 3.1 mmol m −2 d −1 ), and downwelling (0.8 ± 0.9 mmol m −2 d −1 ) 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 −2 d −1 ) 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).

Magnitude and Biogeochemical Relevance of BNF
Our results demonstrate the existence of detectable BNF activity (0.001 ± 0.002-0.095 ± 0.024 µmol N m −3 d −1 ) 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 ( 15 N 2 -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 2 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 , 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 −2 d −1 ) was slightly higher than the estimations for the subtropical NE Atlantic by Mouriño- We cannot discard that our N 2 fixation rates may be underestimated (Mohr et al., 2010;Großkopf et al., 2012) due to the use of the 15 N 2 -bubble addition technique (Montoya et al., 1996) instead of the modified tracer method with 15 N 2 -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 −2 d −1 ) and during relaxation in September 2015 (77 mmol m −2 d −1 ), 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, nonupwelling 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).

Players
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 4 nifH copies L −1 ) 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 4 nifH copies L −1 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. 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(Thompson et al., , 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).
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 2 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 2 fixation rates required to account for BNF rates (see Material and Methods) ranged from 0.8 fmol N cell −1 h −1 (November 2015, 40 m) to 19.3 fmol N cell −1 h −1 (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 −1 h −1 (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.

OUTLOOK
Our study provides a novel assessment of the magnitude and seasonal variability of BNF, its contribution to the budget of new N inputs into the euphotic zone, and the abundance of diazotrophs in a NE Atlantic temperate upwelling region, subjected to contrasting hydrographic regimes. In contrast, most studies focusing on marine BNF have been conducted in the western basin, and in the (sub)tropical North Atlantic. BNF activity and presence of diazotrophs have been demonstrated in a N-rich temperate region, even though their ecological role remains unclear. This will require further spatial-temporal studies in the zone. By calculating the geometric mean of depth-integrated BNF rates estimated during the NICANOR cruises and assuming the surface area of a 3 • × 3 • grid cell (following Luo et al., 2012), we estimated an averaged new N input in the study area of about 0.4 Gg N year −1 from planktonic BNF. In the region of the Canary Current and NW African upwelling, Benavides et al. (2011) estimated an input of 0-20 Gg N year −1 , whereas Mulholland et al. (2012) estimated an input of 0.5-107 × 10 3 Gg N year −1 in the temperate NE American coast (Luo et al., 2012;Benavides and Voss, 2015). Therefore, the current biogeochemical relevance of BNF in our study system is negligible. However, the weakening of the Iberian coastal upwelling in the last decades (Pérez et al., 2010;Pardo et al., 2011), and the enhancement of stratification due to the warming in ocean surface waters might favor an increase of the relative importance of the diazotrophy in the future (Boyd and Doney, 2002;Doney, 2006). In any case, the study of the BNF processes and the responsible organisms in this region shed light on the ecophysiology of the present diazotrophic community. This first report of the presence of UCYN-A and Gammaproteobacteria γ-24774A11 in this region extends the geographic area where these diazotrophs thrive. Moreover, it provides further evidence of the co-occurrence of the UCYN-A sublineages, their wide distribution in the ocean (Farnelid et al., 2016;Martínez-Pérez et al., 2016), as well as their ability to grow in relatively cold waters with high DIN concentrations. Considering the large diversity of marine diazotrophs, and our limited knowledge of their physiology and ecology, sustained observational efforts throughout the ocean will be required to delimit accurately the magnitude and variability of marine BNF.

AUTHOR CONTRIBUTIONS
BM and EM conceived the overall study. VM, BM, EM, AF, AB, and MV designed the experiments and analyzed the data. VM collected the samples, performed the microstructure turbulence measurements and the molecular analysis. VM and AF performed the N 2 fixation measurements. VM and BM analyzed the microstructure turbulence data. AB performed the chlorophyll a and primary production measurements. BM, EM, AB, and MV supplied reagents/materials/analysis tools and equipment. VM wrote the manuscript with contributions from all coauthors.

FUNDING
This work was funded by the NICANOR project (Galician Government, EM2013/021) granted to BM and by the RADIALES project of the Instituto Español de Oceanografía (http://www.seriestemporales-ieo.com). VM was supported by a FPU predoctoral fellowship from the Spanish Ministry of Education, Culture and Sports (FPU13/01674).