Conductivity-temperature-depth (CTD) measurements were made using a Sea-Bird Electronics 911-plus system affixed to a sampling rosette. Nitrate profiles were collected at 7 to 10 stations per transect with an optical nitrate sensor (SUNA V2, Sea-Bird Scientific). To create depth profiles, we aligned the SUNA and CTD data by the recorded time. Water samples from 4-6 depths at 12 representative stations were taken for direct nitrate concentration measurements to calibrate the nitrate sensor. Nitrate concentrations in water samples were measured using an Alpkem RFA continuous flow analyzer following standard colorimetric protocols (Gordon et al., 1993). SUNA nitrate profiles were calibrated by fitting a linear regression to direct measurements from corresponding depths.

To characterize the wind forcing during the expedition, we use the ERA5 reanalysis product from the European Centre for Medium Range Weather Forecasting (ECMWF; Hersbach et al., 2020). The horizontal resolution is 0.25°, and the temporal resolution is 1 hr. The 10 m wind field was used.

Surface water was collected from Niskin bottles mounted to the CTD rosette in dark, 4 L acid-washed (10% HCl) bottles (Nalgene), and promptly filtered (0.22 µm Sterivex™ filters, MilliporeSigma) using a peristaltic pump. Filters were sealed, flash-frozen in liquid N₂ and stored at -80°C until DNA extraction. To prepare samples for extraction, the Sterivex™ filters were first removed from their casing in a laminar flow hood using a clean PVC pipe cutter (1% sodium hypochlorite solution soak) to open the case and an autoclaved scalpel to separate the filter. Filters were then transferred in their entirety to sterile microcentrifuge tubes containing AP1 buffer (Qiagen) and beaten for 2 minutes with silica beads. The beads were removed using a QIAshredder spin-column (Qiagen) and DNA was immediately extracted using the DNeasy Plant Mini Kit (Qiagen) following the manufacturer’s protocol. DNA yields following extraction varied from 1267 to 4810 ng per sample.

To investigate eukaryote community diversity, we performed 18S rRNA amplicon sequencing. First, primers derived from the 18S rRNA Illumina Earth Microbiome Project were used to amplify the SSU rRNA 18S V9 marker gene (Stoeck et al., 2010) via polymerase chain reaction (PCR). The products were then gel-purified (GeneJET Gel Extraction Kit, ThermoFisher), subject to a further round of PCR to attach Illumina indices, Mag-Bead purified using AMPure XP beads (Beckman Coulter), and sequenced on Illumina MiSeq Desktop Sequencer at Old Dominion University using a 2×300 bp kit. The DNA sequences were analyzed using the DADA2 pipeline (Callahan et al., 2016); the average number of reads per sample was 66,000. Quality control steps roughly followed defaults of the DADA2 pipeline with a few minor exceptions. Reads without intact primer sequences were discarded from analysis; intact primer sequences were removed from reads using cutadapt (Martin, 2011). Following primer removal, reads were filtered and trimmed. Amplicon sequence variants were identified via BLASTn comparison (Altschul et al., 1990) to an in-house collection of eukaryote 18S rRNA sequences from the National Center for Biotechnology and Information (USA) and SILVA (German Network for Bioinformatics Infrastructure) databases. Only the subset of 18S rRNA reads that were classified as one of the major phytoplankton lineages (dinoflagellates, diatoms, or haptophytes) are presented here. Only B. bigelowii reads were enumerated to the species level.

