Expeditions to the WTNA were conducted May 23–June 23, 2010 and September 4–October 6, 2011 on board the R/V Knorr and R/V Melville, respectively. Figure 1 shows the cruise tracks and station locations for sample collections during the field expeditions. Sampling for nucleic acids were limited to the second leg (September 9–October 6, 2011) of 2011. Seawater was collected from discrete depths using Niskin bottles (12 or 24 L) arranged on a Conductivity Temperature Depth (CTD) rosette. The suite of parameters recorded at the time of nucleic acid collection and used in the statistical analyses (see below) is summarized in Supplementary Table S1.

Cells for microscopy observations were collected directly from the CTD at the same or hydrographically similar depths as water collected for DNA samples, although at different times. The microscopy samples were always collected at local noon, while the CTD casts for the nucleic acid samples varied from early morning until late afternoon. The entire contents of the Niskin were gravity filtered onto a 47 mm diameter Poretics (Millipore) membrane filter with a pore size of 10 μm. Gravity-filtration time varied from 30 min to 2 h. If the filter clogged by 2 h, the remaining volume in the Niskin was measured using a graduated cylinder, and the volume filtered was noted. The filter was mounted onto an oversize microscope slide (75 mm × 50 mm × 1 mm) and examined at 400X under a Zeiss Axioskop Epifluorescence microscope (Zeiss, Berlin, Germany). Phycoerythrin and chlorophyll a (Chl a) was used to identify the symbiotic Richelia and the phytoplankton community by epifluorescence microscopy using green (510–560 nm) and blue (450–490) nm excitation wavelengths, respectively. The hosts of the DDAs were identified based on cell morphology. Qualitative observations of cell integrity for the DDAs, and general composition of the phytoplankton community were recorded.

Nucleic acids were extracted using the modified DNAeasy plant kit (Qiagen) method described in Moisander et al. (2008) and the final elution volume was 70 μL. The Richelia specific TaqMan (Applied Biosystems) PCR primers and probes described in Church et al. (2005) and Foster et al. (2007) were used to quantify the abundance of het-1 and het-2. The TaqMan PCRs were conducted in a GeneAmp 9700 (Applied Biosystems) or MiniOpticon (Biorad) sequence detection systems with assays and conditions as previously described (Foster et al., 2007). Briefly, the following parameters apply: 50°C for 2 min, 95°C for 10 min, and 45 cycles of 95°C for 15 s followed by 60°C for 1 min. Reactions were performed in quadruplicate with the fourth replicate used to estimate the reaction efficiency (see below). Two microliter of 5 kDa filtered nuclease free water was used for the no template controls (NTCs). No nifH copies for either het group were detected in the NTCs. The detection limit of the assay is between 1 and 10 nifH copies.

Gene copy abundances were calculated from the mean C_t value of the three replicates and the standard curve for the appropriate primer and probe set (see below). Samples where one or two of the three replicates produced an amplification signal, were noted as detectable, but not quantifiable (dnq). The nifH gene copy abundances have been normalized to copy numbers per cell in the Richelia trichome by dividing by 5 and 4 for het-1 and het-2, respectively, as these were the number of cells observed by microscopy in the respective symbiont’s trichome.

For each primer and probe set, duplicate standard curves were made from 10-fold dilution series ranging from 1 to 10⁸ gene copies per reaction. The standard curves were made from linearized plasmids containing the nifH gene. Regression analyses of the number of cycles were analyzed in Excel. The PCR efficiency for each sample was determined as previously described by Short et al. (2004).

Multivariate analysis in the form of piecewise SEM was conducted in the software R¹ using the packages ‘NLME,’ ‘lavaan,’ and ‘piecewiseSEM’ (Lefcheck, 2016). An initial spearman’s rank correlation on the qPCR data was performed using IBM SPSS (ver. 23), which covered correlations of all parameters measured in the system and investigated in the piecewise SEM model.

