New molecular markers for quantifying abundance and gene expression of widespread and often undetected marine N 2 fixing symbionts of diatoms

Decades of research has been devoted to understanding the occurrence, distribution, and activity of different N 2 -fixers in the oceans. Marine N 2 -fixing diatom symbioses involving the heterocyst-forming cyanobacterium Richelia are widespread, yet often go undetected. Some of the diatom- Richelia symbioses have the advantage to be identified by microscopy, however, observations alone cannot inform on the metabolic state of a cell. Here, we developed nine new specific quantitative PCR assays that detect three common symbiotic Richelia strains (ReuHH01, RintRC01, RrhiSC01) of diatoms based on molecular markers for: a high affinity phosphate transporter (pstS), an iron transporter (exbB), and a constitutively expressed gene (rnpA, protein component of ribonuclease P). We first tested the new assays in the lab, including a diel experiment to elucidate the temporal dynamics of gene expression for each molecular marker, then applied the new assays to field collected samples. In the field we also made microscopy observations of the diatom- Richelia symbioses and measured bulk N 2 and Carbon (C) fixation rates with stable isotope labelling incubations. The number of genes encoding the molecular markers varied; typically, fewer in the endobiont strains. Detection in field samples based on qPCR was consistent with microscopy observations. The environmental gene expression for all the new targets, and additionally nifH for N 2 fixation (nitrogenase) were low, and highest expression was detected in the upper water column (0-40 m) consistent with higher densities of the Richelia by qPCR and microscopy observations. Low in situ bulk rates of N 2 and C fixation corroborated the aforementioned low nifH expression. Expression of both nifH and pstS was temporally regulated in the lab experiment of the facultative symbiotic Richelia strain RrhiSC01 with higher expression in the early and late photoperiods, respectively. The nine new assays are an improvement over the previous assays based on the nifH gene as cross-reactivity between Richelia strains is minimum, expression can be normalized to rnpA, and expression is informative of the symbiont nutrient status.

In the open oligotrophic ocean, a few genera of diatoms form stable partnerships, or symbioses, with several strains of the filamentous heterocyst-forming cyanobacteria Richelia. Heterocysts are specialized cells for N 2 fixation, or the reduction of di-nitrogen (N 2) to biologically available nitrogen (N). Two diatom- Richelia symbioses, Hemiaulus-Richelia and Rhizosolenia-Richelia, form large and expansive blooms, which contribute significantly to both new and primary production (Mague et al., 1974, 1977; Venrick, 1974; Subramaniam et al., 2008; Karl et al., 2012). At such high abundances, typically only one diatom- Richelia symbiosis dominates, however, at background densities often the different diatom- Richelia symbioses co-occur (Karlusich et al., 2021). The environmental niches that are occupied by different diatom- Richelia symbioses that allow them to form high density blooms are unknown, as are the ecological and biogeochemical consequences of such changes in dominance. Although these populations are important to study, the collection, detection, and identification of the symbioses are challenged by the host cell morphology (e.g., formation of chains), patchy distributions, low densities in comparison to sympatric populations (e.g., picocyanobacteria), and require a skilled taxonomist. Given the importance of the diatom- Richelia symbioses to the biogeochemical cycles of N and carbon (C), identifying an efficient means to collect, identify and distinguish the different symbioses and simultaneously study their metabolic potential is valuable.

The diatom- Richelia symbioses are a unique system as the symbiont cellular location varies from internal to external, and includes an intermediary, where the symbiont resides in the periplasmic space of the host diatom (between the silicate frustule and plasma membrane) (Norris, 1961; Villareal, 1992; Janson et al., 1995; Caputo et al., 2019; Flores et al., 2022). Thus, the immediate chemical environment varies for each Richelia strain. Recently, the Richelia spp. were renamed and re-classified based on a genome phylogeny (Foster et al., 2022): R. euintracellularis (ReuHH01), R. intracellularis (RintRC01), R. rhizosoleniae (RrhiSC01). Additionally, eight environmental metagenomic assembled genomes (MAGs) have been reported that are highly similar to Richelia (Satinsky et al., 2014; Parks et al., 2017; Poff et al., 2021; Delmont et al., 2021; Paoli et al., 2022; Foster et al., 2022). Only one strain, the facultative external symbiont R. rhizosoleniae RrhiSC01 has been maintained in the lab for extended periods, although without its host diatom (Foster et al., 2010).

Traditionally, the diatom- Richelia symbioses have been identified and enumerated in field samples by microscopy (Heinbokel, 1986; Villareal, 1992). Later, cell-based counting was replaced by quantitative polymerase chain reaction (qPCR) assays based on the symbiont’s nifH gene for nitrogenase, the enzyme for N 2 fixation (Church et al., 2005; Foster et al., 2007). qPCR is an attractive alternative to cell counting as several different gene targets can be assessed from the same template pool. However, in the case of the nifH qPCR assays for detection of the Richelia symbionts, there is a known cross-reactivity between assays for two (ReuHH01, RintRC01) of the three symbiont strains (Foster et al., 2007; Stenegren et al., 2018). Hence the separate populations risk being overestimated. Moreover, currently the results of RNA based qPCR studies for diazotrophs are seldom normalized to a constitutively expressed gene. Thus, developing new assays which can distinguish the various symbiont strains and also provide normalization are necessary and relevant for appropriate quantification.

