The seasonal and temporal diversity of diazotrophs at Stn. ALOHA were characterized using nifH gene HTS as previously described (Turk-Kubo et al., 2015; Cabello et al., 2020). For DNA samples, seawater (~2 L) was collected from 5, 25, 45, 75, 100, 125, 150, and 175 m depths using Niskin^® bottles attached to the CTD rosette and sampled into acid-washed polycarbonate bottles. Water was immediately filtered through 25 mm diameter 0.2 µm pore size polyethersulfone filters using gentle peristaltic pumping. Filters were placed in lysis buffer AP1 (DNeasy Plant Kit, Qiagen, Germantown, MD, USA), flash frozen, and stored at -80°C until processing. DNA was extracted using the DNeasy Plant Kit (Qiagen) following the protocol detailed in Paerl et al. (2008), which includes repeated freeze-thaw cycles, agitation with a bead beater, and a proteinase K digestion to improve lysis of cyanobacterial cell walls.

Partial nifH gene sequences were amplified using a universal nested nifH PCR assay (Zehr and McReynolds, 1989; Zani et al., 2000). The first round of reactions contained 2 µL of template DNA, Platinum™ Taq DNA polymerase (6 units; Invitrogen, Carlsbad, CA), 1X PCR Buffer (-MgCl₂), 4 mM MgCl₂, 400 µM dNTP mix, and 0.5 µM of primers nifH3 and nifH4 in a reaction brought up to 20 µL with RT-PCR grade water (Applied Biosystems, Waltham, MA, USA). Second round reactions contained 2 µL of the first-round reaction, 10 units of DNA polymerase, 1X PCR Buffer, 4 mM MgCl₂, 200 µM dNTP mix, and 1 µM of nifH1 and nifH2 primers in a 15 µL reaction volume. To minimize contamination, all PCR reactions were carried out in an amplicon-free UV-hood as described by Turk-Kubo et al. (2014), and no-template-controls were processed with each set of PCR reactions. Products were amplified using the following thermocycling parameters: 95°C for 3 min, followed by 25 (round 1) or 30 cycles (round 2) of 95°C for 30s, annealing at 55°C (round 1) or 57°C (round 2), and 72°C for 45s. All samples were amplified in duplicate, screened for the correct products using gel electrophoresis, and pooled prior to sequencing. The no-template-controls did not yield any amplification. A targeted amplicon sequencing approach was used to create barcoded libraries as described in Green et al. (2015), using 5’ common sequence linkers (Moonsamy et al., 2013) on second round primers, nifH1 and nifH2. Sequence libraries were prepared at the DNA Service Facility at the University of Illinois at Chicago, and multiplexed amplicons were bidirectionally sequenced (2 x 300 bp) using the Illumina MiSeq platform at the W.M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign. Samples were multiplexed to achieve ca. 40,000 high quality merged reads per sample. Demultiplexed raw sequences are available under BioProject PRJNA913939 in the Sequence Read Archive at NCBI.

De-multiplexed raw paired end reads were merged, quality filtered (trimming reads after two consecutive bases with quality scores <20), and size selected (300-400 bp) using Paired-End reAd mergeR (Zhang et al., 2014). Primers were removed using CLC Genomics workbench (Qiagen). Chimera removal and determination of operational taxonomic units (OTUs; 97% nucleotide identity) were conducted in QIIME 1 (Caporaso et al., 2010) using UCHIME and USEARCH v6.1 (Edgar, 2010; Edgar et al., 2011). Representative sequences from OTUs with greater than 100 sequences were imported into ARB (Ludwig et al., 2004), translated into amino acid sequences, and non-nifH sequences and sequences with frameshifts were removed. After these quality control steps, a total of 4,758,763 nifH sequences remained, ranging from 5,019-33,963 sequences per sample (averaging 18,662). Samples were then rarefied to 5019 sequences using the QIIME script multiple_rarefactions_even_depth.py. OTUs were screened for known reagent contaminants (>92% nucleotide identity using the QIIME script exclude_seqs_by_blast.py) but none were found. Taxonomy was assigned via blastx of representative OTU sequences against a curated nifH database containing full length nifH sequences (www.jzehrlab.com/nifh) where nifH cluster designations have been assigned based on the convention outlined in Zehr et al. (2003). Phylogenetic trees based on partial nifH nucleotide sequences were calculated in ARB.