The results of the spearman rank correlation combined with previous knowledge on the investigated system were used to formulate hypotheses on pathways of interaction between parameters in the model, where all parameters could act as both predictor and response variables. The hypotheses acted as a framework when designing and optimizing the piecewise SEM to fit the generated data. All measured parameters were included in the framework, but not all of them, due to statistical non-significance (p > 0.05), were still part of the model after optimization. Additionally, station and depth were considered random factor variables in the model.

Built into the ‘piecewise SEM’ package is a ‘missing pathways’ command, which provides information, supported by statistical significance to improve the SEM model. Thus, the piecewise SEM model is stepwise and modified to account for missing or incomplete pathways between response and predictor variables for each general linear mixed effects model that makes up the piecewise SEM model. The missing pathways were then tested in parallel with evaluation of the Akaike Information Criterion (AIC), which estimates the robustness of the current model compared to other models of the same dataset, but generated with different pathways (Shipley, 2013). If the tested pathways were statistically significant (p < 0.05) and generated a lower AIC score, they were included and the model was further optimized to account for more variation. Subsequently, non-significant parameters were excluded from the model. It is important to note that AIC in itself does not explain variation but is rather a means to compare models derived from the same dataset.

When the final optimized model was procured (lowest AIC score with the most variables included) chi-square statistics were run to evaluate the model goodness-of-fit (Shipley, 2009; Lefcheck, 2016). If the chi-square was statistically non-significant (p > 0.05) the model was a good fit to the data. In addition to the goodness-of-fit for the model as a whole, an R² value was calculated for each general linear mixed effects model. The R² value is a measure of the variation in the data explained by a general linear mixed effects model for a particular pathway.

Finally, in addition to the direction of the effects found between variables in the piecewise SEM model, a pathway table, with estimate values for all direct effects was reported. The effect estimates were used to calculate and compare the strengths of direct and indirect effects between variables in the system. Indirect effects were described as a predictor variable (Y) having an effect on a response variable (X) through a simultaneous response and predictor variable (Z) and is represented by the following:

For example:

The arrows (1–2) represent the effect estimate values that are multiplied (3) to procure the strength and direction (positive or negative) of the indirect effect of predictor variable Y on response variable X. In equation 3, the numbers are derived from SEM of 2010 dataset.

Following an earlier classification scheme based on sea surface salinities (SSS), three salinity categories were defined: ‘low’ SSS containing stations with SSS < 30, ‘mesohaline’ stations have SSS between 30 and 35, and stations with SSS > 35 were classified as ‘oceanic’ (Foster et al., 2007). A summary of the stations and nutrient concentrations for the three categories is provided in Table 1, and the corresponding station locations and surface nutrient concentrations are shown in Figures 1, 2, respectively.

Dissolved inorganic phosphate and Si in the surface waters (<5 m) were inversely correlated with SSS in 2010 and 2011 (Supplementary Figures S1A,B). The relationship between DIN and DIP (N:P) in the surface waters of both cruises showed deviations from a 16:1 relationship (i.e., Redfield ratio), and the N:P in 2010 compared to 2011 were indicative of an increased N limitation (Supplementary Figures S1C,D). Moreover, the quasi-conservative tracer N^∗, which can be used to estimate the distribution of N₂ fixation (Gruber and Sarmiento, 1997), showed more positive values in areas of low DIN and DIP for both cruises (Figure 2).

Lower SSS were common in the western region (46° – 51° W, 5° – 8° N) of both cruises (Figure 2); one exception was station 23 of 2010 (54° W, 10° N). Warmer sea surface temperatures (SST) coincided with the lower SSS in 2010, whereas in 2011, higher SST (≈30°C) was measured to the north and northeast. In 2010, the low salinity lens varied 5–15 m and was enriched in DIP (0.21–0.59 μmols L^-1) and Si (17.1–56.6 μmols L^-1), however, DIN was often below detection (bd). In 2011, the low SSS lens penetrated the upper 12 m, and had a similar dissolved Si concentrations (19.3–42.1 μmols L^-1) and larger range in DIP (0.18–1.02 μmols L^-1) than in 2010. Similar to 2010, DIN was bd in the surface lens at the low SSS stations of 2011.

