Trichodesmium erythraeum IMS101 (Trichodesmium) was grown on an orbital shaker (150 rpm) under a 12h dark: 12h light cycle (ca. 130 μmol photons m^-2 s^-1) in modified YBC-II medium (Polyviou et al., 2015; Snow et al., 2015a) having a final EDTA concentration of 20 μM (Hong et al., 2017). Experimental cultures at a total volume of 40 ml were initiated from an Fe replete stock culture rinsed with no-added Fe media before inoculation and incubation in triplicate under 4 conditions: (Fe-) with no added Fe, (Fe+) with 400 nM added FeCl₃-EDTA as a source of dissolved soluble Fe (sFe), equivalent to a calculated bioavailable inorganic Fe concentration of 1100 pM Fe’ (Snow et al., 2015a) and two treatments with 0.25 mg ml^-1 of Saharan dust, collected between 04°43’N, 28°55’W and 06°56’N, 28°07’W using a mesh system aboard the RSS Shackleton (Murphy, 1985). Dust was added either directly (Dust+) or inside a barrier to the culture created using 8 kDa MWCO dialysis tubing (DT) (BioDesign, Carmel, NY, United States) ([Dust]). For consistency and to confirm that diffusion of sFe through the membrane was possible the FeCl₃-EDTA was also released into the Fe+ media through the DT. DT was washed by boiling in 120 mM Na₂HCO₃ and subsequently in 10 mM Na₂EDTA and 10 mM NaOH twice before rinsing three times in milli-Q H₂O. Clean DT, handled with 10% HCl washed plastic tweezers, was included in all experimental treatments except for a no-DT control. The no-DT control was grown under no added Fe conditions (i.e., as per the Fe- treatments) with the physiology of Trichodesmium showing no significant differences when compared to the Fe- treatment, suggesting that inclusion of DT did not result in any significant Fe contamination (Supplementary Figure S1). Within a separate control experiment the physiological responses of Trichodesmium were monitored to test whether nutrients from dust dissolved preferentially when dust was free in the culture vessels compared to when it was enclosed in DT (Supplementary Figure S2). Within this control experiment, the same series of treatments (i.e., +Fe, -Fe, Dust+ and [Dust]) were incubated abiotically, i.e., without the addition of Trichodesmium for a period of 14 days before filtration through a 0.2 μm filter into clean culture flasks and subsequent inoculation of the filtrate (i.e., with the dust particles removed) with equal concentrations of cells.

Trichodesmium growth was monitored by cell counts performed using a Sedgewick Rafter counting chamber and a GX CAM-1.3 camera on an L1000A biological microscope (GT Vision Ltd., Suffolk, United Kingdom). Cell and filament lengths were identified using GX capture (GX 14 Optical, Suffolk, United Kingdom) and ImageJ (Schneider et al., 2012). The concentration of cells was calculated for the duration of the experiment (Supplementary Figure S3) by dividing the total filament length by the average cell length and growth rates calculated using the gradient of the natural logarithm of measurements during exponential growth.

Photosynthetic physiology was monitored using a FASTtracka^TM MkII Fast Repetition Rate fluorometer (FRRf) integrated with a FastAct^TM Laboratory system (Chelsea Technologies Group Ltd., Surrey, United Kingdom). Measurements were made 2.5 h after the beginning of the photoperiod following dark adaptation of samples for 20 min and consequent exposure to a background irradiance of 29 μmol photons m^-2 s^-1 for 2–5 min (Richier et al., 2012). F_v/F_m was used as an estimate of the apparent PSII photochemical quantum efficiency (Kolber et al., 1998). Data presented is the average of three technical replicate measurements for each of the three biological replicates.

Trichodesmium cultures were filtered onto GF/F filters on day 17 of the experiment and snap frozen in liquid nitrogen. RNA was subsequently extracted using the RNeasy Plant Mini Kit (Qiagen, Manchester, United Kingdom), treated with TURBO DNA-free^TM DNase (Life Technologies Ltd., Paisley, United Kingdom) and stored at -80°C. Absence of genomic DNA contamination was confirmed by PCR.

Paired end libraries for Illumina sequencing were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina), following the manufacturer’s low-throughput protocol. This was modified for use with bacterial RNA samples by replacing the initial poly-A-based mRNA isolation step with ribosomal RNA (rRNA) depletion using a bacterial Ribo-Zero rRNA Removal Kit (Illumina Inc., San Diego, CA, United States). Libraries were pooled and sequenced on an Illumina MiSeq instrument using paired-end sequencing, with a read length of 151 bp.