The abundances of UCYN-A1, UCYN-A2 and B. bigelowii were determined in triplicate via qPCR. Samples were first diluted 1:10 with nuclease-free qPCR-tested ultrapure water (Millipore). Reactions were then prepared by combining 10 µl 2x TaqMan Fast Advanced Master Mix (Applied Biosystems), 4 µl nuclease-free water, 1.6 µl 5 µM custom forward and reverse primers, 0.8 µl 5 µM custom probe, and 2 µl diluted DNA extract (2.5 – 10 ng DNA) for a total reaction volume of 20 µl. These assays targeted nifH sequences (a structural gene associated with the N₂ fixation enzyme) specific to two UCYN-A sublineages, UCYN-A1 (Church et al., 2005) and UCYN-A2 (Thompson et al., 2014), and an 18S rRNA sequence associated with a B. bigelowii identified in association with UCYN-A2 (Thompson et al., 2014). Reactions were carried out using a StepOnePlus Real-Time PCR system (Applied Biosystems) following manufacturer’s guidance. Assay efficiencies for this study ranged from 92.4 to 100.4%. Samples, serially-diluted (5 × 10^-2–10^-9 ng DNA µl^-1) synthetic plasmids (GeneWiz), and no-template controls were all run in triplicate; mean values are presented here. The gene concentration range for synthetic plasmid dilution curves ranged from 9.98 × 10¹-10⁸ gene copies L^-1 for B. bigelowii 18S rRNA, and 3.56 × 10¹-10⁸ gene copies L^-1 for UCYN-A1 and UCYN-A2. Effective limits of detection and quantification were calculated as described in Selden et al. (2021) by assuming that the minimum detectable and quantifiable concentrations are 3 and 10 copies reaction^-1, respectively (Bustin et al., 2020). These limits varied due to slight differences in sample filtration volume, ranging from 170-220 gene copies L^-1 and 575-735 gene copies L^-1, respectively.

One-way analysis of variance (ANOVA) tests were performed to assess differences among transect occupations in diazotroph and host mean absolute abundances and host relative abundance among all phytoplankton sequences. Calculated test statistics were compared to those computed from 10,000 random permutations of the data to assess significance (Manly, 2007), using an alpha criterion of 0.05. This approach does not require the assumption that the residuals are normally distributed.

As the upwelling event commenced, the water column in the vicinity of the Beaufort shelfbreak was weakly stratified. Nitrate concentrations in the upper water column were low or undetectable ( Figure 2A ), and haptophytes represented only a small fraction of the eukaryotic phytoplankton present (≤1%; Supplementary Table 1 ) on the shelf and slope. The absolute abundance of B. bigelowii was low across the shelf and shelfbreak (mean = 1.44 × 10⁴ gene copies L^-1, n = 4 stations; Supplementary Table 1 ; Figure 2A ), but increased at the most seaward station (2.6 × 10⁶ gene copies L^-1). UCYN-A1 and UCYN-A2 mirrored this pattern ( Supplementary Table 1 ; Figure 2A ). As the abundance of sublineage-specific nifH genes is low compared to the 18S rRNA sequences used to quantify B. bigelowii, UCYN-A1 was only detected at the most offshore station and UCYN-A2 could not be detected at one of the mid-shelf stations ( Supplementary Table 1 ).

Upwelling along the Alaskan Beaufort Sea shelf typically occurs when easterly winds exceed 4 m s^-1 (Lin et al., 2019). This reverses the normally eastward-flowing shelfbreak jet, and brings water at depth from the basin onto the shelf (Pickart et al., 2011; Schulze and Pickart, 2012). A sub-surface mooring situated just seaward of the shelfbreak, approximately 35 km to the west of the R/V Sikuliaq repeat transect, provided context for the shipboard measurements ( Supplementary Figure 1 ). The along-coast winds became upwelling favorable just as the first transect was occupied (August 30-September 1). Following this, the near-bottom salinity increased and the shelfbreak jet reversed to the west. The second transect (September 1-2) was occupied just past the peak of the event, while the third transect (September 3-4) was carried out after the upwelling had subsided. The final transect (September 5) was done at the start of another event that was considerably weaker.

The first upwelling event caused isopycnal shoaling and drove nitrate-rich waters onto the shelf in a thin stratified layer roughly 20 m above the bottom ( Figure 2B ). Nitrate concentrations on the shelf decreased as upwelling subsided ( Figure 2C ), likely due in part to uptake by growing phytoplankton communities ( Figure 3 ) as well as deepening of isopycnals. On the final transect occupation (September 5; Figure 2D ), nitrate was again drawn onto the shelf despite the weaker wind event. This is likely because the water column had not fully relaxed from the earlier event, and therefore the previously upwelled water was still located somewhere close to the shelf edge.