The mesohaline surface salinity stations of 2010 were localized in the northeast region (50° – 56° W, 5° – 11° N), while in 2011, mesohaline SSS was observed more broadly (Figure 2). Surface dissolved nutrient concentrations at the mesohaline SSS stations during both cruises were characterized as low in DIP (bd-0.60 μmols L^-1), modest concentrations of Si (2.19–17.7 μmols L^-1), and DIN was often bd.

Stations with surface salinities typical of oceanic conditions in 2010 and 2011 were located in the eastern region during both cruises and in addition to the north in 2011 (Figure 2). DIP and DIN were bd at most stations of 2010 with the exception of stations 6, 7, 20, and 27 where DIP measured 0.1–0.7 μmols L^-1 and at station 27 the DIN measured 0.02 μmol L^-1. Dissolved Si was measureable at all oceanic stations of 2010 and 2011 and ranged 0.78–5.99 μmol L^-1. DIP and DIN were higher in concentration in the oceanic stations of 2011, where DIP measured as high as 0.26 μmol L^-1 and DIN measured 0.11 μmol L^-1 at station 26.

There were several low SSS stations that experienced fluctuations in their SSS and dissolved nutrient concentrations increased, decreased or remained unchanged while on station. An example was at stations 9 and 16 in 2010 and station 13 in 2011 where the SSS decreased by 3–6 PSU and DIP and Si each increased up to a twofold higher nutrient concentration. Only at station 9 did the DIN increase from bd to 0.27 μmol L^-1.

The low SSS stations of 2010 (stations 4, 10, 14, 16, and 23) were composed of non-symbiotic diatoms; in particular, the surface sample at station 4 was largely composed of Coscinodiscus spp. (Supplementary Figure S2A), while observations from depth (21.6, 41.9, and 71.5 m) found higher densities of Pseudo-nitzchia spp. and Skeletonema spp. diatoms. The latter were also common at stations 10, 14, and 16 of 2010. Station 23, another low SSS, was dominated by non-symbiotic Chaetoceros spp. diatoms, and at 10 m, the high densities of Chaetoceros spp. were observed exuding a gelatinous matrix. Below the freshwater lens at several of the low SSS stations of 2010 (e.g., stations 4, 16, 23), symbiotic H. hauckii and R. clevei were observed.

Similar to the qPCR results (see below), high densities of the H. hauckii-Richelia symbioses were observed and limited to the meso and high SSS stations of 2010. Although highly abundant, the cell integrity of both symbiont and host varied greatly, such that both degraded cells (e.g., symbiont comprised of only a terminal heterocyst, host with reduced chloroplast biovolume or empty frustules) were observed mixed with fully ‘intact’ chain forming cells (e.g., 8–12 cells and < 50 cells) or recently divided cells (e.g., pairs of cells). Degraded or moribund looking H. hauckii symbioses were more common below the depth of maximum abundance and usually coincided with a higher abundance of R. clevei-Richelia symbioses with good cell integrity. The latter was most obvious at stations 5, 6, 7, and 20. Consistent with a poor detection by qPCR, the R. clevei-Richelia symbiosis was rarely observed during the 2011 cruise. The H. membranaceus-Richelia symbioses was less common on the 2010 cruise with the exception of station 1 where short chains were observed and the Richelia had no vegetative cells (Supplementary Figure S2B), while on the 2011 cruise it was the dominant DDA albeit at low densities.

Cells similar in size and morphology to Crocosphaera watsonii (hereafter referred to as C. watsonii-like cells) were infrequently observed at low densities in 2010 (station 4, 5, 6, 7, 17) and in 2011 (station 13). At several stations, and typically at depths below 30 m, Trichodesmium spp. colonies and/or free-filaments were observed in a state of degradation reminiscent of viral lysis or as short trichomes (e.g., hormogonia; for example, at stations 4, 6, 7, 16, 20, 25, 27 in 2010). A few observations (stations 19, 20 in 2010) were noted of large Trichodesmium colonies embedded with radiolaria, dinoflagellates, diatoms, or other plankton, including DDAs (Supplementary Figure S2C). Trichodesmium was less common in 2011 and observed as mostly as free filaments.