The raw sequence reads were pre-processed with CutAdapt (v1.8.1, Martin, 2011) to remove TruSeq adapter sequences and low quality bases from the 3′ end of reads, using a quality threshold of 15 and minimum read length of 36 bp. Trimmed reads were mapped against the Trichodesmium erythraeum IMS101 reference genome assembly (Ensembl Bacteria database, accession GCA_000014265.1) using TopHat v2.0.14 (Kim et al., 2013), and reads mapping to annotated genes in IMS101 were counted using HTSeq, version 0.6.0, (Anders et al., 2015). The –library-type fr-firststrand and -stranded = reverse options were used in TopHat and HTSeq, respectively, to correctly account for strand-specific read mapping. Libraries were normalized based on the total library size using the median-of-ratios method that is implemented in DESeq2 (Love et al., 2014).

Differential gene expression analysis was carried out using the Bioconductor package DESeq2 version 1.8.2 (Love et al., 2014), running on R version 3.2.0. A set of genes with expression significantly affected by growth condition was identified using an ANODEV approach implemented in the DESeq2 package. The fit of read count data against a one-factor negative binomial general linear model, where replicates were grouped by growth condition, was compared (using likelihood ratio testing) to a reduced model where condition information was removed. Pairwise contrasts between conditions were calculated using a Wald test. In both cases, genes with an adjusted P-value of < 0.05 following Benjamini–Hochberg correction for multiple testing were classed as significantly differentially expressed.

The results of the RNA sequencing (RNAseq) analysis were validated via quantitative reverse transcription (RT) Polymerase Chain Reaction (qPCR) using the cytochrome c₅₅₃ gene (cyt c₅₅₃, Tery_2561) the iron-stressed-induced protein-A gene (isiA, Tery_1667) and the Mo-dependent nitrogenase-like-protein gene (nif-like, Tery_4114). Patterns of expression from the two methods closely resembled each other and regression analysis indicated a good correlation (R² = 99%) (Supplementary Figure S4).

To group differentially transcribed genes by transcription profile across treatments, the Pearson product moment correlation coefficient (PPMCC) was calculated for each pair of genes from regularized log (rlog) transformed read count data. Hierarchical clustering was then carried out using the WPGMA method, with 1-PPMCC as the distance metric. Similarly, to cluster samples based on overall transcription profile, pairwise sample distances were calculated using correlation coefficients and hierarchical clustering was carried out using the average linkage method (Minitab 17.3.1, Minitab Inc., Coventry, United Kingdom).

The effect of Saharan desert dust addition on the physiology of non-axenic Trichodesmium erythraeum IMS101 (hereafter Trichodesmium) when in the direct vicinity (Dust+) or separated ([Dust]) from the cells within otherwise low Fe medium was compared to iron deplete (Fe-) and iron replete (Fe+) treatments (as described in Materials and Methods). Growth of Trichodesmium was significantly faster under Fe+ (0.31 day^-1) and Dust+ (0.28 day^-1) compared to Fe- (0.16 day^-1) and [Dust] (0.17 day^-1) conditions (Figure 1A and Supplementary Figure S3) (GLM and Tukey test; P < 0.05). Additionally, while photosynthetic efficiency (F_v/F_m) declined rapidly in Fe- and [Dust] treatments from day 6 until the end of the experiment, it was maintained at significantly higher levels in Fe+ and Dust+ treatments (Figure 1B). In contrast, although the Fe+ media similarly supported a higher growth rate within the abiotic control experiment, growth rates within media preincubated with dust either inside or outside the DT were statistically indistinguishable from both each other and the -Fe treatment (Supplementary Figure S2).

Samples collected on the last day of the main experiment (day 17) were analyzed using RNAseq to ascertain the gene transcription profiles relating to each experimental condition. The Trichodesmium erythraeum IMS101 genome is composed of 4451 annotated genes, ∼54% of which are annotated with a Gene Ontology (GO) classification number (GO annotated) (Figure 2). The majority of Trichodesmium genes (∼98%) were identified (read counts/million > 0) to be transcribed in at least one treatment and 4076 (∼92%) were identified in all four treatments. Of these, 3257 genes had read counts above the statistical threshold required for inclusion in differential gene transcription analysis, as determined by the DESeq2 Bioconductor package (Love et al., 2014). Approximately one third of those genes (1050) were significantly affected by growth condition (adjusted P-value < 0.05 following Benjamini–Hochberg correction for multiple testing, DESeq ANODEV) (Figure 2).