Following the onset of upwelling (September 1), chlorophyll-a fluorescence increased ( Figure 3 ) and peaks broadened, suggesting enhanced phytoplankton biomass in the well-lit mixed layer. Concomitantly, we observed a marked increase in the abundance of haptophyte 18S rRNA sequences relative to the sum identified as eukaryotic phytoplankton (diatoms, dinoflagellates, haptophytes) across the transect from a mean of 1.1% to 3.0% (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 1} = 5, F = 9.34, p = 0.019), though haptophytes never made up more than 4% of the community at any given station ( Supplementary Table 1 ).

The largest sub-group among haptophyte reads were most closely aligned to sequences from B. bigelowii (>99% identity with coverage of 161/163; Hagino et al., 2013). Prior to upwelling, reads identified as B. bigelowii were only detected in the low-salinity waters (<28 psu) of the offshore end-member ( Figure 3 ; Supplementary Table 1 ), consistent with prior reports from the region (Shiozaki et al., 2018). The mean relative abundance of B. bigelowii sequences among those of haptophytes, and among those of eukaryotic phytoplankton broadly, increased at all stations during upwelling except at the most coastal station where it remained too low to be observed among 18S rRNA reads ( Supplementary Table 1 ; Figure 3 ). The highest relative abundances of B. bigelowii were observed at slope stations as upwelling subsided and at the start of the second weaker event, where it represented >1% of eukaryotic phytoplankton ( Figure 3 ) and >40% of haptophyte 18S rRNA reads. While the increase in B. bigelowii as a proportion of the haptophyte population at stations during the onset of upwelling and during peak upwelling were not statistically significant using a standard alpha criterion of 0.050 (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 1} = 5, F = 3.92, p = 0.054), likely due to our small sample size, the final occupation (at the beginning of a second, weaker event) bore significantly higher values than the occupation during the onset of upwelling (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 5} = 5, F = 6.58, p = 0.041).

The absolute abundance of B. bigelowii 18S rRNA reads correlated with its relative abundance among eukaryotic phytoplankton ( Supplementary Figure 2 ). During upwelling, the absolute abundance of B. bigelowii and both UCYN-A sublineages increased at all stations except the shallowest, where abundance of B. bigelowii and UCYN-A2 decreased ( Figure 2 ; Supplementary Table 1 ). Nitrate was undetectable in the water column at this station ( Figure 2B ). B. bigelowii absolute abundance across the entire transect was greater during upwelling than pre-upwelling (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 1} = 5, F = 9.16, p = 0.046). Significant differences between pre-upwelling and upwelling conditions were not observed for UCYN-A1 (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 2} = 5, F = 4.51, p = 0.094) or UCYN-A2 abundance (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 1} = 5, F = 5.03, p = 0.086). As upwelling began to relax, a clear gradient in the abundance of the three organisms appeared with the highest abundances observed at the most offshore station ( Figure 2C ). By the fourth occupation (at the onset of a weaker upwelling event), abundances of UCYN-A1 (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 5} = 5, F = 4.00, p = 0.049) and UCYN-A2 (one-way ANOVA, n_{Aug. 30} = 5, n_{Sep. 5} = 5, F = 3.90, p = 0.009) had significantly increased relative to pre-upwelling conditions.

Given the Eulerian nature of this study, and the physical dynamics of the region, some of the observed enhancement in B. bigelowii and UCYN-A abundances may be attributable to biological patchiness in shelf-waters. However, it is highly unlikely that the observed enhancement resulted solely from alongshore transport of shelf waters because B. bigelowii reads were too rare to be detected among 18S rRNA amplicon sequences collected on the shelf to the east of the study site at the onset of upwelling (Einarsson, 2021). Indeed, they were only detected to the east in waters seaward of the shelf-break, as was the case for the transect described in this study.