Picocyanobacteria, or cells with similar morphology and cell diameters (>0.2 mm) to Synechococcus, and pico-eukaryotes were also observed, despite their smaller cell diameters than the pore size of the collection filter. Both populations were observed as free cells, attached to marine snow particles, as components of fecal pellets, and colonial picocyanobacteria were also frequently observed at depth (Supplementary Figures S2D,E). Densities were higher for the pico-plankton in the deeper depths (e.g., DCM), and in particular at station 20 of 2010 it appeared to be ‘rainfall’ of Synechococcus (Supplementary Figure S2F).

In 2010, the het-2 was highly abundant (detected in 76%, 73 of 96 samples) at all stations sampled with the exception of the low SSS stations (Figure 2 and Supplementary Table S1). Modest densities (<10² nifH copies L^-1) for het-2 were detected at low SSS stations 4 and 16, but at 2 depths below the low SSS lens. Highest densities for het-2 were at the mesohaline followed by oceanic SSS stations (Figure 2). Although het-1 was present and co-varied in depth distribution to het-2, it was less abundant (bd-10⁵ nifH copies L^-1) and only represented 0.06% of the total nifH copies detected compared to the 99.94% representative of het-2. Het-1 consistently co-occurred with het-2, and unlike het-2, the abundances were higher at the oceanic SSS compared to the mesohaline SSS stations. Interestingly, the average depth of maximum abundance for het-1 and het-2 were consistently parallel, and shallower at the oceanic SSS stations (8 m) compared to the mesohaline and low SSS stations (20.2 and 20 m, respectively).

In 2011, het-1 abundance was bd in all samples assayed. Het-2 was present, although at a lower detection (40%, 23 of 57 samples assayed) and limited to the mesohaline SSS stations; maximum abundances never exceeded 10³ nifH copies L^-1. Given the limited number of samples taken (e.g., one low and one high SSS) and the low detection of het-2 in 2011, our comparison to 2010 is restricted. Unlike the 2010 expedition, the depth of maximum abundance varied greatly for het-2, from the surface (station 25) to 70 m (station 21) (Supplementary Table S1).

Given the primary interest of applying the SEM model was to predict the influence a parameter or sets of conditions have on DDA abundances, we summarize only variables that were significantly directly/indirectly affecting het-1 or het-2 abundance. Thus, we exclude the majority of influence of the environmental parameters on each other and any cascading effects in the system.

Overall, the SEMs were adequate fits to the data in both 2010 (χ²= 58.95, p = 0.30) and 2011 (χ² = 68.11, p = 0.54) (Figure 3). Turbidity was the most significant parameter (-0.76, p < 0.0005), directly influencing the abundance of het-2 in 2010 (Figure 3A). Moreover turbidity mediated the indirect and negative effects of a number of other parameters (DIN, Si, fluorescence, salinity and temperature: -0.28, p < 0.001, -0.27, p < 0.0001, -0.75, p < 0.0001, -0.47, p < 0.0001, -0.13, p < 0.001, respectively) on het-2 abundance and in addition, had an indirect negative effect on het-1 (-0.60, p < 0.0005) mediated through het-2 abundance. Het-1 abundance was directly affected by both salinity (0.14, p < 0.05) and DIP concentration (0.28, p < 0.0001). One of the strongest relationships identified in the 2010 dataset was the positive direct effect het-2 abundance had on het-1 abundance (0.79, p < 0.0001).

The SEM model applied to the 2011 dataset identified a different set of conditions and decreased number of paths influencing the het-2 abundance (Figure 3B). DIP was the most prominent factor positively (0.50, p < 0.02) influencing het-2 abundance. Moreover, DIP mediated the indirect and negative effects of temperature (-0.46, p < 0.0001) and salinity (-0.35, p < 0.0001) on het-2 abundance. Nitrite concentration was also indirectly, but positively influencing het-2 abundance (0.16, p < 0.005) and was mediated through DIP. Beam transmission was identified as a positive variable on het-2 abundance in 2011 although the relationship was direct it was less robust (0.29, p < 0.03).