The validity of experimental biological replicates was confirmed through hierarchical clustering of gene transcription profiles which revealed closer intra-treatment compared to inter-treatment similarities (Figure 3). In addition, two distinct groups of treatments were resolved, highlighting similarities in gene regulation between [Dust] and Fe- treatments (group 1) and between Fe+ and Dust+ treatments (group 2) (Figure 3A).

A heat map illustrating all 1050 differentially regulated genes identified in the transcriptomic analysis revealed 10 distinct patterns of gene transcription from the clustering analysis (Figure 3B and Supplementary Figure S5) demonstrating significant changes in gene expression across all treatments. The majority of genes (42%) fall in clusters regulated similarly between Fe+ and Dust+ treatments with either increased (Cluster 2, n = 185) or reduced (Cluster 10, n = 259) transcription relative to Fe- and [Dust] (Supplementary Figure S5). Differential gene regulation in Fe+ compared to Fe-, Dust+ and [Dust] treatments (Supplementary Figure S5) is also observed for a large number of genes (Cluster 4, n = 109 and Cluster 6, n = 135).

The expression of a set of 12 genes (Supplementary Table S1) of known importance during Fe stress conditions and/or previously identified to be differentially expressed (at the transcript or protein level) under variable Fe concentrations was examined in order to access the molecular response to Fe specifically. The genes include the chlorophyll-binding antenna (isiA, iron stress induced protein A gene, Tery_1667) (Bibby et al., 2001; Shi et al., 2007; Richier et al., 2012; Snow et al., 2015a), the flavodoxin genes (fld1, Tery_1666; fld2, Tery_2559) which are an Fe free alternative to ferredoxin (LaRoche et al., 1996; Chappell and Webb, 2010), the interchangeable plastocyanin (copper binding) and cytochrome c₅₅₃ (Fe binding) genes (petE, Tery_2563 and petJ,Tery_2561, respectively) (Wood, 1978; Peers and Price, 2006; De la Cerda et al., 2007), genes for the Fe storage proteins bacterioferritin (bfr, Tery_2787) and ferritin (ftn, Tery_4282) (Keren et al., 2004; Marchetti et al., 2009), the ferric uptake regulators (fur1, Tery_1958; fur2- Tery_3404 and fur3, Tery_1953) (González et al., 2012, 2016), the Fe binding nitrogenase NifH subunit (nifH, Tery_4136) (Shi et al., 2007; Richier et al., 2012; Snow et al., 2015a) and fructose bisphosphate aldolase class II (fbaA, Tery_1687) (Snow et al., 2015a).

Four of these twelve genes, petJ, petE, fur1 and fld2, show significant regulation both in response to direct dust addition (Dust+ compared to [Dust]) and increased availability of dissolved Fe (Fe’) (Fe+ compared to Fe-) while a further four genes, fbaA, isiA, fur2 and fld1, were only significantly regulated through variable Fe’ (Figure 4A). Genes bfr, ftn, fur3 and nifH were not found to be differentially regulated by any experimental condition.

Analysis of a further set of genes (Supplementary Table S2) was targeted to better understand mechanisms of Fe acquisition under the differing growth conditions. Suggested mechanisms for Fe reduction by cyanobacteria include the transfer of electrons to extracellular Fe³⁺ through the plasma-membrane-located alternative respiratory terminal oxidase (ARTO) (Kranzler et al., 2014) and potentially electrically conductive pili (Lamb et al., 2014). Genes ctaC (Tery_0278), ctaD (Tery_0277), and ctaE (Tery_0276) encoding subunits of ARTO displayed elevated transcription in Fe- compared to Fe+ and [Dust] compared to Dust+ (Figure 4B) revealing regulation by Fe’ and dust additions within the vicinity of the cells. Interestingly, transcripts of these genes were significantly downregulated in Dust+ compared to all other treatments. The pilA (Tery_2388) gene encoding the pili protein A was not differentially transcribed as a function of Fe’ (Fe+ compared to Fe-) and although transcription was reduced by dust (compared to Fe-), this was only significant when dust was separated from the cells (Figure 4B).

Well-characterized cyanobacterial Fe uptake pathways include the FutABC and FeoABC systems for transport of ferric (Fe³⁺) and ferrous (Fe²⁺) iron, respectively, across the inner membrane (Kammler et al., 1993; Katoh et al., 2001a,b). The futA (Tery_3377) gene, annotated as idiA (iron-deficiency induced) and encoding the Fe-binding subunit of the Fut Fe³⁺ transporter (Polyviou et al., unpublished), was downregulated in response to increased Fe’ (Figure 4B). Similarly, the transcript abundance of feoB (Tery_2878) was reduced in Fe+ cultures compared to the other treatments although the difference was not statistically significant.