Here we developed nine new qPCR assays which can detect the three common Richelia strains, including assays for a constitutively expressed gene (protein component of ribonuclease P, rnpA) for normalization. The two other gene targets were a high-affinity ABC transporter phosphate-binding gene (pstS) and a gene involved in iron (Fe) transport (exbB). The latter were selected as both were found in all three Richelia strains, and when expression is detected, these markers have the potential as bioindicators of nutrient status (pstS for P limitation, exbB for Fe demand) for the symbionts (Harke et al., 2012; Dyhrman and Haley, 2005; Pereira et al., 2016, 2019; Held et al., 2020). Optimized assays were applied to both lab and field collected samples along a transect in the North and South Atlantic where conditions of low dissolved inorganic nitrogen (DIN) and co-limitation of Fe and phosphorus (P) were expected based on previous studies (Mills et al., 2004; see references in Benavides and Voss, 2015). Bulk incubation experiments were also performed along the transect to measure N 2 and C fixation rates. Finally, a lab based diel experiment analyzed the temporal regulation of the target genes in Richelia strain RrhiSC01, the only symbiont strain currently available in culture.

Field expedition

Samples were collected on board R/V Atlantis (8Feb-11March 2018) at 20 stations in the North and South Atlantic on the TriCoLim research cruise AT39-05 (Figure 1; Supplementary Table 1) using a conductivity temperature depth (CTD) rosette fitted with 12 L niskin bottles.

Nutrient and chl a analyses

Samples (5–15 mLs) for nutrient analyses were taken from 4–6 depths in the upper 100 m of the water column and immediately frozen until laboratory analysis. Dissolved inorganic phosphate (DIP) and ammonium (NH 4 +) were measured as previously described (Noble et al., 2012) using a Technicon AutoAnalyzer IITM components. Dissolved silicate (DSi) and total inorganic nitrogen (TIN: nitrate + nitrite) were measured on an Alpkem rapid flow analyzer (RFA) 300TM. The limits of detection (LOD) were 0.011 μM for DIP, 0.063 μM for TIN, 0.186 μM for DSi, and 0.133 μM for NH 4 +.

Samples (1 L) for chl a determination were collected from 5 depths (5–90 m) and directly filtered onto a 25 mm diameter glass fiber filter (GFF) using vacuum filtration. Chl a was extracted overnight in 1.4 mL of 90% acetone at -20°C and measured on a Turner Trilogy model with a non-acidified chlorophyll module using previously published equations (Welschmeyer, 1994). The LOD was 459.66 relative fluorescence units (RFU) which corresponds to an approximate concentration of 0.03 μg Chl a /L.

Shipboard microscopy observations

At each station, a plankton net (75 μm mesh size) fitted with a flowmeter was deployed for 20 min at the surface by hand. One mL of the sample was enumerated using a Sedgewick rafter counting chamber. Epi-fluorescent microscopy was used to identify the various diatom- Richelia symbioses by blue (450–490 nm) and green (510–560 nm) excitation. The stations selected for setting up the bulk water stable isotope incubation experiments (see below) were those that had cell densities of at least 50 symbiotic cells L -1. Later these observations were used to select stations for estimating gene abundances and gene expression.

Incubation experiments

The incubation experiments were amended with stable isotopes (15 N 2 and 13 C-bicarbonate) and set up with bulk seawater collected from the CTD or a trace metal water sampler from 5 or 15 m, respectively (Figure 1; Supplementary Table 1). For each experiment, transparent polycarbonate bottles (2.75 L) fitted with septa caps were filled without air bubbles: 3 replicates each for time zero (T 0), control, and for bottles amended with DOP. For the DOP amended bottles, 1 mL of 20 μM DOP solution (1:4 mixture of 2-aminoethylphosphonic acid, a phosphonate and beta-glycerophosphate disodium salt hydrate, a phosphomonoester; Sigma Aldrich, Stockholm, Sweden) was used. These represent two common DOP compounds in the ocean and the ratio is based on the proportion of phosphonates (25%) and phosphomonoesters (75%) in the total DOP pool (Hong et al., 1995; Clark et al., 1999; Kolowith et al., 2001). Bottles were amended with 1.5 mL of 15 N 2 gas (98% 15 N 2) and 0.5 mL of a 500 μM 13 C-bicarbonate solution (NaH 13 CO 3 -) (Cambridge Isotope Laboratories, MA USA) using gas tight syringes. Bottles were kept in an on deck transparent plastic chamber with continuously flowing surface seawater and appropriate shading (25% shading). Triplicate bottles were filtered at each time point (time 0, 24 hrs) onto pre-combusted 25 mm diameter GFF filters. The GFF samples were dried at 65°C and acidified overnight in a desiccator containing a petri dish with 37% hydrochloric acid (HCl) to aspirate un-incorporated 13 C. Subsequently, the filters were folded and prepared into tin capsules for analyzing 15 N and 13 C values using a Delta V Advantage Mass spectrometer coupled to a Costech Instrument Elemental Combustion system at the University of Hawaii Biogeochemical Stable Isotope Facility. Labeling percentages and N 2 and C fixation rates were calculated as described previously (Montoya et al., 1996).

Nucleic acid sample collection and extractions

Seawater samples (2.5L) for nucleic acids were collected from 4–5 depths in the upper water column (0–95 m) into pre-washed (10% bleach) polycarbonate bottles. Samples were immediately filtered onto 0.2 μm pore size Supor filters (Pall Corporation) held in a 25 mm diameter filter holder using a peristaltic pump. Samples for RNA contained 350 μl RLT buffer amended (Qiagen RNeasy mini plant kit; Qiagen, Hilden, Germany) with 10% beta-mercaptoethanol; all samples were flash-frozen in liquid N 2 immediately and stored at -80°C. DNA samples collected at stations 1, 6-7, 11 and 16–17 were extracted using the modified DNAeasy Plant kit (Qiagen) method (Foster et al., 2010) and the final elution volume was 50 μl. Total RNA from the environmental samples and the RrhiSC01 diel experiment (see below) was extracted using the Qiagen RNeasy mini plant kit (Qiagen) following the kit instructions and previously described modifications, including a 1 hr DNAse step (Foster et al., 2010). The field collected RNA samples were re-eluted in 20 μl and the RNA extracts from the diel experiment were re-eluted in 45 μl. All RNA extracts were stored at -80°C.

Lab experiment for assessing temporal dynamics of gene expression in Richelia RrhiSC01