Similar to earlier work in the WTNA, the influence and extent of the AR plume was distinguishable by decreased surface salinities during both cruises. Monthly MODIS (or Moderate Resolution Imaging Spectroradiometer) composites showed a broad influence on surface chl a and entrainment of the plume in the North Brazil current (NBC) and flow northwest toward the Caribbean Sea in 2010, while the composite from October 2011 showed evidence of the seasonal and narrow retroflection (Muller-Karger et al., 1988) at approximately 10°N where surface waters flowed easterly (Supplementary Figure S3). Consistent with earlier investigations is the immediate increase in chl a, and subsequent distinct community composition of a non-symbiotic diatom bloom in the AR plume waters. Whereas the surface abundances of both DDAs in 2010 appeared to increase with the plume’s entrainment northwest, in 2011, DDAs rarely exceeded background densities (e.g., <10² nifH copies L^-1). Moreover, only het-2 was detected by qPCR in 2011 and appeared less affected by the plume’s directional flow or the plume was overall less influential on the DDA abundances.

The dissolved nutrient concentrations measured during our expeditions were similar to previous investigations in the WTNA (Carpenter et al., 1999; Foster et al., 2007; Subramaniam et al., 2008; Goebel et al., 2010). Surface DIN was not correlated with SSS for either cruise, and was often below detection. Moreover, in 2010, DIN deficits were also coincident with low maximum quantum yield measurements (e.g., Fv/Fm) indicative of nutrient stress, including measurements made at the low SSS, or plume stations and in areas of high DDA densities (e.g., Fv/Fm = 0-0.2) (Goes et al., 2014). Given that DIN is not limiting for the DDAs due to the diazotrophic symbiont, other stressors or limitation must have influence over the DDA cell abundances and cell activity. Moreover, our microscopy observations of cell integrity of both symbiont and host were also indicative of non-optimal conditions at many of the high DDA abundance stations.

Both DIP and Si in the surface waters were inversely correlated with SSS during both expeditions suggesting the AR as major source of DIP and Si. Similar relationships between SSS and Si have been reported in the surface waters near Barbados in the Caribbean Sea (Steven and Brooks, 1972) and near Puerto Rico (Froelich et al., 1978). Despite measurably high concentrations of Si even in the low SSS lens, Goes et al. (2014) claimed Si deficiency throughout the 2010 cruise, especially upstream of the plume which was coincident with lower Fv/Fm measures and maximum abundance of non-symbiotic diatoms, many of which possess large cell diameters (Coscinodiscus, Chaetoceros) and thus higher Si requirements. A higher expression of silicon transporters was also identified in a metatranscriptomic study of the 2010 expedition (Zelinski et al., 2016) in the outer plume, which was coincident with the higher DDA bloom densities. Despite similar surface Si concentrations in 2011 and slightly elevated DIP concentrations, neither DDA reached high densities in 2011.

Others have reported evidence of a fresh water lens distant from the AR mouth with unique chemical and phytoplankton signatures (Borstad, 1982a,b) such as the early observations made in the Guinea Current (10°N) of the NBC (Ryther et al., 1967; Cochrane, 1969; Gibbs, 1970). These lenses, which are sometimes referred to as meanders, originate from the AR and carry distinct phytoplankton populations (Borstad, 1982a,b). The meanders can be 300 km in diameter and 20 to 50 m deep with salinity anomalies up to 9 PSU (Flagg et al., 1986). Notably, several stations in both 2010 and 2011 that were occupied for longer periods of time had considerable changes in salinity and dissolved nutrient concentrations within a short period of time. For example, stations 9 and 16 in 2010 experienced decreases in the SSS (4-6 PSU) with relatively sharp increases in dissolved nutrients (e.g., twofold increase in DIP and Si) within hours (e.g., 2.5 and 7 h, respectively) while abundances of Richelia remained relatively stable below the surface (8–16 m). However, in the immediate surface at these two stations, the abundance of Richelia associated with H. hauckii decreased strongly from 1.27 × 10⁴ nifH copies L^-1 to near the detection limit (dnq). Both stations (9 and 16) are near one another (∼898 km) and on a northwest trajectory from the AR source. The decreasing salinity and increasing nutrients suggests a possible influence of a less saline AR meander, however, it appears that the near surface symbiotic diatom populations were either displaced by the meander and/or negatively influenced by its presence.