Translocation of organically complexed Fe to the periplasm in cyanobacteria happens through TonB dependent transporters (TBDTs) at the outer membrane which seem to require energy transduction by TonB and ExbBD (Larsen et al., 1999; Brinkman and Larsen, 2008; Ollis et al., 2009; Ollis and Postle, 2012; Ollis et al., 2012), although ExbB and ExbD have also been previously linked to inorganic Fe uptake in Synechocystis sp. PCC 6803 (Jiang et al., 2015). Transport of organically complexed Fe across the inner membrane employs Fec/FhuD systems (Krewulak and Vogel, 2008; Stevanovic, 2015). The only Trichodesmium homolog to fhu/fec genes, fhuD (Tery_3943) (encoding for a periplasmic binding protein homolog), was downregulated by increased Fe’ and dust when added to the immediate cellular environment (Figure 4B). Homologs of tonB (Tery_1560, Tery_2593) and exbBD (exbB: Tery_1868, Tery_4448; exbD: Tery_4449) were not identified to be significantly differentially regulated, although transcription appeared higher in the Fe+ and Fe- treatments compared to the [Dust] and Dust+ treatments (with the exception of Tery_1868 for which read counts were below the threshold for inclusion in the analysis) (Figure 4B).

Members of a previously identified putative siderophore production and uptake pathway, proteins Tery_3824–3826 (Snow et al., 2015a) were observed to be significantly downregulated under increased Fe’ (significant difference for Tery_3825 and Tery_3826) (Figure 4B). Finally, transcription of the heme oxygenase homolog (hmox, Tery_0335), involved in extraction of heme-bound Fe was stimulated by the dust treatments compared to Fe+ and Fe- irrespective of the separation from the cells using the DT (Figure 4B). Some indication of increased transcript abundance due to increased Fe’ was also observed although this was not statistically significant.

One third of the differentially regulated genes in the analysis were not GO annotated. Amongst them was a series of genes encoding for proteins with Haemolysin-Type Calcium Binding (HTCaB)-like domains (Supplementary Table S3) members of which were identified amongst the most differentially regulated genes in this study. Of 10 genes annotated as having a HTCaB-like domain and identified as regulated by some experimental condition in our analysis, Tery_0419, Tery_0424, Tery_2055 and Tery_1355, are regulated by Fe’ and when contact with dust is permitted and Tery_3467 and Tery_2710 by Fe’ only (Figure 5A). Bioinformatic analysis indicates that HTCaB genes include multiple CHRD domains (identified in chordin) expected to have an immunoglobulin-like β-barrel structure based on some similarity to superoxide dismutases (Hyvönen, 2003). Among these genes, Tery_3467 which was regulated by Fe’ has been previously identified as the zinc binding alkaline phosphatase (APase) gene phoA (Orchard et al., 2003).

Also among the uncharacterized genes is a cluster which spans eight genes (Tery_0843–Tery_0850) all of which were strongly differentially regulated in response to Fe’ availability and dust addition to the direct cellular vicinity (Figure 5B) with five (Tery_0850, Tery_0845, Tery_0846, Tery_0849, Tery_0847) amongst the top twenty most differentially transcribed genes in the RNAseq analysis. Tery_0848 and Tery_0849 are annotated as cell surface proteins while Tery_0843 is recognized by the UniProtKB Automatic Annotation pipeline as a membrane protein (Supplementary Table S3). In addition, Tery_0844 is predicted to contain an iron-sulfur binding site and Tery_0845 belongs to the heme oxygenase superfamily. The fifth gene in this cluster (Tery_0847) encodes for MetE (5-methyltetrahydropteroyltriglutamate-homocysteine S-methyltransferase) which catalyzes the transfer of a methyl group to and from methionine (Supplementary Table S3).

Trichodesmium has a large genome, with abundant non-protein coding regions, and has recently been shown to exploit sophisticated gene regulation (Pfreundt et al., 2014; Pfreundt and Hess, 2015; Walworth et al., 2015), as well as being polyploid (Sargent et al., 2016). The Trichodesmium genome harbors 17 group-II introns interrupting a total of 11 genes, 7 of which (Tery_0428, Tery_4732, Tery_4799, Tery_3633, Tery_2080, Tery_1635, and Tery_0008) are differentially regulated in one or more of the experimental treatments (Figure 6). Genes Tery_0428, Tery_4732, and Tery_4799 are regulated only when dust is added directly to the cell environment, while Tery_2080 is expressed more highly in both Dust+ and [Dust] compared to the Fe- and Fe+ treatments. Significant Fe’ regulation is only observed for Tery_3633.