Triplicate (800 mL) RrhiSC01 cultures were maintained in SO media made with a seawater base from the Sargasso Sea (Waterbury et al., 1986) at 26°C and illumination of 30 µmol photons m -2 s -1 and a 14:10 hour light: dark cycle. The initial cell density was determined by Chl a concentration as 2.46 ± 0.84 µg/mL. Samples (25 mls) were taken every 2–4 hours for 48 hours (Supplementary Table 2) from each replicate RrhiSC01 culture by gentle vacuum onto a 1 μm pore size filter (25 mm diameter). Filters were placed in pre-sterilized 2 mL tubes containing zirconium dioxide ceramic beads (Precellys ®, VWR) and 350 µl RLT buffer (Qiagen) amended with β-mercaptoethanol (1%), flash frozen in liquid nitrogen and kept in - 80°C until extraction.

Oligonucleotide design

Oligonucleotides (Supplementary Table 3) that target exbB, pstS and rnpA for the three Richelia strains: ReuHH01, RintRC01, RrhiSC01 were designed using Applied Biosystem’s (AB) Primer Express software. The respective sequences were derived from the available draft genomes (Hilton et al., 2013; Hilton, 2014) (Supplementary Table 4) and are as follows for ReuHH01, RintRC01, and RrhiSC01 for exbB: IPCPKDNE_00723, JHFCPCOE_01496, IACDABIC_00827; for pstS: IPCPKDNE_02124, JHFCPCOE_04511, IACDABIC_03299; for rnpA: IPCPKDNE_02195, JHFCPCOE_01428, IACDABIC_04460. Probes were either 5’ labeled with a fluorescent reporter FAM (6-carboxyfluoresceom) or VIC, and all probes were 3´ labeled with TAMRA (6-carboxytetramethylrhodoamine) as a quenching dye. The nifH oligonucleotides were as previously described (Church et al., 2005; Foster et al., 2007). Lab based cross-reactivity tests were performed for each newly designed oligonucleotide and TaqMAN ® assay.

Quantitative polymerase chain reaction assays

Total RNA (2-3.5 μl or 3.5-90.7 ng/μl) was reverse transcribed (RT) using a commercially available kit (SuperScript ® III First-Strand Synthesis System) as previously described (Foster et al., 2010). Triplicate 20 μl RT reactions were performed and pooled prior to running the qPCR assays. Reactions without reverse transcriptase (No RT) and negative controls (water as template) were also run. cDNA was stored at -20°C until used in the qPCR assays. For the diel experiment with RrhiSC01 isolate, 2 μl of total RNA (12.5–170 ng/μl) template was used in the triplicate RT reactions, which were pooled prior to running the qPCR for the various target genes.