In 2010, warmer SST were localized to the stations of low SSS, which is also indicative of the AR influence. In contrast, the increased SST observed in 2011 was localized to oceanic stations far from the AR and therefore more suggestive of higher SST caused by increased stratification common to the open ocean. Moreover, the stations in 2011 with the highest SST also had the least dissolved nutrients (DIN, DIP, and Si) also indicative of decreased mixing from below and strong stratification. Thus, hydrographic conditions in the WTNA largely derived from the AR, favor DDAs, and in particular het-2 and Hemiaulus spp. hosts.

Abundances of het-1 and het-2 from both cruises corroborate and complement a number of earlier studies in the North Atlantic (Villareal, 1994; Carpenter et al., 1999; Foster et al., 2007; Subramaniam et al., 2008; Goebel et al., 2010; Ratten et al., 2015). Although collections were from different years and few studies use similar units for reporting abundances, the combined earlier datasets represent a near complete yearly investigation of DDA abundance and distribution in the North Atlantic, and also underscore the widespread and persistent distribution of DDAs, in particular the Hemialus-Richelia symbiosis throughout the year in the WTNA (Supplementary Table S3).

In spite of different seasonal occurrences, the blooms of the H. hauckii-Richelia symbioses in the WTNA have a few common observed hydrographic features that seem to prime the region for high densities. For example, highest abundances are consistently recorded at mesohaline SSS, while the salinity and temperature at the depth of maximum abundance vary greatly from 31.4 to 36.17 PSU (mesohaline-oceanic) and 25.7–28.5°C, respectively. High densities of het-1 and het-2 have been reported in other less oceanic habitats with similar salinity and temperature optima, e.g., New Caledonia coast (Turk-Kubo et al., 2015); North American Coast (Ratten et al., 2015); Congo and Niger Plumes (Foster et al., 2009); S. China Sea (Bombar et al., 2011). Wider ranges in temperature optima appear to be more the norm (Moisander et al., 2010), however, unlike the het groups, other diazotrophs, e.g., Trichodesmium spp., have less of a range in salinity optima (e.g., 34–36 PSU).

Most blooms in the WTNA have been reported from early spring and are within the same range as the densities reported here in May–June of 2010. Moreover, maximum abundances in mesohaline SSS are often in the subsurface (20–25 m), while depth maxima for DDAs appear in the near surface at the oceanic stations, potentially favored by higher light intensities. Our correlation matrix (Spearman’s rho) also identified several highly significant (p < 0.001) relationships between het-2 (and het-1) abundance favoring higher light intensity and quality, including PAR and beam transmission (Supplementary Table S2). Similar correlations between PAR and het-1 and het-2 abundance was determined in a dataset from the SW Pacific in the Melanesian archipelago, where het-1 maximum abundance (15 m) was similar to that reported here in the subsurface (Stenegren et al., unpublished).