For the Trichodesmium and UCYN-A microdiversity analyses, OTUs representing >0.1% of the total nifH sequence dataset were included. Trichodesmium OTUs were resolved to the species level if the OTU cluster representative sequence shared >98% nucleotide (nt) identity over the full length of the nifH gene fragment to nifH from a Trichodesmium isolate.

UCYN-A1, Crocosphaera spp. (UCYN-B), Trichodesmium spp., and the DDA groups Het-1, Het-2, and Het-3 were enumerated as part of a prior study using qPCR targeting nifH (dataset doi: 10.5281/zenodo.4477269). Details on qPCR assays are summarized in Church et al. (2009), with commonly used primers/probes (Church et al., 2005a; Church et al., 2005b). Two additional NCD groups were enumerated using ddPCR as part of this study based on previously developed qPCR assays, one targeting a cluster III sequence type (Church et al., 2005a) and the other targeting γ-24774A11 (or gamma A; Moisander et al., 2008). All aspects of the ddPCR protocol, including reaction conditions, droplet generation, thermocycling parameters, thresholding, and detection limits, are detailed in Gradoville et al. (2021).

The nifH gene HTS dataset was generated from the same DNA extracts used for qPCR and ddPCR allowing for useful comparisons of the nifH gene HTS- and nifH qPCR-based approaches. For comparisons between nifH gene HTS- and qPCR-based abundances, the dominant OTU for a given target group was selected, e.g. UCYN-A1_otu0 relative abundances vs. UCYN-A1 qPCR-based abundances, after verifying the OTU had no significant mismatches to the qPCR primers/probe ( Table S1 ).