Gene abundances (DNA, cDNA) were quantified by TaqMAN qPCR assays on a StepOnePlus system (Applied Biosystems) using default conditions for the reaction parameters (50°C for 2 min, 95°C for 10 min, 45 cycles at 95°C for 15 s and 60°C for 2 min. The reaction volume was 25 μl (12.5 µl 2X buffer, 0.5 µl probe (10 μM), 1 µl each of each primer (10 μM), 8 µl nuclease-free water for each reaction and including 2 µl of template DNA (or cDNA). Standard curves of synthesized gBlocks (IDT, Leuven, Belgium) (exbB, pstS, and rnpA) or a plasmid containing the target (nifH) were made in dilution series ranging 10 8 to 10 1 gene copies per reaction and used to estimate the gene copies in the samples. Samples were run in triplicate, and additionally four no template controls (NTCs) and the No RT controls were run. Inhibition controls with 2 µl of 10 5 standard and 1 µl of sample as template were run on a few initial samples to verify that the reactions were not inhibited. No inhibition was detected.

Statistical analysis

Statistical analyses were done using R (version 4.3.2). The periodicity and identification of peak gene expression by RrhiSC01 for exbB, nifH and pstS normalized to rnpA expression in the diel experiment was tested using the R package RAIN (Rhythmicity Analysis Incorporating Non-parametric Methods; Thaben and Westermark, 2014). The effect of station and DOP amendment on the bulk N 2 and C fixation rates were tested using the non-parametric Kruskal-Wallis test. The gene copy abundance of nifH and rnpA from the environmental samples was compared with a non-parametric Wilcoxon test for the respective strains ReuHH01 and RintRC01. A non-parametric Spearman’s rank correlation was conducted to test possible correlations between expression of exbB, nifH, and pstS by the Richelia strains (ReuHH01, RintRC01) and dissolved nutrients and depth of sampling. The resulting correlation matrices were visualized in a heat map of hierarchical clustering in R (ver. 3.2.2) using packages ‘hmisc’ and ‘gplots’. For samples where expression was below detection (bd), expression was set to 0.

BlastN analyses were performed to identify the cross-reactivity for the oligonucleotides. The cross-reactivity tests used a database of all cyanobacterial genomes in GTDB (version R08-RS214), the four Richelia draft genomes, and eight MAGs (MAG 1-8) (Foster et al., 2022), and a threshold e-value < 0.1. BlastN was used to identify the number of gene copies (exbB, nifH, pstS, rnpA) and estimate the sequence similarity between copies (in genomes that had multiple copies). A pairwise comparison was performed between the multiple copy genes in the draft genomes (ReuHH01, RintRC01, RrhiSC01; Hilton et al., 2013; Hilton, 2014) using the EMBOSS program (v.6.6.0.0; Rice et al., 2000). To assess the phylogenetic relationship between the multiple copies of a particular gene, the multi-copy gene sequences were aligned and neighbor-joining phylogenies were constructed using Molecular Evolutionary Genetic Analysis version 11 (MEGA11). The phylogenies were visualized in Tree Viewer (Tamura et al., 2021; Bianchini and Sánchez-Baracaldo, 2024).

Gene copies/genome

The selected molecular markers were identified in the various Richelia draft genomes and eight environmental MAGs, with a few exceptions (Table 1). In general, the highest number of gene copies for each target gene were identified in the facultative symbiont RrhiSC01, followed by the endobionts RintRC01 and ReuHH01, and the respective environmental MAGs (Table 1). Those genomes which lacked the molecular markers had poor completeness scores by checkM2 (e.g. <90%).

The multiple copies of the target genes (exbB, nifH, pstS) present in the draft Richelia genomes (ReuHH01, RintRC01 and RrhiSC01) were considered highly divergent based on a pairwise comparison (Supplementary Table 5). For example, the two copies of nifH in ReuHH01, RintRC01, and RrhiSC01 had sequence identities of 42.0, 43.7, and 44.0%, respectively. Similarly, the two copies of pstS in RintRC01 had a sequence identity of 42.3%; while the four copies of pstS in RrhiSC01 had sequence identities 41-52.6%. The pairwise comparison of the three copies of exbB in RrhiSC01 were also dissimilar (e.g., 41.1-66.4% sequence identity). Additionally, we ran a BLASTN of the newly designed oligonucleotides against the multiple copies to test for possible cross-reactivity. We found little to no sequence similarity between the newly designed oligonucleotides and the multiple copies e.g., 42 -75% sequence similarity across the expected amplicon region (Supplementary Table 6). Moreover, phylogenies based on the target genes present in more than one copy per genome showed that the single copies of pstS and exbB targeted by our assays are orthologs, originating from the same ancestral gene (Supplementary Figure 1). As for the nifH genes, the single copies targeted in ReuHH01 and RrhiSC01 are also orthologs but the nifH copy targeted in RintRC01 lays on a separate branch and is therefore of paralogous origin (Supplementary Figure 1).

Oligonucleotide design

In order to check the specificity of the newly designed qPCR assays, the oligonucleotide sequences were tested by BLASTN analyses and cross-reactivity assays. Targets that had significant alignment for each oligonucleotide were identified and corresponded to a sequence similarity of 92.6-100% across the entire fragment of the expected amplicon (Figure 2; Supplementary Table 6). The newly designed oligonucleotides were 100% identical to their respective Richelia strain and to the environmental MAGs that cluster with each strain (Figure 2). The only exceptions were the two largely incomplete (73.3 and 89.3%, respectively) genomes: R. euintracellularis (ReuHM01; GCF_00350125.1) and one of the R. euintracellularis MAGs (MAG6: MARD_SAMEA2272302_REFG_MMP2272302), which contained 1 mismatch in the probe and resulted in a lower sequence similarity for exbB (98.8%). Additionally, both rnpA and pstS were absent in these genomes (Table 1), likely explained by both genomes having a low genome completeness. The MAG Candidatus Richelia exalis (MAG 8) had a high sequence similarity to the oligonucleotides targeting ReuHH01. We identified that the previously designed oligonucleotides targeting ReuHH01 nifH had only 98.4% sequence similarity to all R. euintracellularis strains (including ReuHM01) due to a mismatch on the third position from the 5’ end of the probe.

Apart from Candidatus Richelia exalis (MAG 8), there is no expected cross-reactivity with other organisms (e.g., environmental strains) or between the Richelia strains for the rnpA, exbB and pstS assays. Three freshwater cyanobacterial strains of Trichormus variabilis had a 95.2-96.8% sequence similarity to the nifH oligonucleotides targeting RrhiSC01 and one unclassified Nostocales cyanobacterium had a 95.3% sequence similarity to the nifH assays targeting RintRC01 (Figure 2). Lab based cross-reactivity assays indicated that the new assays for rnpA, pstS, and exbB had little to no cross-reactivity between the strains (Supplementary Figure 2).

Temporal regulation of gene expression in Richelia RrhiSC01

To understand the diel regulation in gene expression, we studied gene expression for exbB, pstS, and nifH in the facultative symbiont RrhiSC01 grown in N deplete media over two consecutive light and dark periods (Figure 3; Supplementary Table 2). The expression of rnpA (cDNA gene copies L -1) was used for normalizing. RAIN analysis identified a significant periodicity in pstS and nifH gene expression (p-value = 6.80e -10 and p-value = 1.17e -07, respectively). The peak expression of nifH occurred in the early photoperiod, or after 4 hours of light (L4); maximum expression of pstS was detected in the late photoperiod, or after 13 hours of light (L13). Although not significant, exbB expression peaked similar to nifH early in the photoperiod (L7). In general, all genes had lower expression levels during the dark period (Figure 3; Supplementary Table 2).

Hydrographic conditions

The sea surface temperature (SSTs) ranged 21.1-28.0 during the transect in the N. Atlantic; warmest (>27.5) at stations 6–14 and cooler (<24.1) at stations 1-3 (Supplementary Figure 3). Depth profiles for TIN, DIP, and dissolved silicate (DSi) followed expected nutrient profiles with depletion in the surface and higher concentrations at depth (e.g., below 50 m) (Supplementary Figure 5; Supplementary Table 7) with a few exceptions. For example, in the near surface samples (20 m) elevated DSi and DIP were measured at stations 6–8 and 16–17 for DSi and at stations 6–8 for DIP.

Shipboard microscopy observations

In general, observations from the net tows showed that the two symbiotic diatoms, Hemiaulus and Rhizosolenia, co-occurred at similarly low densities (10–20 symbiotic cells L -1) (Figure 1). The highest densities of both symbioses were observed at station 17 (> 100 symbiotic cells L -1) and fewer of both symbioses at stations 8, 16, 18-19, while densities of Rhizosolenia-Richelia symbioses were elevated at stations 1-3, 7, 11-12 (Figure 1).

Gene abundances of wild Richelia spp.

Based on the higher densities of the Rhizosolenia-Richelia and Hemiaulus-Richelia symbioses by microscopy observations, stations 1, 6-7, 11 and 16–17 were chosen to compare abundances of the endobionts RintRC01 and ReuHH01 based on nifH and rnpA genes. The Chaetoceros-Richelia symbioses were never observed, and therefore the facultative symbiont RrhiSC01 were not assayed in the field samples. In general, detection of RintRC01 was higher than ReuHH01 in the depth profiles. Gene copy abundances of nifH were significantly higher than that of rnpA for both endobionts (Wilcoxon test; p-value = 0.016 and p-value = 0.02, respectively) (Figure 4; Supplementary Table 8). On average, the nifH assay resulted 5.7- and 6.5-times higher gene copy abundances than the rnpA assay for ReuHH01 and RintRC01, respectively. The nifH abundances were also 2–3 orders of magnitude higher than the densities of ReuHH01 and RintRC01 by microscopy observations, while abundances based on rnpA copies were typically similar or 1–2 orders higher than the microscopy counts (Figure 1b; Supplementary Table 8). Moreover, the gene copy abundances of RintRC01 were higher in the surface water (5–15 m) except at two stations (1 and 16) where above 40 m there was no detection. The depth profiles of ReuHH01 had a higher variability compared to profiles of RintRC01 and gene copy abundances were higher in the samples from depth. For example, at several stations (1, 7, 16, and 17), the highest abundances were from the samples collected between 40–60 m (Figure 4; Supplementary Table 8).

Environmental gene expression of exbB, nifH, pstS for Richelia strains

In general, expression levels of exbB, rnpA and pstS were low and often at the detection limit (1–10 copies) for both ReuHH01 and RintRC01 (Supplementary Figure 5; Supplementary Table 9);. Expression of rnpA by ReuHH01 was detected in all three technical replicates in 53% of samples (18 out of 34), compared to 32% of samples (11 out of 34) for RintRC01. In samples that estimated expression of rnpA of ReuHH01 and RintRC01, the expression (cDNA gene copies L -1) ranged between 1.20 x 10 3 -9.13 x 10 3 and 9.35 x 10 2 -1.02 x 10 3, respectively.

For the other molecular markers, nifH expression by ReuHH01and RintRC01 was in general higher and detected in more samples than exbB and pstS (Supplementary Table 9) Moreover, at deeper depths (>40 meters), expression of nifH by RintRC01 was detected but expression of rnpA and/or the other target genes was not detected (Supplementary Table 9). For ReuHH01, pstS showed a higher and more consistent expression in comparison to exbB transcripts, which was only detected in a few samples and typically in the upper water column (e.g., 5-35m) (Supplementary Figure 5; Supplementary Table 9). There was one exception at station 1, where exbB expression was detected at 70 m. Pearson correlation analysis of the normalized gene expression (exbB, nifH and pstS) by ReuHH01and RintRC01 with environmental variables showed that there was a significant and negative correlation between pstS expression and depth (Pearson correlation coefficient: -0.46, -0.48, p-value: 0.036, 0.046), for ReuHH01and RintRC01, respectively (Figure 5). The nifH and exbB expression by RintRC01 was also significantly negatively correlated with depth (p-value = 0.0004 and p-value = 0.0148, respectively; Figure 5). Moreover, the normalized pstS and exbB expression were positively correlated for both ReuHH01and RintRC01. (Pearson correlation coefficient: 0.70, 0.54, p-value: 4.00 x 10– 4 and 0.02, for ReuHH01and RintRC01, respectively) (Figure 5). Several nutrients were positively correlated with depth; expression of pstS by ReuHH01 was negatively correlated to nitrate; expression of all three gene targets by RintRC01 was significantly negatively correlated with phosphate while only nifH and pstS by RintRC01 was significantly negatively correlated with TIN (Figure 5). In general, expression was very low by both RintRC01 and ReuHH01, and as such few data were available for statistical analyses (Supplementary Table 9).

Bulk N 2 and C fixation rates

Bulk N 2 and C fixation rates were assayed at select stations (stations 1, 5, 6, 11, 15, 17) and at five of these stations, experiments included DOP amendments as a treatment (see methods) (Figure 6; Supplementary Table 10). Consistently, N 2 fixation rates were low (0.01-0.024 nmol N mL -1 h -1), but comparable with previous measurements in the North and South Atlantic (Figure 6; Supplementary Table 10; Supplementary Table 11). Highest N 2 fixation rates (0.24 ± 0.05 nmol N mL -1 h -1) were measured at station 17, which were also coincident with the highest densities of the symbiotic diatoms (Figures 1, 6; Supplementary Table 10). Carbon fixation was highly variable and in general appeared to decrease along the transect (8.85 ± 2.0 nmol C mL -1 h -1 to 0.24 ± 0.15 nmol C mL -1 h -1), which was consistent with the lower chl a concentrations measured in the upper water column at stations 11, 15 and 17 (Supplementary Figure 4; Supplementary Table 10). There was a significant difference in both N 2 and C fixation rates between the stations (Kruskal-Wallis test: p-value = 1.65 x 10– 3 and p-value = 2.27 x 10 -4, respectively). The amendment of DOP resulted in a significant decrease in N 2 fixation in all experiments, however no significant impact on the measured C fixation rates (Supplementary Table 10; Kruskal-Wallis test p-value = 0.042 and p-value = 0.17, respectively).

The diatom- Richelia symbioses are widespread and typically highly active N 2 -fixing populations in the global oceans, and therefore important for primary and new production (Carpenter et al., 1999; Subramaniam et al., 2008; Karl et al., 2012; Luo et al., 2012; Karlusich et al., 2021; Shao et al., 2023). Despite their pivotal roles, their presence and activity are often undetected due to difficulties in collection, low background densities and patchy distributions. Additionally, the diatom- Richelia symbioses can co-occur with other diazotrophs making it difficult to distinguish their activities and contributions to new and primary production. Currently it is unknown which conditions favor one N 2 -fixer over another, especially why one symbiosis tends to bloom over another. This work aimed at providing new qPCR methods to easily detect and distinguish between the three Richelia symbionts and when combined with gene expression assays and environmental monitoring (e.g. dissolved nutrients) can potentially provide additional information on nutrient status of the symbionts, e.g., P limitation and demand for Fe. The diatom hosts are highly dependent on their respective symbionts for growth (Foster et al., 2011, 2021), thus understanding how the symbionts are limited is important for the overall state of the symbioses.

Newly designed assays are specific

A common method of enumerating the relative abundance of a particular diazotroph, including Richelia, is to use its nifH gene in qPCR assays. These assays are typically single target and use DNA or RNA templates for detecting abundance and expression, respectively. However, few environmental RT-qPCR studies for diazotrophs have addressed the known cross-reactivity between the two closely related endobiont Richelia strains (ReuHH01, RintRC01) (Foster et al., 2007; Stenegren et al., 2018). This can be problematic when the two strains co-occur at densities equivalent to 10 3 copies L -1, hence one aim of this work was to develop new qPCR assays that could distinguish between the closely related strains of Richelia. Additionally, we chose gene targets that would be indicative of nutrient status for P and Fe, as the latter are important elements for phytoplankton, and in particular diazotrophs (Orchard et al., 2009; Schoffman et al., 2016; Hutchins and Sañudo-Wilhelmy, 2022).

The Richelia spp. genomes vary in size and content which is also related to the symbiont cellular location (Hilton et al., 2013); the endobionts possess smaller, more streamlined genomes that differ in transporter content (Nieves-Morión et al., 2020). Thus, we initially selected molecular markers (pstS, exbB, rnnpA) for the qPCR assays that were present in all three Richelia strains. Both the lab based cross-reactivity tests and the in silico (BLASTN) analyses showed a minimal cross-reaction among the newly designed oligonucleotides. In the lab cross-reactivity tests, amplification was non-linear and poorly detected. Moreover, cross-reaction was only found at the higher standard concentrations (10 7 and 10 8 gene copies) which correspond to co-occurring densities that have not been reported in nature. The new assays are a great improvement over the nifH based assays which have been shown to cross-react and overestimate the abundances of Richelia in environmental samples (Foster et al., 2007; Stenegren et al., 2018). In fact, direct evidence of the cross-reaction was shown in the environmental samples. At stations 1, 11, and 17, gene copy estimates based on nifH were approximately one order of magnitude higher than rnpA, and rnpA copy numbers were about 10 3 copies L -1 which is the approximate abundance that is known to cause the assays to cross-react (Foster et al., 2007; Stenegren et al., 2018). An overestimation of nifH can also be expected in the environmental gene expression. For example, expression of nifH was detected for ReuHH01 and RintRC01 (2.15 x 10 4 – 2.76 x 10 4 and 4.4 x 10 3— 8.6 x 10 3 cDNA copies L -1, respectively) while the expression of rnpA by RintRC01 was either undetected or below detection limit. Likewise, the detection (DNA) and expression level (cDNA) of nifH was consistently higher than that of exbB and pstS. Even though the discrepancy is due to the cross-reaction of the nifH assays, it can in part also be attributed to a general higher expression of nifH than the other target genes. For example, in the lab based diel experiment with RrhiSC01, nifH expression was higher than that of exbB, rnpA and pstS at all measured time points.

The cross-reactivity analyses also highlighted that a recently reported MAG, Candidatus Richelia exalis (Delmont et al., 2021), had potential cross-reaction in the assays that were designed to target ReuHH01. Candidatus Richelia exalis forms a sister clade to R. euintracellularis ReuHH01 (Foster et al., 2022), and is thus closely related. Candidatus Richelia exalis was also obtained from metagenomic data from the larger size fraction (20-180 μm) (Delmont et al., 2021) and has a modest average nucleotide identity (ANI) (87.3%) and similar genome content to ReuHH01 (Delmont et al., 2021). Thus it is likely that Candidatus Richelia exalis has a symbiotic lifestyle. Unexpectedly, this MAG was reported from a region off South America within the equatorial Pacific Ocean, an area not well-sampled for diazotrophs, including the diatom- Richelia symbioses. The in-silico specificity analyses also identified a 95-97% sequence similarity between the nifH oligonucleotides that target RintRC01 and RrhiSC01 and a few other strains of Nostocales. However, the isolation sources according to the database (GTDB) for the other Nostocales were from terrestrial and freshwater environments, and hence not likely to co-occur in marine waters where Richelia strains exist. In summary, the new assays for rnpA, exbB, pstS are specific for the three Richelia strains and are thus an improvement over the frequently used nifH assays.

Temporal dynamics in nifH, pstS, and exbB gene expression in lab and wild Richelia populations

The application of the new assays requires an understanding of how the gene expression for each target is regulated with regards to in situ conditions, e.g., time of day. This is particularly important when working with wild populations. RAIN analyses identified that both nifH and pstS were significantly temporally regulated in the facultative symbiont RrhiSC01. RAIN is a nonparametric method that can be used to detect biological rhythms and requires a routine sampling strategy (Thaben and Westermark, 2014). For example, a RAIN analysis was applied to a transcriptomic study of the Rhizosolenia-Richelia symbioses in the North Pacific, where samples were taken ever 4 hours over three consecutive days in order to identify if certain gene transcription by both the host diatom and cyanobacterial symbiont were temporally regulated (Harke et al., 2019). Although, our field sampling in the North Atlantic also included samples collected at various times of day, the detection was limited and samples were collected from multiple depths and stations without consistency. Therefore RAIN analysis was not possible. However, we still take advantage here of our environmental gene expression dataset for the other two Richelia symbionts (ReuHH01, RintRC01) and compare with other similar works on how environmental conditions might affect gene expression.

The nifH gene encodes the nitrogenase iron protein and expression of nifH in some cyanobacteria exhibits diel regulation (Postgate, 1998). Nitrogenase is also sensitive to and becomes deactivated in the presence of oxygen. Heterocyst-forming cyanobacteria can perform nitrogen fixation under aerobic conditions by differentiating heterocysts to house and protect nitrogenase from oxygen (Haselkorn, 1978). Our results presented here on RrhiSC01 agree with an earlier lab-based study that nifH expression is highest during the early to mid-part of the photoperiod (Foster et al., 2010). Similarly, field-based studies have shown highest nifH expression in the early light period for the other Richelia symbionts: RintRC01 in the North Pacific (Harke et al., 2019; Church et al., 2005b) and midday for both ReuHH01 and RintRC01 in surface waters surrounding Cape Verde (Turk-Kubo et al., 2012), an area in close proximity to stations 1–3 of our study. Thus, it seems the nifH expression of Richelia symbionts, including the facultative RrhiSC01, is consistent with active N 2 fixation during the day, which is expected for heterocyst-forming cyanobacteria as energy derived from C fixation drives N 2 fixation (López-Igual et al., 2010). High expression during the day is also consistent with previous work showing high daytime activities for N 2 fixation during diatom- Richelia blooms (Venrick, 1974; Mague et al., 1974, 1977; Carpenter et al., 1999; Subramaniam et al., 2008; Shao et al., 2023).

Both photosynthetic activity and growth require a high amount of P. The pstS gene encodes for a high-affinity ABC Phosphate transporter, and often pstS expression is used as an indicator for P limitation in bacteria, including N 2 -fixing microorganisms (Harke et al., 2012; Pereira et al., 2016, 2019). Previous studies on diazotrophs have shown that pstS has a diel transcription pattern. For example, in the unicellular cyanobacterium C. watsonii WH 8501, and additionally in wild populations of RintRC01, pstS expression is highest approximately 10 hours into the photoperiod (Shi et al., 2010; Pereira et al., 2016; Harke et al., 2019). Thus, the high pstS expression by RrihSC01 in the late light period is consistent with these other diazotrophs and the expected higher P demand from photosynthesis (Orchard et al., 2009; Lin et al., 2016). Note that in the lab-based study presented here, the RrhiSC01 is growing in a free-state without its symbiotic diatom which could influence the gene transcription patterns. DNA replication also contributes to a higher P requirement (Zhao et al., 2022). A temporal pattern was only true for one of the two pstS gene copies found in the transcriptomic study of wild RintRC01 symbiont in the North Pacific, where the second copy did not have a periodic expression (Harke et al., 2019). The pstS gene copy that did not have a periodic expression is the copy targeted by the assay presented here (gene id: JHFCPCOE_04511). The divergent temporal expression patterns could reflect a difference in function of these two gene copies, which is common for paralogous genes (Tatusov et al., 1997).

In aquatic ecosystems, Fe (as Fe II or Fe III) occurs at nano- to sub-nanomolar ranges (Moore et al., 2001), and a majority (99%) of Fe is organically bound (Wu et al., 1995). Iron uptake is particularly important for photosynthetic diazotrophs as they possess higher Fe: C ratios (Kustka et al., 2002). The nitrogenase enzyme for N 2 fixation has a high Fe requirement (Schoffman et al., 2016; Hutchins and Sañudo-Wilhelmy, 2022), and thus Fe availability often limits N 2 fixation (Moore et al., 2009). Fe acquisition and transport in eukaryotic and prokaryotic plankton is complex and largely uncharacterized except for a few model systems. We identified and chose a putative Fe transporter, exbB, that is part of the exbB-exbD gene cluster that encodes the plasma membrane ExbB-ExbD complex (Hantke and Zimmermann, 1981; Stevanovic et al., 2012; Jiang et al., 2015), as at least one copy was present in all three of the Richelia symbionts. We expected the exbB expression to be up-regulated with N 2 fixation. Indeed, in the North Atlantic transect, expression for exbB and nifH was significant and positively correlated in both Richelia strains (RintRC01, ReuHH01).

Gene copy number and expression as indicators of nutrient status

In general, cyanobacteria (and other phytoplankton) have evolved different strategies for storing and surviving periods of low concentrations of nutrients (Hudek et al., 2016). One such mechanism involves possessing multiple systems for nutrient transport and/or import, which can be up or down regulated depending on the external nutrient concentration. For example, four copies of the pstS gene were identified in the facultative (RrhiSC01) symbiont and likely reflects a strategy to deal with the low P concentrations in the ambient water. For comparison, the genomes of several strains of the unicellular diazotroph Crocosphaera possess multiple and varying numbers of pstS (3–7 copies) (Bench et al., 2013). However, the large colonial filamentous diazotroph Trichodesmium erythraeum IMS101 which also co-occurs, only has one copy of pstS (Orchard et al., 2009). Similar to pstS, exbB exists in multiple copies in the genome of RrhiSC01. Although it is possible that the additional copies found in RrhiSC01 are due to duplications that occurred after speciation. The more likely scenario is that the three and four gene copies of exbB and pstS found in RrhiSC01, respectively, is the ancestral state. The genomes (albeit draft) of the other Richelia symbionts RintRC01 and ReuHH01 are reduced and lack multiple genes in comparison to RrhiSC01. Thus, the lower gene copy numbers in RintRC01 and ReuHH01 are not unexpected and likely reflects a redundancy of the lost genes. Residing in the cytoplasm of the host (Caputo et al., 2019; Flores et al., 2022), the endobiont ReuHH01 is not exposed to the ambient water and thus limited to nutrient concentrations within the host cell. Likewise, the intermediate endosymbiont RintRC01 is restricted to availability of nutrients in the host’s periplasmic space (Villareal, 1992). The lower retention of genes related to Fe and P scavenging could suggest that the endobiont ReuHH01 is acquiring nutrients, either through competition or through direct exchange with the host.

In environmental studies, functional gene expression is often used as indicators of activity (e.g., nifH for N 2 fixation) and limitation. For example, in some organisms pstS is up-regulated at a threshold concentration of dissolved inorganic P (DIP) (Pereira et al., 2019). Similar evidence presented here showed that pstS expression by both RintRC01 and ReuHH01 was highest in surface water samples (5–20 meters) when DIP was at its lowest measured concentration (e.g. 0.19- 0.20 µmol L -1 at 20 m). Additionally, pstS expression by RintRC01 was negatively correlated with DIP concentration. Likewise, nifH expression and consequently N 2 fixation is tightly coupled to low or often immeasurable concentrations of DIN (Knapp, 2012); in heterocystous cyanobacteria, heterocyst differentiation is initiated by the depletion of fixed N (Peterson and Wolk, 1978). Likewise, the environmental nifH expression by ReuHH01 and RintRC01 was lowest at depth where DIN was highest and for at least RintRC01 was negatively correlated with TIN concentration. When subjected to low Fe concentrations, the cyanobacteria strains Anabaena sp. PCC 7120 and Synechocystis PCC 6803 show higher expression of exbB (exbB3 and exbB1 - exbB3, respectively) and inactivation of the exbB results in Fe starvation (Stevanovic et al., 2012; Jiang et al., 2015). Moreover, in Anabaena sp. PCC 7120, expression of exbB1 is also regulated by growth stage and N limitation with a 12- fold increase in expression when grown without N, and higher expression in early and exponential growth stages (Stevanovic et al., 2012). In comparison to our diel experiment, where RrhiSC01 was growing exponentially and under constant N limitation, the expression of exbB was highly variable and did not show a significant periodicity. Wild populations of RintRC01 also express exbB without significant periodicity (Harke et al., 2019).

In summary, both pstS and exbB are potentially good candidates for indicating nutrient status, e.g., P and Fe limitation, respectively, for the various Richelia spp. symbionts. Additionally, future work could consider to compare the ratio of expression to presence (e.g. RNA: DNA) for a particular gene target to assess the overall activity of a strain.

Low rates of N 2 fixing activity consistent with low abundances of diazotrophs

The North and South Atlantic have been consistently surveyed over the past 50 years for bulk N 2 (and C) fixation. These represent regions with ideal conditions for diazotrophy: low DIN, sufficient P, and seasonal inputs of Fe by aerosol deposition (e.g., Saharan dust storms) and/or riverine (Orinoco, Amazon) inputs (Benavides and Voss, 2015; Luo et al., 2012; Shao et al., 2023). Prior to the field campaign, unusually high deposition of dust was recorded (Matzenbacher et al., 2024), yet other conditions, e.g., measurable DIN, were not conducive to high diazotrophic activity. The latter is also consistent with the negative correlation of nifH expression by both Richelia strains and DIN concentrations. Furthermore, both the nifH expression by Richelia strains and measured bulk N 2 fixation was consistently low, and observations of other common cyanobacteria N 2 fixers (Trichodesmium, Crocosphaera) were additionally sparse (personal observation). Despite measuring low N 2 fixation, our data contribute new rate measurements for the Atlantic Ocean and are consistent with low rates reported for non-cyanobacterial diazotroph dominated communities (Moisander et al., 2017).

The highest rates were recorded at station 17 in the western region of the transect and are consistent with the microscopy observations and qPCR detections for co-occurrence of both diatom- Richelia symbioses, albeit at lower densities (10–20 cells L -1). Both Hemiaulus- and Rhizosolenia-Richelia symbiosis have high N 2 fixation rates (e.g., 0.1-4.29 and 4.1-4.4 fmols N cell -1 h -1, respectively) (Foster et al., 2021), and therefore even at low densities can contribute to the bulk rates that were measured.

Several experiments were amended with a DOP cocktail to assess the impact on the N 2 fixation community. According to genome content and growth studies of cyanobacterial diazotrophs, some can metabolize different DOP (e.g., phosphonate) compounds as a source of P (or C or N) (Dyhrman and Haley, 2005; Dyhrman et al., 2006). Of the different Richelia symbionts, only the genome of the external facultative symbiont RrhiSC01 encodes a phosphonate transporter and some homologues for genes encoding the C-P lyase pathway which is necessary for hydrolysis and assimilation of phosphonate (Nieves-Morión et al., 2020). It was however unexpected that the amendment would have a significant negative impact on the bulk N 2 fixation. A similar decrease in bulk N 2 fixation activity was reported previously in experiments carried out in the Northern Red Sea with another DOP (glycerol phosphate) addition (Foster et al., 2009). One explanation for the decrease is that the amendment of DOP increased competition in the experimental bottle leading to low N 2 fixing activity, or perhaps contained trace amounts of N and thus negatively impacted diazotrophy.

Decades of research has been devoted to understanding the occurrence, distribution, and activity of different N 2 -fixers in the global oceans. With this work we provide new methods for specific and sensitive detection of three biogeochemically important marine N 2 -fixing Richelia spp. symbionts of diatoms. The nine new qPCR assays were experimentally tested in lab and field samples. The assays target three gene orthologs that when expressed can be indicative of metabolic activities between various Richelia spp.: phosphorus transport (pstS), iron transport (exbB), and a constitutively expressed gene for normalization (rnpA). Our results show that these assays are specific, can detect temporal dynamics in gene expression of lab cultures, and importantly are potentially indicative of nutrient limitation (e.g., P, Fe). The next obvious steps would be to apply the assays in controlled experiments under varying nutrient concentrations in order to identify the threshold nutrient concentration for upregulation of gene transcription. However, the latter requires cultivation of the various symbioses, which is currently unsuccessful. Finally, the assays are an improvement over the previously used nifH assays as estimations of abundance (DNA based) and metabolic activities (RNA based) were without cross-reaction.