Waters depleted in DIN, but measurable DIP (0.04–0.13) and dissolved Si ranging from 1.0 to 10.7 μmol L^-1 also appears to favor high abundances of DDAs in the WTNA (Supplementary Table S3). The higher range of DIP and Si is consistent with expected higher seasonal outflow of the AR, and is consistent with reports from early spring (Supplementary Table S3). In fact, the 2010 bloom reported here occurred during the maximum discharge of the AR for 2010 (Yeung et al., 2012). The low to background densities of the DDAs in 2011, when nutrient concentrations were seemingly similar to 2010, suggests that there is a seasonal pre-requisite for early dry season, however, similarly high abundances (>10⁶ cells m²) for the H. hauckii-Richelia symbioses have been reported in winter (Feb 2001) and fall (October–November 1996) (Carpenter et al., 1999; Subramaniam et al., 2008). Therefore, it seems that the high abundances for the DDAs, and Hemiaulus-Richelia in particular are not directly linked to a season or single environmental condition, but rather large ranges in salinities (mesohaline to oceanic), temperatures, irradiances, and nutrient concentrations. The latter highlights the difficulty in identifying a single parameter that drives DDA abundance.

Draft genomes for two of the Richelia strains (het-2) which associate with Hemiaulus spp. and the Richelia symbiont of R. clevei (het-1) have been sequenced and their genomic capacity indicates that Richelia cannot utilize a number of DIN sources (e.g., nitrate, nitrite, urea) and lacks ammonium and urea transporters (Hilton et al., 2013). The het-2 strains, in addition, lack alkaline phosphatase, phoA (Foster, unpublished). Alkaline phosphatase is used to cleave free phosphate groups through the hydrolysis of phosphomonoesters (Kuenzler and Perras, 1965) and is indicative of P stress. A closer comparison of the genes for P acquisition in the het strains identifies a number of disparities. For example, the het-1 symbiont draft genome contains in addition to phoA, several genes for dissolved organic P (DOP) acquisition and transport (e.g., phnD, phosphonate binding; phoD, phosphodiesterase; ppK, ppX polyphosphate kinase). Thus, the genome content of the Richelia het-2 symbiont requires measurable DIP, while the Richelia het-1 has a broader genetic repertoire for scavenging P and, perhaps can alleviate P limitation by utilizing DOP. Although we have no information on DOP concentrations, higher abundances of het-1 were enumerated and perhaps favored at the oceanic stations, areas with the lowest DIP and DIN concentrations. Evidence of DOP positively influencing abundance, and nifH gene transcription by het-1 was recently reported in a mesocosm experiment in the subtropical eastern North Atlantic (Meyer et al., 2016). It should be noted that the Richelia genomes are draft, and therefore lacking a gene or pathway could be a result of incomplete genome sequencing. The differences in the het strains genome content does highlight how the DDA symbionts differ in potential resources which could easily drive their distribution, abundances, and niche adaptation. Similar adaptations have been identified in other diazotrophs (Dyrhman and Haley, 2006).

Piecewise SEM (hereafter referred to as SEM) is still new to the field (Lefcheck, 2016), however, it has been successfully applied in a diverse suite of studies to determine the links between ecosystem multifunctionality and biodiversity, the role of biodiversity on top-down controls, and the role of multitrophic functional diversity on estuarine ecosystem function (Duffy et al., 2015; Jing et al., 2015; Lefcheck and Duffy, 2015). SEMs are often referred to as probabilistic models where a number of predictor and response variables are connected within a network of direct, indirect and cascading effects (Lefcheck, 2016). Moreover, the SEM differ from more traditional model approaches by incorporating hypothesized causality between variables, and allowing variables to be both predictor and response variables in the same model (Lefcheck, 2016). The causality is based on previous knowledge of the investigated system, meaning that the SEM is built on informed choices and hypotheses.

Historically, blooms of the Hemiaulus-Richelia symbioses are common in the WTNA, and correlated with seasonal discharge and SSS (Subramaniam et al., 2008). Thus, we expected to find correlations between het-2 abundance and variables indicative of the AR. Second, we hypothesized that high abundances of one DDA, e.g., het-2 as observed in 2010, would have a negative impact on the other DDA since resources could have become limiting or competitive. Finally, we hypothesized that there are conditions or sets of conditions that favor one DDA over the other given that host diatoms differ (genus, size), as do the symbiont strains (e.g., trichome length, location, diversity and genome content).