Environmental data from Stn. ALOHA ( Table S2 ) were obtained from the Hawaii Ocean Time-series Data Organization & Graphical System (HOT-DOGS; hahana.soest.hawaii.edu/hot/hot-dogs/interface.html), except for sea level anomaly (SLA) data, which are described below. Environmental data included in this analysis are from Niskin bottle collections (hahana.soest.hawaii.edu/FTP/hot/water) and included the following parameters: temperature (°C); salinity; dissolved oxygen (µmol kg^-1); pH; depth (m); particulate carbon (PC, µmol kg^-1); particulate nitrogen (PN, µmol kg^-1), particulate phosphorus (PP, nmol kg^-1); low-level nitrate+nitrite (N+N, µmol kg^-1), low-level soluble reactive phosphorus (SRP, µmol kg^-1), dissolved organic N (DON, µmol kg^-1); dissolved organic C (DOC, µmol kg^-1); total dissolved N (TDN, µmol kg^-1); heterotrophic bacteria, Prochlorococcus, Synechococcus and photosynthetic picoeukaryote abundances (x 10⁵ mL^-1); and chlorophyll a concentrations (chl a, ng L^-1). Missing data were imputed by depth via the R package bcv (https://arxiv.org/abs/0908.2062) using multivariate singular value decomposition ( Table S2 ). The mixed layer depth (MLD) for each cruise was defined by a potential density offset of 0.03 kg m^-3 relative to 10 m (de Boyer Montégut et al., 2004). For each HOT cruise, the mean potential densities at each pressure level over all CTD casts were used.

Measurements of photosynthetically active radiation (PAR) at the sea surface were conducted using a shipboard LI-COR I-1500 Light Sensor Logger (https://hahana.soest.hawaii.edu/FTP/hot/light/licor/). The resulting measurements were used to derive the mean daily-integrated PAR fluxes for each cruise. The coefficient describing the vertical attenuation of PAR (K_PAR) was computed based on measurements of downwelling PAR through the upper ~150 m using a profiling HyperPro (hp, Satlantic, Halifax, Novia Scotia) (hahana.soest.hawaii.edu/hot/hot-dogs/interface.html). Downwelling PAR was converted to scalar PAR by multiplying by 1.2 (Wozniak et al., 1992). For each cruise, we calculated the average daily integrated PAR at discrete depths (PAR_depth, mol quanta m^-2 d^-1) as the product of daily integrated incident PAR (PAR_surf, mol quanta m^-2 d^-1) and the derived PAR isopleths (% PAR), similar to Letelier et al. (2004).

Multimission satellite altimetric observations of sea level anomalies were obtained from the Copernicus Marine and Environment Monitoring Service (http://www.marine.copernicus.eu). Modifications were applied as published in Barone et al. (2019) that corrected for both the long-term linear trend and seasonal cycle; consistent with that study, the resulting SLA product is referred to here as SLA_corr.

Weighted gene correlation network analysis (WGCNA) was used to identify OTUs with highly correlated relative abundances over the time series, using the approach of Langfelder and Horvath (2008) as implemented in the WGCNA R package (v1.70.3). Of the 2916 total OTUs, the analysis included 131 non-rare OTUs (>10 reads in ≥10 samples) which accounted for 87.5% of all the reads in the rarefied dataset. Due to the observed bias between relative abundances and qPCR absolute abundances in lower euphotic zone samples (when comparisons can be made; Figure S1 ), network analysis focused on data from the 128 upper euphotic zone samples (≤ 75 m), which contained 128 non-rare OTUs (44.3% of all rarefied reads). On the log2 abundances (+1) for the 128 OTUs, we used blockwiseConsensusModules() to create a signed network with minModuleSize = 20 and power = 4, selected based an evaluation of different soft thresholds using pickSoftThreshold(). Two modules were identified which had 43 OTUs (M1) and 37 OTUs (M2). Pearson correlations between these OTUs and environmental variables were considered significant for p-values < 0.05, calculated using corPvalueStudent() with nSamples=128. These correlations were hierarchically clustered for OTUs as well as environmental variables with the pvclust R package v2.2 (correlation-based distances, average linkage clustering, 2000 bootstraps, cluster significance alpha = 0.95; Suzuki and Shimodaira, 2006).

Over the course of the study (Jun 2013 - Jul 2016), biological, chemical, and physical conditions were largely consistent with previously reported climatological records at Stn. ALOHA (Karl and Church, 2017; Table S2 ). The upper euphotic zone (≤75 m) was characterized by high-light and low-nutrient concentrations, while the lower euphotic zone (100-175 m) was characterized by lower-light and higher-nutrient concentrations. Variation in mixed layer depths (MLD) were typical for Stn. ALOHA, with deeper mixing in the winter and spring (mean of 60 m) and shoaling in the summer and fall (mean of 25 m). Over the study period, uncharacteristically shallow mixing for the winter season was observed in 2014 ( Table S2 ). SLA_corr ranged from -16 to 14 cm, with 75% of measured values falling between -4 to 6 cm.

In the upper euphotic zone (≤75 m), temperatures ranged from 22.4-27.5°C ( Table S2 ), with the highest temperatures occurring in summer-early fall (Aug-Oct). Mean daily integrated PAR at the sea surface (PAR_surf) was highest (>50 mol quanta m^-2 d^-1) in May-August, which is consistent with Stn. ALOHA historical climatology (Letelier et al., 2004). Cooler temperatures in the upper euphotic zone occurred during winter months, when PAR_surf was at a minimum (Jan, 26.3 mol quanta m^-2 d^-1). In the lower euphotic zone (100-175 m), temperatures ranged from 17.7-24.7°C.

Nutrient and chlorophyll a concentrations also followed a seasonal cycle, most pronounced in the lower euphotic zone ( Table S2 ). In the upper euphotic zone, N+N remained low year-round, ranging from 0.001-0.013 µmol kg^-3. In the lower euphotic zone, N+N ranged from 0.001-4.6 µmol kg^–1 and generally peaked in late fall to early winter, consistent with patterns described by Venrick (1988) and Letelier et al. (2004). In the lower euphotic zone, N+N concentrations were inversely correlated to SLA_corr (Pearson correlation, r(122) = -0.348, p<0.001). SRP concentrations ranged from 0.02-0.15 µmol kg^–1 in the upper euphotic zone and from 0.03-0.37 µmol kg^–1 in the lower euphotic zone. Like N+N, SRP concentrations were inversely correlated to SLA_corr (Pearson correlation, r(127) = -0.232, p<0.01) in the lower euphotic zone. The resulting N+N:SRP values were well below Redfield stoichiometry in both the upper and lower euphotic zone, averaging 0.05 and 0.65, respectively. Chl a ranged from 43-415 ng L^-1 in the upper euphotic zone and from 8-535 ng L^-1 in the lower euphotic zone, with concentrations of chl a demonstrating a consistent maximum near 100 m ( Table S2 ).

The most highly recovered nifH sequence types at Stn. ALOHA included both cyanobacteria (nifH cluster 1B) and NCDs (nifH clusters 1G and 3) ( Figure 1A ). In the rarefied dataset, 12 OTUs accounted for 79.1% of the total nifH sequences; 7 of these were cyanobacteria (59.5% of total), and 5 were NCDs (19.5%) ( Table S3 ). Five of the major OTUs have been enumerated using qPCR/ddPCR for this time series, enabling direct comparisons of patterns observed using both approaches.

The most highly recovered diazotroph was UCYN-A1_otu0 (31.9% total rarefied dataset), which had 100% nucleotide identity with the haptophyte symbiont UCYN-A1 (Candidatus Atelocyanobacterium thalassa, CP001842.1) originally discovered at Stn. ALOHA (Zehr et al., 1998). This sequence type is targeted with the UCYN-A (grpA) qPCR assay described originally in Church et al. (2005a) and used to enumerate this group at Stn. ALOHA. UCYN-A1_otu0 had highest relative abundances in the upper euphotic zone, with a subsurface maximum ( Figure 1B ) that was most prevalent in the spring and summer (Mar-Jul; Figure S2 ). These depth and seasonal patterns are broadly consistent with those determined using qPCR ( Figures S1, S3 ).

UCYN-A3 (UCYN-A3_otu5) and UCYN-A2 (UCYN-A2_otu11) sublineages were also observed, but at much lower relative abundances (2.9% and 1.0% total sequences, respectively; Table S3 ). These groups cannot currently be differentiated using qPCR that targets nifH gene fragments (Farnelid et al., 2016) thus no comparisons can be made. UCYN-A3_otu5 had peak relative abundances between 0-75 m ( Figure 1B ) but persisted in the lower euphotic zone during the fall and spring ( Figure S2 ). Relative abundances of UCYN-A3_otu5 were generally low throughout the year ( Figure 1C ) but were slightly higher in the euphotic zone between Feb-May ( Figure S2 ). Depth distribution patterns of UCYN-A2_otu11 were distinct from the other two UCYN-A sequence types, with generally low relative abundances through the upper euphotic zone and peak abundances in the lower euphotic zone ( Figure 1B ), although this is driven largely by high relative abundances at 150 m in April and October ( Figure S2 ).

Trichodesmium sp. (Tricho_otu1) was the next most abundant of all sequence-types (13.0% total sequences) and identical to T. thiebautii H9-4 within Clade I of Trichodesmium. Tricho_otu1 is enumerated by the qPCR assay described in (Church et al., 2005a; Table S1 ). The seasonality of Tricho_otu1 was distinct from that of UCYN-A1_otu0, with peaks in relative abundances beginning in the late summer through early winter (Jul-Dec, Figure 1C ), consistent with qPCR patterns ( Figure S3 ). Relative abundances of Tricho_otu1 were highest in the upper 45 m of the water column and lowest at depths ≥ 125 m ( Figure 1B ).

Additional cyanobacterial sequence types that were among the most highly recovered included three unicellular cyanobacterial groups: Crocosphaera (Croco_otu9, 2.0% total sequences), a Crocosphaera-like group (Croco_otu3, 5.9% total sequences), and the newly described spheroid body of the marine diatom Epithemia pelagica strain UMH3203 (EpSB_otu8, 2.8% total sequences; Schvarcz et al., 2022). Croco_otu3 has >99% nucleotide identity to nifH sequences reported at Stn. ALOHA (DQ088688; J. Zehr, unpublished) as well as the South Pacific (HQ229012; Moisander et al., 2014) and the Arabian Sea (JX064489; Bird & Wyman, 2012). Both Croco_otu3 and EpSB_otu8, neither of which were targeted by qPCR assays routinely used at Stn. ALOHA previously, were distinct in having higher relative abundances in the lower euphotic zone (125-175 m; Figure 1B ). This is consistent across all seasons at Stn. ALOHA; relative abundances for both groups were low in euphotic zone waters year-round, with maxima in the lower euphotic zone in the late winter, spring, and early summer (Feb-Jul; Figure S2 ). Croco_otu9, which clusters with Crocosphaera watsonii ( Figure 1A ) and is targeted by the UCYN-B qPCR assay (Church et al., 2005a), generally co-occurred with Tricho_otu1 ( Figure 1C ) consistent with qPCR patterns ( Figures S1, S3 ). However, Croco_otu9 was underrepresented in the nifH HTS dataset when compared to the qPCR dataset, where UCYN-B (grpB) nifH-based abundances were of similar magnitude to UCYN-A (grpA) and Trichodesmium ( Figure S1 ), suggesting the nifH gene HTS PCR primers used in this study do not amplify Crocosphaera nifH as well as other diazotroph groups.

All three major groups of DDAs (Het-1, Het-2, Het-3) were detected but had comparatively low relative abundances. Het-1 had the highest relative abundances of the DDAs (0.6% total sequences, Het-1_otu12), Het-3 were present at low relative abundances (0.1% total sequences, Het-3_otu27), and Het-2 sequences were recovered but not retained in the rarefied dataset (data not shown).

Although a majority of the sequences were from cyanobacterial diazotrophs, NCDs were also present. Four NCD groups affiliated with nifH cluster 1G (primarily γ-proteobacteria) had high relative abundances ( Figure 1 ). The putative γ-proteobacterium “gamma A” was recovered from nearly every sample, and a single OTU accounted for 8.5% of the total sequences (gammaA_otu2; Table S3 ). GammaA_otu2 clusters within the Marine 1 group of γ-proteobacteria (defined in Langlois et al., 2015) that also contains AO15 (Zehr et al., 1998), proteo_1 (Fong et al., 2008), γ-proteo_1 (Church et al., 2008), γ-24774A11 (Moisander et al., 2008), as well as gamma B and gamma C sequence types (Langlois et al., 2015) ( Figure 2A ). This group is targeted by the γ-24774A11 qPCR assay (Moisander et al., 2008). GammaA_otu2 was a persistent member of the diazotroph community, with high relative abundances found throughout the year ( Figure 1C ) predominantly in the upper euphotic zone ( Figure 1B ). Gamma A depth and seasonal patterns from the nifH HTS are largely consistent with qPCR-based patterns, particularly in upper euphotic zone waters ( Figures S1, S4 ).

Three additional γ-proteobacterial NCDs had high relative abundances across the dataset. 1G_otu7, which accounted for 1.8% of the total sequences, is identical to sequences recovered in the Eastern Tropical North Pacific (ETNP) Oxygen Minimum Zone (OMZ) actively transcribing nifH (e.g. KY968013.1; Jayakumar et al., 2017) and shares 99% nucleotide identity with a previously described γ-proteobacterial group Gamma 3 (Halm et al., 2012). 1G_otu6 accounted for 1.7% of the total sequences and is identical to sequences recovered from the NPSG (DQ062533.1; Zehr et al., 2007) and clusters with gamma E (Langlois et al., 2015). Finally, 1G_otu10 (0.9% total sequences) is identical to sequences reported from the South China Sea (HQ456044.1; Kong et al., 2011) and clusters with sequence types originally reported in the vicinity of Stn. ALOHA, proteo_2 (Fong et al., 2008) and γ-proteo-2 (Church et al., 2008) ( Figure 2A ). These 1G phylotypes were found throughout the water column and, except for 1G_otu7, did not have strong depth ( Figure 1B ) or seasonal patterns ( Figure S2 ).

An additional NCD affiliated with cluster III (which contains mainly anaerobes including δ-proteobacteria), cIII_otu4, was also among the most highly abundant sequence types. Accounting for 6.5% of the total nifH sequences ( Table S3 ), cIII_otu4 had 99.7% nucleotide identity to a cluster III sequence originally described and targeted with qPCR by Church et al. (2005a). In the nifH HTS dataset, this group had a distinct depth profile, with peak relative abundances in the lower euphotic zone ( Figure 1B ), but ddPCR-based abundances were generally low and did not exhibit distinct depth patterns ( Figures S1, S3 ). However, both nifH HTS and ddPCR datasets captured a strong seasonal pattern, with highest abundances in the fall (Sept-Dec; Figures 1C , S4 ).

This dataset provided a unique opportunity to characterize the nifH-based microdiversity of Trichodesmium and UCYN-A in the NPSG ( Figure 3 ). Tricho_otu1 (T. thiebautii) was the dominant Trichodesmium OTU (representing 13% of total sequences); however, >100 Trichodesmium OTUs were present in the dataset ( Table S3 ). Among these, three additional OTUs represented >0.1% of the total nifH gene sequences: Tricho_otu46 (undefined species, 0.6% of total sequences), Tricho_otu42 (T. erythraeum, 0.5% of total sequences), and Tricho_otu16 (undefined species, 0.3% of total sequences) ( Figure 3A ). These four Trichodesmium OTUs had different seasonal and depth-distribution patterns ( Figure 4 ). Tricho_otu16 was significantly, positively correlated to Tricho_otu1 (T. thiebautii) ( Figure S5 ). Relative abundances of Tricho_otu16 were highest in the fall and lowest in the spring and decreased with depth throughout the year ( Figure 4 ). In contrast, T. erythraeum (Tricho_otu42) was negatively correlated to T. thiebautii ( Figure S5 ), with highest relative abundances in the spring and summer, and often a subsurface maximum in relative abundance ( Figure 4 ). Tricho_otu1 and Tricho_otu16 sequences have no mismatches to the Trichodesmium qPCR primer/probe set routinely used at Stn. ALOHA (Church et al., 2005a). Tricho_otu46 and Tricho_otu42 have one and two mismatches to the forward primer, respectively, and both have one mismatch to the probe ( Table S1 ).

As with Trichodesmium, the dataset contained many UCYN-A OTUs (>100), however microdiversity analyses focused on the 11 UCYN-A OTUs each representing >0.1% of the total nifH gene sequences ( Figure 3A ). These included 5 OTUs with >96% nt similarity to UCYN-A1 (UCYN-A1_otu0, UCYN-A1_otu2301, UCYN-A1_otu35, UCYN-A1_otu49, and UCYN-A1_otu30) which together comprised 32.2% of total sequences, single OTUs that clustered with UCYN-A2 (UCYN-A2_otu11; 1.0% of total sequences), UCYN-A3 (UCYN-A3_otu5; 2.9% of total UCYN-A), and UCYN-A5 (UCYN-A5_otu21; 0.2% of total sequences), as well as 2 OTUs from unknown UCYN-A sublineages (UCYN-A_otu17, UCYN-A_otu18; both 0.3% of total sequences) ( Figure 3A ; Table S3 ).

Distinct seasonal and depth distribution patterns were seen within UCYN-A1 sublineages. In contrast to UCYN-A1_otu0, which had high relative abundances in euphotic zone waters in the late winter-early summer (Jan-Jul; Figure S2 ), UCYN-A1_otu30 was found to have highest relative abundances beginning in the summer and continuing through the fall (June-Dec), and a shallower maximum in the upper euphotic zone (≥50 m; Figure S6 ). These differences are reflected in correlations between relative abundances and temperature: UCYN-A1_otu30 had a positive correlation with temperature, whereas UCYN-A1_otu0 was negatively correlated ( Figure S7 ). Additionally, several UCYN-A1 OTUs had highest relative abundances in the lower euphotic zone, including UCYN-A1_otu35 and UCYN-A1_otu49 ( Figures 3B , S6 ). Notably, qPCR primers targeting UCYN-A1 enumerate all the UCYN-A1 OTUs (as well as UCYN-A_otu17 and UCYN-A_otu18; Table S1 ).

WGCNA and correlation analyses were conducted for the upper euphotic zone samples (≤75 m) to gain insight into shifts in diazotroph community composition in response to changing environmental conditions. Network analysis produced 2 modules of nifH gene OTUs based on highly correlated abundances over the time series ( Table S4 ; samples 5-75 m; 128 non-rare OTUs comprising 88% of all reads). Module 1 (M1; 43 OTUs) and was comprised of Trichodesmium (Tricho_otu1, Tricho_otu16, Tricho_otu46), Crocosphaera (Croco_otu9, Croco_otu40), in addition to several UCYN-A OTUs (UCYN-A1_otu18, UCYN-A1_otu30) and NCDs (cIII_otu4, 1G_otu7). Module 2 (M2; 37 OTUs) contained mostly UCYN-A (UCYN-A1_otu0, UCYN-A3_otu5, UCYN-A2_otu11, UCYN-A_otu17, UCYN-A1_otu21, UCYN-A1_otu49, UCYN-A1_otu504, UCYN-A1_otu2301) and γ-proteobacteria affiliated with gamma A (gammaA_otu2, gammaA_otu20, gammaA_otu31) OTUs, but also notably EpSB (EpSB_otu8). The remaining 48 OTUs had abundance patterns that correlated poorly to OTUs in M1 or M2 and were assigned by WGCNA to the catch-all module M0 ( Table S4 ). Some OTUs in M0 were rare in the upper euphotic zone but had high relative abundances in the lower euphotic zone (depths which were not included in these analyses), including Croco_otu3 and 1G_otu6.

In general, OTUs in M1 and M2 modules had opposite correlations to environmental variables ( Figure 5 ), reflecting different seasonal patterns which appear to be driven by the cyanobacterial OTUs ( Figure S8 ). Many of the OTUs in M1 had significant positive correlations to temperature, particulate matter concentrations (i.e., PC, PN, PP), and heterotrophic bacteria abundances. M1 OTUs also had significant negative correlations to SLA_corr, as well as to photosynthetic picoeukaryote (peuk) and Synechococcus abundances. M2 OTUs had largely opposite correlation patterns with significant positive correlations to SLA_corr, photosynthetic picoeukaryote and Synechococcus abundances, along with significant negative correlations to temperature and concentrations of chl a, PN, PP, and PC, as well as Prochlorococcus and heterotrophic bacteria abundances. Beyond these general trends, this analysis revealed other correlations between specific OTUs and environmental conditions. For example, in M1, both cIII_otu4 and Het-1_otu12 had significant positive correlations to chl a, reflecting the high relative abundances of these organisms in the lower euphotic zone. Likewise, in M2, UCYN-A2_otu11 correlation patterns were quite different than other M2 OTUs, with significant positive correlations to PAR_depth, K_PAR, and SRP, and inverse correlations with depth suggesting that UCYN-A2 lives in the upper sunlit layers.

Diazotroph community composition at Stn. ALOHA has been previously explored using two main approaches that target known organisms. Large, conspicuous diazotrophs, such as Trichodesmium and DDAs, have been enumerated using microscopy (Letelier and Karl, 1996; White et al., 2018). Additionally, diazotroph community composition has commonly been assessed using qPCR targeting the nifH gene from six cyanobacterial diazotrophs (Church et al., 2009; dataset doi: 10.5281/zenodo.4477269) based on their presence in several early clone libraries constructed from nifH gene transcripts (suggesting active synthesis of nitrogenase) from a small number of samples (Church et al., 2005a; Church et al., 2005b). The present study shows that using a newer approach, nifH gene HTS, enables analysis of the seasonality, depth distributions, and environmental drivers of potentially important diazotrophs that are not (or cannot be) targeted by microscopy or qPCR assays.

Abundance patterns derived from nifH gene HTS agree with those from qPCR/ddPCR-based abundance trends for UCYN-A1, Trichodesmium, gamma A, and cIII ( Figures S1, S3, S4 ), but not for other important diazotroph groups. DDAs (Het-1, Het-2, Het-3), which have long been established as quantitatively significant N₂-fixers at Stn. ALOHA (White et al., 2007; Karl et al., 2012) have low amplification efficiency using the universal nifH primers used in this study (discussed in Caputo et al., 2018). As such, seasonal and depth patterns for Het-1 and Het-3 are not consistent with qPCR-based patterns ( Figure S1 ), indicating this nifH HTS approach is not suitable for DDAs-specific inquiries. Furthermore, Crocosphaera is unexpectedly underrepresented in the nifH gene HTS, despite being detectable in the same DNA extracts using qPCR, which suggests a methodological artifact such as low PCR amplification efficiency.

The nifH gene HTS dataset also has a skew in depth distributions, that is most evident for some taxa, e.g. cIII_otu4 in samples from >100 m depths ( Figure S1 ). This likely arises from applying end-point PCR to samples having low starting nifH gene copy numbers, as might be expected in the lower euphotic zone. Although relative abundance data need to be interpreted with caution, this approach does provide qualitative information on how the overall community changes with depth and which organisms persist into the lower euphotic zone or have a preferred niche at these depths.

Furthermore, the nifH gene HTS dataset also highlights the extensive nifH diversity in this region, revealing several potentially important diazotroph taxa and potential environmental controls (discussed below). Of particular interest is the new Crocosphaera-like group (Croco_otu3) that were highly abundant in this dataset yet are not captured by current UCYN-B qPCR assays. Furthermore, this amplicon-based study allows characterization of diazotroph microdiversity of major taxa at Stn. ALOHA, including the UCYN-A/haptophyte symbiosis sublineages and Trichodesmium (discussed below). Collectively, this work illustrates the complementary role that HTS analyses can and should play in studying microbial diversity. Furthermore, coupling qPCR/ddPCR with nifH gene HTS offers insights into absolute gene abundances of specific, often relatively abundant, diazotrophs, while also providing a detailed look at diversity (including microdiversity) and relative abundances of the “entire” diazotroph community.

It is well-established that cyanobacterial diazotrophs are important in the NPSG, and Stn. ALOHA is one of the best characterized N₂-fixing ecosystems. While the nifH gene HTS dataset supports the prevalence of previously identified cyanobacterial diazotrophs (discussed above), this study provides insight on the seasonality, depth distributions, and environmental factors influencing the distribution of several diazotroph taxa that are poorly-characterized in the North Pacific.

There were diverging abundance patterns of diazotrophs at the species- or clade-level. For instance, Trichodesmium are typically enumerated at the genus level (Letelier and Karl, 1996; Church et al., 2005a; White et al., 2018), but this group includes species from multiple clades that differ in their physiology and ecology (Hynes et al., 2012; Hutchins et al., 2013) and cannot be reliably classified morphologically (Hynes et al., 2012). In the present time-series, the dominant Trichodesmium OTU (Tricho_otu1) was closely related to T. thiebautii (Clade I), agreeing with previous nifH sequencing-based studies near Stn. ALOHA (Gradoville et al., 2014; Gradoville et al., 2017b). The especially high relative abundance of Tricho_otu1 sequences from Sept-Dec 2015 corresponds to a period of high Trichodesmium cell abundances and an abnormally large number of colonies (as opposed to free filaments) reported in a microscopy-based study that used samples from the same HOT cruises (White et al., 2018), suggesting that this OTU may be linked to colony formation. Many other Trichodesmium OTUs were present, and the most abundant representatives appear to differ in their seasonality and depth distributions ( Figure 4 ). Relative abundances of Tricho_otu1 were negatively correlated with relative abundances of Tricho_otu42 (T. erythraeum) ( Figure S5 ), agreeing with previous suggestions of niche partitioning of these two species (Rouco et al., 2014). The mechanisms driving the different abundance patterns of these two species remain unknown. The vast majority of culture studies have used laboratory strains of Clade III (e.g. T. erythraeum IMS101); however, a study by Hutchins et al. (2013) showed that T. thiebautii and T. erythraeum isolates respond differently to pCO₂, suggesting that Trichodesmium clades/species respond differently to variations in the environment. More studies on the environmental drivers of Trichodesmium Clade I are needed, since T. thiebautii appears to be the dominant species in open-ocean environments (Chappell et al., 2012; Rouco et al., 2014) and perhaps globally (Marumo and Nagasawa, 1976; Carpenter and Price, 1977).

Amplicon- or metagenomic/metatranscriptomic- based studies are currently required to characterize abundance patterns that incorporate the known microdiversity of the UCYN-A/haptophyte symbiosis sublineages, which are composed of numerous sequence types that display niche and biogeographical partitioning (Turk-Kubo et al., 2017; Turk-Kubo et al., 2021). Similar to other open ocean regions, the UCYN-A symbiosis assemblage at Stn. ALOHA is primarily composed of major (globally abundant) and minor (rare) OTUs that cluster with UCYN-A1 sublineages ( Figure 3 ). The UCYN-A1 symbiosis is characteristically found in cooler waters than other diazotrophs (Moisander et al., 2010), with temperature optima <25°C. However, nifH gene HTS revealed two OTUs, UCYN-A_otu18 and UCYN-A1_otu30, that have opposite patterns, clustering with M1 ( Figure 5 ) and with highest relative abundances in waters >25°C ( Figure S7 ). The contrasting seasonal and depth distribution patterns for these minor UCYN-A1 OTUs suggest they occupy different niches and have different temperature optima. Given that the UCYN-A cyanobacterium is hypothesized to be an obligate symbiont, these patterns are likely to be driven by the ecophysiology of the prymnesiophyte hosts, which are not well understood. A recently obtained isolate of the B. bigelowii genotype III, which hosts the UCYN-A2 sublineage, is leading to important insights into B. bigelowii biology including validating its insensitivity to the presence of nitrate described in Mills et al. (2020) (Suzuki et al., 2021). However there are multiple B. bigelowii strains that have been morphologically and genetically characterized by Hagino et al. (2009) that have not yet been linked to nifH-based sublineages, limiting interpretations based solely on nifH or microscopic observations of B. bigelowii.