As a thought experiment, we calculated doubling times for sequence abundance at each station ( Supplementary Table 1 ). These values were largely reasonable when compared to doubling times for cells (1+ days)—a comparison which assumes that the rate of DNA replication is equal to the rate of cell division. However, B. bigelowii 18S rRNA doubling times at the three intermediate shelf stations during the second occupation (upwelling; September 1-2), and the most coastal station at third occupation (relaxation; September 3-4), showed extreme rates of increase which were unrealistic when compared to doubling times for cells. This trend can be explained by changes in ploidy. Many eukaryotes undergo cyclic changes in ploidy related to their growth cycles and cell size (Parfrey et al., 2008). High levels of DNA are observed in large vegetative cells (Kondorosi et al., 2000). Increases in haptophyte ploidy are thus consistent with the hypothesis that upwelling relieved growth limitation of the diazotroph-haptophyte symbiosis. We note that UCYN-A also displayed rapid DNA doubling times at certain stations; recent work indicates that some cyanobacteria do contain multiple copies of their genome and that this number can greatly vary at different growth stages (Griese et al., 2011; Oliverio and Katz, 2014). As discussed above, community patchiness, too, represents a confounding variable in these calculations.

Throughout the study, the abundances of B. bigelowii 18S rRNA and UCYN-A2 nifH sequences remained tightly correlated at a ratio of ~5 18S rRNA:1 nifH gene copy (linear regression, R² = 0.966, p < 10^-11). In conjunction with recent work from Mills et al. (2020), who found that B. bigelowii relies on diazotroph-derived N to meet its N demand even when other forms of N are available, this finding suggests that UCYN-A was continuously fixing N for its host. It is thus reasonable to expect that volumetric N₂ fixation rates increase with abundance of the symbiosis.

Taken together, the results presented here indicate that upwelling favored growth of UCYN-A and its host, B. bigelowii, in surface waters across the Beaufort shelfbreak. In so doing, the upwelling event likely drove increased N₂ fixation as well. We propose two potential explanatory mechanisms: (1) the growth of UCYN-A and its symbiont were indirectly stimulated by nitrate (e.g., by increased availability of plankton exudates) as reported from the California Current System by Mills et al. (2020), and/or (2) the UCYN-A/haptophyte symbiosis proliferated because non-diazotrophic phytoplankton drew down nitrate in advance of other upwelled nutrients and the community became N-limited while phosphate and iron remained available, creating a niche for the diazotrophic symbiosis. As samples for DNA analysis were only collected from surface waters, where nitrate concentrations were always low due to rapid biological utilization, and because samples for phosphate determination were not collected, we unfortunately cannot speculate as to mechanism at this time.

Regardless of the exact mechanism, the response of the UCYN-A/haptophyte symbiosis to upwelling observed here may be occurring in upwelling systems globally. Subramaniam et al. (2013) observed enhanced N₂ fixation rates following upwelling in the equatorial Atlantic where UCYN-A is known to be highly active (Martínez-Perez et al., 2016). Additionally, Sohm et al. (2011) observed high N₂ fixation rates in the Benguela upwelling system. Though they assessed the abundance of a range of diazotrophs, these authors were unable to identify the diazotroph responsible. However, the primers now commonly used for UCYN-A2 and its host had not yet been developed at the time of their study.

Finally, these results suggest that Arctic Ocean N₂ fixation may be highly variable as local diazotrophs respond to dynamic physical conditions. This variability is likely not reflected in the only published estimate of N₂ fixation on the Arctic shelf (3.5 Tg N y^-1 or ~2.7% of global N₂ fixation; Sipler et al., 2017), which is based on relatively few data points with low spatial and temporal coverage. Considering upwelling-driven N₂ fixation could thus significantly augment basin-wide fixed-N input estimates in the Arctic, as well as in other upwelling systems.