Contrary to our predictions and previous findings that DDA abundance, in particular H. hauckii-Richelia symbioses, is largely influenced by salinity and localized to regions of fast DIN utilization relative to DIP, the SEM did not find significant direct interaction between salinity or the latter nutrients and het-2 abundance. Despite high abundances quantified by qPCR, our observations of cell integrity for the H. hauckii-Richelia symbioses suggested that as a whole, the populations appeared senescent. The latter may explain the poor correlation between salinity and/or nutrients and het-2 abundance. In fact, most often mixed populations of moribund and healthy cells were observed. However, if we consider the indirect effects of nutrients and salinity, the SEM model did in fact identify that both parameters were indirectly influencing the het-2 abundance mediated through turbidity. Moreover, salinity and DIP directly influenced het-1 abundance. Thus, the SEM model supported our initial hypothesis that the AR is a predictor of DDA abundance.

Turbidity is also an important parameter influencing species abundance and activity and related to the AR plume. Likewise, turbidity was identified in the 2010 SEM as a robust negative parameter directly and indirectly influencing het-2 and het-1 abundances, respectively. The AR discharge is equally high in particulates and colored dissolved organic matter (CDOM), as in dissolved nutrients, and as such limits light availability for photosynthesis (Nittrouer et al., 1986; Smith and DeMaster, 1996; Del Vecchio and Subramaniam, 2004). In the low SSS and mesohaline SSS with higher particulate loads, DDAs reside below the plume, or penetrate deeper depths favored by sufficient light for photosynthesis. Indirectly, DIN and Si negatively impacted het-2 through turbidity (predictor), which is consistent with the findings of Goes et al. (2014) that the phytoplankton community structure during the 2010 cruise was in a Si deficiency and largely controlled by salinity and nutrient limitation. Hence, much of the environmental influence on the DDAs seems to be due to the indirect effects of turbidity, or the AR plume, mediating light quality and nutrient concentrations.

One of the more surprising outcomes of the SEM model was the robust positive influence (0.79, p < 0.0001; R²= 0.71) on het-2 had on het-1 abundance. Het-1 and het-2 co-varied with depth; similar distributions and strong correlations between the two symbiotic strain abundances have been previously documented in diverse habitats, including the Western Tropical South Pacific (WTSP), South China Sea, North Pacific Gyre, and earlier in the WTNA (Stenegren et al., unpublished). However, in correlation analyses, the direction of influence cannot be resolved, as in our SEM model. The specific mechanism of het-2 influence on het-1 abundance remains unknown.

Still unresolved is why the Hemiaulus-Richelia symbioses tend to dominate in this region of the ocean rather than the Rhizosolenia-Richelia symbioses. The host taxonomy and morphological characters differ greatly, e.g., cell diameter (Sournia, 1986; Villareal, 1992), which equates to different Si requirements. The difference in host frustule dimension is one obvious precondition that could drive differences in parameter influence. Si was not directly influencing het-1 or het-2 abundance, however, Si was considered an indirect negative predictor of het-2 abundance in 2010. The latter supports the conclusion by Goes et al. (2014) of an efficient Si drawdown mediated by both non-symbiotic and DDA populations leading to an overall Si deficiency. Moreover, the SEM model shows direct positive effect of salinity on het-1 and an indirect negative effect of salinity on het-2. Therefore, het-1 seems to thrive more at higher salinities as detected by our qPCR, and explains the higher abundances of het-1 over het-2 in the open ocean habitats (Church et al., 2008; Stenegren et al., unpublished).

In conclusion, our study represents the first application of SEM to determine the direction of environmental parameter impact on DDA abundances, which is a departure from the commonly used univariate statistics (e.g., correlation based analyses). The resemblance in our SEM model predictions of the two cruises (bloom vs. background densities) demonstrated that several environmental parameters (temperature, salinity, turbidity, and dissolved nutrients) acted directly and in concert to influence DDA abundances in a highly dynamic ecosystem. The SEM model output also establishes a framework of new and testable hypotheses by lab based experiments and/or similar application of the SEM to other datasets.