The relationship between toxicity and mixotrophy in bloom dynamics of the ichthyotoxic Prymnesium parvum

Konstantinos Anestis^1*, Sylke Wohlrab^1,2, Elisabeth Varga³, Per Juel Hansen⁴ and Uwe John^1,2*

¹ Alfred Wegener Institute, Helmholtz Center for Polar and Marine Research, Am Handelshafen 12, 27570 Bremerhaven, Germany

² Helmholtz Institute for Functional Marine Biodiversity at the University of Oldenburg (HIFMB), Ammerländer Heerstraße 231, 26129, Germany

³ Department of Food Chemistry and Toxicology, Faculty of Chemistry, University of Vienna, Währinger Straße 38-40, 1090 Vienna, Austria

⁴ Marine Biology Section, University of Copenhagen, Strandpromenaden 5, 3000 Helsingør, Denmark

*Corresponding authors:

Konstantinos Anestis (kanestis@awi.de)

Uwe John (uwe.john@awi.de)

Keywords: HABs; haptophytes; marine toxins; polyketides; microbial ecology; species interactions

ABSTRACT

Toxin production in planktonic protists is widespread and comprises an effective mechanism to eliminate competitors or grazers. Toxins released into the water can mediate the immobilization of prey for subsequent consumption or can mediate the lysis and uptake of the released nutrients. For the worldwide-distributed ichthyotoxic and mixotrophic haptophyte Prymnesium parvum , the relationship between toxin production; impact on co-occurring species, and mixotrophy remains blurred. In the current study, we show that changes in salinities (5 vs 30), phosphorus (P) availabilities (P-replete vs P-deplete) and cell densities affect growth, toxicity and mixotrophy in P. parvum . Cell density positively affected cellular toxin content by a factor of up to 10. Low salinity resulted in a higher mortality of the cryptophyte prey Teleaulax acuta , with ~80% of cells being lysed after 2 h of incubation. However, phagotrophic rates were higher in P-deplete conditions, independent of the salinity. Transcriptomic analysis of the monocultures revealed the up-regulation of genes involved in endocytosis under either low salinity and phosphorus, suggesting that this process is evolutionarily conserved, triggered by environmental stressors and independent of prey presence. Polyketide synthase genes, potentially involved in toxin biosynthesis, exhibited distinct expression patterns, depending on the physiological status, toxicity and with generally higher expression under the high cell density conditions. Overall, our study contributes to a better understanding of the dynamics between the two critical processes of toxin production and mixotrophy, and has important implications for bloom formation and its maintenance in this ecologically important species.

INTRODUCTION

Harmful algal blooms (HABs) are characterized by a rapid proliferation with detrimental effects on the ecosystem. The causes of such plankton blooms are variable and can be of both anthropogenic and environmental origin (Lewitus et al. , 2012). Among the best studied driving factors is eutrophication due to the inflow of dissolved nutrients, which can be subsequently used by phototrophic plankton (both prokaryotic and eukaryotic) for rapid growth (Heisler et al. , 2008). The increasing global occurrence of HABs (Anderson et al. , 2012; Gobler et al. , 2017) has additionally been suggested to be a consequence of climate change (Gobler, 2020).

Haptophytes are a diverse group of (nano)plankton with a worldwide distribution and important contributions to primary production and biogeochemical cycles (Edvardsen et al. , 2016). Known bloom-forming haptophytes include Chrysochromulina leadbeateri ,Prymnesium polylepis , Prymnesium parvum andPhaeocystis spp .. HABs often cover large coastal areas that cause mortalities of fish and other marine fauna and/or damage to ecosystem function, e.g., disruption of food webs or oxygen depletion, or loss of recreational opportunities due to biofouling of beaches and coastal waters (Karlson et al. , 2021). Massive fish mortalities may be attributed to known ichthyotoxins or unknown compounds (Andersenet al. , 2015; Mardones et al. , 2019). Several fish killing HAB events have been recorded world-wide in recent years (Hallegraeffet al. , 2021) and are often associated with haptophytes (Bresnanet al. , 2021; Karlson et al. , 2021; John et al. , 2022). In most cases in Scandinavia, massive fish-killing events have been directly linked to blooms of marine haptophytes, particularly members of the genera Prymnesium and Chrysochromulina . Many HAB forming species, including many haptophytes, are mixotrophic, i.e., they are able to combine the two trophic modes of phagotrophy/osmotrophy and phototrophy in order to cover their nutritional needs (Burkholder et al. , 2008; Unrein et al. , 2013; Flynn et al. , 2018). This is especially the case in eutrophic ecosystems, where growth rates of mixoplanktonic species (Flynn et al. , 2019) can be affected by both the increased dissolved nutrients and the availability of algal and bacterial prey (Burkholder et al. , 2008). Additionally, mixotrophy can be highly advantageous under low or imbalanced nutrient conditions, as it can be an efficient mechanism to compensate for limiting nutrients (Stoeckeret al. , 2017).

The ichythyoxic HAB species P. parvum is known to produce toxic compounds, collectively called prymnesins. These cause lysis/death of competitors and grazers, as well as lysis of fish gill cells. Prymnesins are large ladder-frame polyketide compounds, of which three types are described, A-, B- and C- types, each differing in the length of the carbon backbone (Igarashi et al. , 1999; Rasmussen et al. , 2016). The production of lytic compounds by P. parvum , whether studied using bioassays or by actual measurements of prymnesins, are known to be influenced by phosphorus (P) availability, light and temperature ((Beszteri et al. , 2012; Qin et al. , 2020; Taylor et al. , 2020, 2021; Medić et al. , 2022). Most probably, prymnesins are produced by polyketide synthases (PKS) of the modular type I (Anestis et al. , 2021). The release of lytic compounds by P. parvum into the surrounding water has been shown to assist feeding, via immobilization and lysis of the prey, which allows uptake of nutrients through either osmotrophy or phagotrophy (Skovgaard and Hansen, 2003; Tillmann, 2003).

P. parvum typically form blooms in estuaries at low salinity (2-8), despite the fact that this species grows fine in the salinity range of 0.5-45 (Edvardsen and Paasche, 1998; Barone et al. , 2010). Why it preferably blooms at low salinities is unknown, but could be due to input of organic material from freshwaters, or simply because very few species can thrive at these salinities (Brand, 1984). Additional factors that influence the bloom dynamics of this species include the mediation of viruses that interact with components of their cellular membrane (Wagstaff et al. , 2018, 2021).

In this study, P. parvum strain UIO223 was used to investigate the combinatory effects of bloom formation supporting conditions i.e., low salinity, P limitation, and cell density. More specifically, we used two salinity (salinity of 5 vs salinity of 30) and P (P-replete vs P-deplete) treatments to examine cellular toxin content and toxin profile, phagotrophy, and metabolic processes related to adaptation under the different treatments. The metabolic processes were related to growth, potential to mixotrophy (expression of phagotrophy related genes in monocultures) and biosynthesis of polyketides. Additionally, we sampled at two cell densities (low vs high) to investigate the processes that support the maintenance of a bloom with high cell density. We hypothesized that A) P. parvum will increase its toxin content and production at low salinity, and low P concentrations in combination with increased phagotrophy. B) cellular and metabolic processes deduced from gene expression analyses will depict the corresponding cellular adjustment; C) P. parvum needs to reallocate cellular resources to maintain growth under these limiting conditions under high cell density (bloom concentration), which has an impact on levels of phagotrophy.

MATERIALS AND METHODS

Experimental set-up, sampling and incubations with prey

Phosphorus replete cultures of P. parvum strain UIO-223 were grown in standard K-medium with PO₄^3-concentration of 36 μM, while for the P deplete conditions, PO₄^3- was added at a final concentration of 2.4 μM. Culture medium with a salinity of 5 was obtained by diluting sterile filtered seawater with MilliQ water. The concentration of inorganic carbon in the low salinity medium was restored by adding 1 M NaHCO3 in 1 mL L^-1 of medium. Prior to establishing the experimental cultures, the culture was rendered axenic using a cocktail of antibiotics (165 μg mL^-1 ampicillin, 33.3 μg mL^-1gentamicin, 100 μg mL^-1 streptomycin, 1 μg mL^-1 chloramphenicol, 10 μg mL^-1ciprofloxacin). The treatment with antibiotics lasted 4 days and was performed twice with an interval of one week. The axenicity of the cultures was examined by fluorescence microscopy after staining with 4′,6-diamidino-2-phenylindole (DAPI).

Four replicate cultures were established for each combination of salinity (30 and 5) and phosphorus (P replete and P deplete; 2x2 factorial design). All cultures were kept in 17 °C in 16:8 light:dark cycle and under a photon flux density of 80 μmol photons m^-2 s^-1. The cultures were bubbled to avoid carbon deficiency. The pH of the cultures was monitored daily, and it never exceeded 8.5. All cultures were sampled at both low cell density and high cell density.

Cell enumeration was performed using a Multi-Sizer III particle counter (Beckman-Coulter, Fullerton, USA). Phagotrophy was estimated usingTeleaulax acuta (SCCAP K-1486) as prey, grown under the same condition as P. parvum . The initial cell concentration for the incubation was 30*10³ cells mL^-1 and the ratio between P. parvum and prey was 1:1, following the recommendation of the existing literature (Lundgren et al. , 2016). To estimate phagotrophy, 2 mL of culture were fixed with Lugol (final concentration of 2%, v:v), and inspected in a Sedgewick Rafter chamber under an Axio Vert A1 Microscope equipped with a Colibri 7 (Zeiss) light source.

Physiological parameters sampling and processing

Samples for inorganic nutrient measurements were taken by filtrating 15 mL of culture using 0.2 μM Millipore filters in order to eliminate cells. The nutrients that were measured included nitrate, nitrite, and phosphate, and were analyzed with a continuous-flow autoanalyzer (Evolution III, Alliance Instruments, Freilassing, Germany). The protocols that were used are standard for the quantification of nitrite/nitrate (Armstrong et al. , 1967) and phosphate in seawater (Eberlein and Kattner, 1987).

Samples for POC and PON analysis were collected by filtering 30 mL of culture through precombusted glass microfiber filters (Whatman GF/F, Maidstone, UK; nominal pore size: 0.7 µm) and were immediately frozen in pre-combusted glass vials until further analysis. The filters were dried at 50 °C overnight, acidified with 300 µL 0.2 N HCl, and again dried overnight at 50 °C. The acidified and dried filters were packed in tin foil and analyzed on a Euro Elemental Analyzer 3000 CHNS-O (HEKAtech GmbH, Germany).

Samples for Chl-a concentrations were taken by filtering 15 mL of culture through 22 mm glass microfiber filters (Whatman GF/F, Maidstone, UK; nominal pore size: 0.7 µm). The sample filters were frozen at −80 °C until laboratory analysis. Extraction was performed by sonication of filters in 10 mL 90% acetone and then incubated overnight in darkness at 4 °C. The extract was centrifuged at 3020 × g for 10 min, and the fluorescence of the supernatant was determined at 665 nm (TD-700 fluorometer, Turner Designs, Sunnyvale, USA).

Toxin extraction and quantification was performed according to a previously described protocol (Svenssen et al. , 2019), with small modifications as described in Anestis et al. , 2021. In brief, cells were collected on 22 mm glass microfiber filters (Whatman GF/F, Maidstone, UK; nominal pore size: 0.7 µm). The biomass on each filter was extracted two times with 20 mL MeOH each using an ultrasonic batch for 30 min with a centrifugation step in between (4300 ×g for 15 min at 4 °C). The combined extract (40 mL) was evaporated to dryness using a CentriVap Benchtop Vacuum Concentrator (Labconco Corporation, Kansas City/ MO, USA) at 35 °C. The samples were reconstituted with 1 mL methanol:H₂O (90:10, v:v) and short-time ultrasonic bath treatment. HPLC-FLD measurements were performed after derivatization with the AccQ-Tag Fluor Reagent Kit (Waters Cooperation, Milford/MA, USA) with a 1200 HPLC system (Agilent Technologies, Waldbronn, Germany) using fumonisins B1 and B2 as external calibrants due to the lack of standards and the obtained results are an approximation of the prymnesin content in the samples. To confirm the presence of prymnesins and identify the specific prymnesin analogues, HPLC-HRMS-measurements were performed using a 1290 UHPLC system coupled to a 6550 iFunnel QTOF LC/MS (both from Agilent Technologies). Chromatographic separation was achieved with a Kinetex F5 (2.1 × 100 mm, 2.6 μm, Phenomenex, Aschaf-fenburg, Germany) column using a water-acetonitrile gradient (eluent A: H₂0, eluent B: acetonitrile: H₂O (90:10, v:v)), both eluents contained 1 mM ammonium formate and 0.1% formic acid. The mass spectrometer was operated in the positive ionization mode in a scanning range of m/z 50 to 1700 with 3 scans per second.

RNA extraction, library construction and sequencing

Cells for RNA extraction were harvested by centrifugation at 1500 × g for 10 min and the cell pellet was transferred to 1 mL pre-chilled TriReagent mixed with glass beads. RNA isolation was performed as described in Wohlrab et al. , (2017). Libraries for sequencing were prepared using the Truseq Stranded mRNA Samples Prep LS Protocol (Illumina GmbH, Berlin, Germany) and 1 μg of RNA as input. The paired-end cDNA libraries (2x150 bp) were sequenced on the Illumina Nextseq 500 machine (Illumina, San Diego, USA) and a high-output kit v2 (2 x 150 cycles).

Quality control of sequencing data, de novo assembly and annotation

Prior to assembly, the raw reads were pre-processed using Trimmomatic v.0.39 (Bolger et al. , 2014) and reads contaminated with adapters or quality scores of <5 were trimmed using the default settings and only paired-end reads were retained. Inspection with FastQC assured the high quality of the reads and absence of adapter con-tamination before proceeding to the construction of the de novo assembly (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The reference transcriptome was constructed using Trinity (v2.11.0) with minimum contig length of 300 bp (Haas et al. , 2013). To reduce redundancy due to the assembling method, the output contigs were clustered together with a similarity threshold of 0.95 using the cd-hit-est command of CD-hit and the longest representative for each cluster transcript was kept (Fu et al. , 2012).

Gene annotation was performed using the Trinotate functional annotation suite (version 3.2.1; Grabherr et al., 2011). Transdecoder (v5.5.0) was used for detecting open reading frames (ORFs) for each gene (https://github.com/TransDecoder/TransDecoder). The gene sequences were aligned against the UniprotKB/Swiss-Prot (release 2021_03) using diamond with BLASTX and BLASTP as an option (blast version 2.10.1) and the best hit for each query was retained with an e value of <1 e^-9. Further annotation information, like entries from the KEGG and the assignment of Gene Orthology (GO) terms were retrieved from the UniprotKB database.

Gene expression analysis

Quantification of gene expression was performed by mapping the paired-end reads in CLC Genomics Workbench 20, using the default settings. The resulting gene count matrix was normalized using the variance stabilizing transformation function of the Deseq2 package (Loveet al. , 2014) and multiple comparison tests with three-way analysis of variance (ANOVA) were performed in order to define the effect of salinity, phosphorus and cell density in gene expression. For genes with significant interactions only, single t-tests were applied to identify the conditions causing expression changes. For every gene, a factor was considered to have a significant effect on the expression variance when satisfying an adjusted p value threshold of <0.01 and fold change >1.5. P values were adjusted by applying the Benjamini-Hochberg false discovery rate for multiple test correction (Benjamini and Hochberg, 1995). For the visualization of the intersecting sets of differentially expressed genes, the R package UpSetR was used (Conway et al. , 2017).

Statistical analysis

Statistical analysis of the physiological parameters was performed using three-way analysis of variance (ANOVA) and the factors salinity, phosphorus and cell density. The normality of the distribution of the residuals was checked with the Shapiro-Wilk test, and an appropriate transformation of the data (logarithmic or square transformation) was applied to achieve p > 0.05. Pairwise comparisons were examined using Tukey’s HSD (honestly significant difference) test. All analyses were performed using the R software.

RESULTS

Growth physiology

For the first 3 days of the experiment, growth rates for cultures at low salinity varied from 0.86 to 0.92 μ d^-1 in comparison to 0.72 to 0.75 μ d^-1 for cells growing at high salinity, irrespective of initial P treatment (p-value < 0.01; Table S1). Subsequently, growth rates gradually decreased for all treatments, with 0.16 to 0.19 μ d^-1 at the final day of the experiment for the non-limited cultures, while complete inhibition of growth was observed in the P-limited cultures (Figure 1a; Figure S1). In the low phosphorus (LP) treatments growth stopped at day 7, and the cell concentration remained stable for the following 2 days, indicating P limitation of growth (~450,000 cells mL^-1; Figure 1a-b). Extra P was added at day 9 to the LP treatments, which led to an immediate increase in cell concentration. The addition of extra P in this treatment was done to validate that P was the limiting factor and that growth was affected only by its depletion . In non-limited P cultures, cell concentration increased until the end of the experiment at day 13 and reached to 3.5 10⁶ cells mL^-1 and 3.9 10⁶ cell mL^-1 for high salinity and low salinity respectively (Figure 1a).

Particulate organic nutrients (C:N) were measured to gain understanding of the elemental composition of the cells under the different treatments. Under low cell density, the per cell quota in POC was affected by salinity with 59.9-61.2 pgC cell^-1 and 42.1-43 pgC cell^-1 for high and low salinity cells respectively (Figure 1c). When P depletion was reached these values were lower for the P replete cells with 16.5-21.6 pgC cell^-1, in contrast to 29.8-28.7 pgC cell^-1 observed in the P depleted cells (Figure 1c). All salinity, phosphorus and cell density treatments had statistically significant effects on the cellular particulate organic carbon (POC) content (Table S2). In general, cellular POC was higher in low cell density cultures for all treatments (pairwise comparisons with p < 0.005; Table S2). At low cell densities, the low salinity treatments had the lowest cellular POC content, while at the high cell densities, P-starved cells had a significantly higher cellular POC content. In accordance with the decrease in cellular POC content, a cell density dependent decrease in cellular particulate organic nitrogen (PON) content was also observed (Figure 1d). Salinity and P as single factors seemed to have no statistically significant effect on cellular PON content; significant differences were, however, observed in combinations with salinity and P. The decrease in cellular PON was higher than in cellular POC, indicated by the cellular C:N stoichiometry; high cell density cultures had a significantly higher C:N ratio (p <0.001; Table S2) compared to low cell density cultures (Figure 1e). Cells grown at low salinity, had consistently lower C:N ratios compared to those grown at high salinity. Cellular chlorophyll-a (Chl-a ) content was reduced by ~2.5-fold at the highly dense cultures as compared to low cell density cultures, this was consistent for all treatments (p <0.001; Table S2; Figure 1f). The combination of P-depletion and high cell density resulted insignificantly higher cellular Chl-a (p<0.001; Table S2).

Prey ingestion

Teleaulax acuta was provided as prey and its mortality and phagotrophy by P. parvum was measured as the percentage of cells containing food vacuoles after 2 h, 6 h, 24 h, 48 h and 72 h of incubation with the prey. T. acuta cells showed higher mortality when incubated with low salinity P. parvum cultures, and independently of phosphorus condition, with ~80% being lysed 2 h after the beginning of the incubation with P. parvum(Figure 2d). Decrease in T. acuta concentration at 2 h were observed also under the high salinity condition, but was considerably less extent than in low salinity (Figure 1c). For the low salinity/P replete treatments, complete prey removal was observed 24 h after the beginning of the incubations, while for the low salinity/P deplete treatments complete prey removal was observed at 72 h (Figure 2c). For the high salinity treatments, T. acuta mortality was lower, but complete prey removal was observed at 72 h in all cases. The growth of monocultures of T. acuta growing under the same culturing media used for P. parvum was monitored (Figure S2).

P. parvum cells containing a food vacuole were first observed after 6 h (Figure 2d). Phagotrophy was significantly induced in P starved conditions, independently of salinity, with an average of 11 ± 1 % of the total cells containing a food vacuole from 6 h until 48 h. In the P replete treatments, P. parvum phagotrophy was less with a maximum of < 6 % of the total cells feeding at some time point. The number of cells with a food vacuole considerably decreased for all treatments from 48 h to 72 h, consistent with the depletion of prey at 72 h. Interestingly, even though the highest mortality rates ofT. acuta were observed at the low salinity, in the P-replete condition, phagotrophy rates were very low. Growth rates recorded forP. parvum grown with prey were lower than the growth rates of monocultures for the high salinity treatments, which was not the case at the low salinity treatments, where growth rates of monocultures were similar to the cultures with prey addition (Table S1, Figure S3).

Toxin content and profile

Cellular prymnesin contents in P. parvum were significantly higher at high cell density compared to at low cell density for all treatments (p < 0.0001, Table S2). Cell density and P depletion both had a statistically significant effect in explaining the variance in cellular prymnesin contents (p < 0.001), as both factors led to increases in prymnesin content by ~8.3-fold and ~2.4-fold respectively (Figure 2a). Salinity had statistically significant effects only in combination with P concentration and cell density (Table S2). In low cell density culture conditions no significant differences in toxin contents were observed. The highest amount of cellular toxin content was observed among high density and P starved conditions (p <0.005).

The percentage contribution of the prymnesin analogs to the overall prymnesin quantity, varied among the treatments, and was affected by all factors of the experiment; salinity, P and incubation time (Figure 2b). The four analogs that were detected differ in presence or absence of attached sugars including a pentose, a hexose or both of these two. The main analog for the majority of time points and treatments was the prymnesin B-type containing one chlorine and one pentose moiety, with the exception of the low salinity/P-deplete cells for which the prymnesin without attached sugar accounted for about half of the overall toxins at high cell density. The prymnesin analog without sugar showed a general increasing pattern in the low salinity treatment and especially in high density. The prymnesin analog containing both pentose and hexose showed consistently low contribution to the overall prymnesin content, however, it accounted for 9 ± 1% and 7 ± 2% for the P starved cells under high salinity and low salinity, respectively. The hexose containing prymnesin analog was mainly present in the low salinity and P starved cells with a contribution of 11 ± 5%.

Gene expression analyses

The transcriptome of P. parvum consisted of 103,051 contigs, which resulted after their clustering with 95% similarity (Figure 3a). Gene expression differences were ana-lyzed via Principal Components Analysis (PCA) for assessing the driving experimental factors for changes in gene expression (Figure 3b). Salinity induced the big-gest change in the transcriptome under low cell density. For the high cell density cultures, P-replete and P-deplete culture samples formed two distinct clusters. Salinity explained the expression variance of 2233 transcripts, followed by cell density with 1622 transcripts, the interaction of salinity and P with 1041 and the interaction of salinity and cell density with 1125 transcripts (Figure 3c). Phosphorus as a sole parameter influenced the expression of 203 transcripts, while the number of transcripts increased to 988 for the interaction between P and cell density. The number of differentially ex-pressed genes (DEGs) with a Kyoto Encyclopedia of Genes and Genomes (KEGG) identifier for which salinity had an effect was 324, followed by cell density with 309 and P with 225 (Figure 4a-c).

Growth and cell density effects

Carbohydrate metabolism, including anabolic and catabolic processes, are significantly influenced by the salinity treatment, with genes involved in catabolic processes being significantly up-regulated under low salinity (Table 1). Significantly overexpressed genes under low salinity included glyceraldehyde 3-phosphate dehydrogenase (contig51264_c0_g1_i1) and pyruvate dehydrogenase E1 component (contig50014_c0_g1_i1), which showed 31-fold and 30-fold increase respectively. Glyceraldehyde 3-phosphate dehydrogenase catalyzes the conversion of glyceraldehyde 3-phosphate to 1,3-bisphosphoglycerate which can subsequently form pyruvate. Pyruvate in turn can be converted to acetyl-CoA and the production of energy in the form of NADH. Acetyl-CoA can subsequently enter the citric acid cycle and produce further energy. Additional overexpressed genes that are involved in glycolysis included pyruvate dehydrogenase E2 component (contig6111_c0_g1_i19) and glucose-6-phosphate 1-epimerase (contig1336_c0_g1_i33) with 2-fold and 6-fold increase respectively. Enzymes involved in valine, leucine and isoleucine degradation, such as isovaleryl-CoA dehydrogenase (contig51612_c0_g1_i1) and malonate-semialdehyde dehydrogenase (contig37437_c0_g1_i1) were up-regulated by 77-fold and 46-fold respectively.

On the other hand, malate synthase (contig7939_c0_g1_i20), which participates in the anabolic process of glyoxylate cycle, showed a 28-fold decrease under low salinity, highlighting the re-arrangement of cellular metabolism in favor of energy production rather than energy storage. Down-regulation under low salinity was also observed for genes encoding for glucokinase (contig9551_c0_g1_i5), phosphoglucomutase (contig10485_c0_g1_i71) and phosphoglycerate kinase (contig4027_c0_g1_i11).

P-starvation and cell density induced changes in the expression of genes involved in energy metabolism (Figure 4a). In P depleted cells, ten genes involved in energy metabolism were down-regulated while five genes were up-regulated. Processes related to translation (38 down-regulated vs 9 up-regulated) and transcription (18 down-regulated vs 3 up-regulated) were strongly down-regulated under P depletion (Figure 4a). The cell density effect in growth was also depicted in energy metabolism, with 27 genes being down-regulated in contrast to 3 that were up-regulated (Figure 4c).

Phagotrophy

To study the potential for phagotrophy in monocultures of P. parvum we examined the expression pattern of genes involved in endocytosis, phagosome, lysosome and peroxisome formation (Figure 5). Genes involved in endocytosis were mainly up-regulated under low salinity and P starvation (Figure 5). More specifically, ten genes were up-regulated under either low salinity or P starvation. Concerning cell density, two genes were differentially expressed and all up-regulated in low cell density. The expression of 16 genes involved in lysosomes, resulted in two distinct clusters, with eight genes presenting a higher expression under either low salinity or phosphorus depletion (Figure 5). No clear patterns in relationship to the treatments could be observed for genes involved in the phagosome or peroxisome formation.

Toxin production

The expression dynamics of PKSs were examined in order to check their expression profile under the different treatments. The expressed PKSs formed 10 distinct clusters (clusters A-J) (Figure 6). The expression of PKSs showed higher expression values under the high cell density (bloom conditions) conditions, with the exception of the low salinity and P-replete condition. Clusters D, E, I and J showed the higher over-all expression values, and mainly consisted of modular Type I PKSs. The length and functional domain organization of contigs belonging to the clusters D and I are provided in Table 2. The expression of contig UIO223_409_c0_g1_i1 was consistently high across all treatments, but the highest expression value was observed under P-starvation and independently of salinity. This contig is characterized by the presence of a polyketide-type polyunsaturated fatty acid synthase (pfaA), which are involved in the biosynthesis of omega-3 polyunsaturated fatty acids. The effect of P starvation and higher expression of PKSs under this condition was more visible in clusters F and H, and especially for the high salinity treatment. In cluster D, the expression of three contigs (UIO223_7525_c0_g1_i2, UIO223_34765_c0_g1_i1 and UIO40_c0_g1_i4) was higher under P starvation and low cell density.

DISCUSSION

Cell growth and physiological responses to phosphorus depletion and low salinity

P. parvum is an euryhaline species with the ability to grow in a wide range of salinities from 0.5 to 45 (Edvardsen and Paasche, 1998; Beszteri et al. , 2012; Granéli et al. , 2012). Much research has been conducted on the relationship between salinity and blooming potential in P. parvum . In the field, blooms of P. parvum usually occur in estuaries with low salinity (2-8), while in some cases salinity can be even <2 (Richardson and Patiño, 2021). In the lab, growth rates under different salinities are variable and no clear relationship can be observed between these two parameters (Roelke et al. , 2011; Patiño et al. , 2014; Rashel and Patiño, 2017), thus, highlighting the importance of culturing conditions and the different responses of P. parvum isolates. In the current study, growth rates for P. parvum growing under low salinity presented statistically significant higher growth rates.

Cellular particulate organic carbon and nitrogen contents were negatively influenced by cell density, irrespective of salinity and inorganic phosphorus concentration. The POC per cell measurements for the P deplete cells were similar to the quota of 24.2 pg cell^-1 previously reported for P depleted P. parvum (Lundgren et al. , 2016). It is a possibility that the general decrease in POC content in high cell density cultures (not nutrient limited) could be attributable to lower dissolved inorganic carbon (DIC) availability due to removal of DIC at high cell density. Even though DIC was not measured in this study, the pH of the cultures was controlled by bubbling and we have no indications that cells were carbon limited. Also, the decrease in cellular POC was accompanied by a decrease in cellular PON, suggesting a general decrease in cell size. These findings support the concept that cellular stoichiometry and growth rates are considered to be strictly related, with decrease in carbon content when cells present lower growth rates (Garcia et al. , 2016; Moreno and Martiny, 2018)

Genes involved in carbohydrate metabolism seemed to be mainly affected by salinity and cell density, suggesting a general re-adjustment of carbon-related pathways under these conditions. The expression patterns of carbohydrate metabolic genes suggest increased carbon turnover related to the growth rates of cells under the different salinity and cell density conditions. To fulfill the energy and carbon units demands key enzymes involved in glycolysis and pyruvate metabolism were significantly up-regulated. Furthermore, the glyoxylate cycle, a way of bypassing the two oxidative decarboxylation reactions of the TCA cycle and directly converts isocitrate through isocitrate lyase and malate synthase into malate and succinate, showed a decrease in gene expression. Carbohydrate rearrangement seems to be a critical factor under different salinity treatments, as it is also related to the production of osmolytes for the balance of osmotic pressure in cells (Harding et al. , 2016, 2017).

Prymesin production, cell lysis and phagotrophy

Prymnesins are currently considered to be the toxins in P. parvumresponsible for the lytic effects observed on other organisms (Rasmussenet al. , 2016; Binzer et al. , 2019). Recent advances in prymnesin quantification have facilitated our understanding of toxin dynamics, production and release in the environment. Svenssen et al. , 2019 developed an indirect method for the estimation of prymnesin type B and found that the majority of the toxin is intracellular. Moreover, the combination of bioassays and biochemical approaches have reinforced the idea that prymnesins are the actual factors behind acute toxicity.

In the present study, the intracellular prymnesin toxin content (cell quota) was significantly higher in high cell density cultures, suggesting cell density or growth rate (accumulative) dependent mechanisms that define the cellular toxin content. This apparent stimulation of prymnesin production as the cultures grow dense has never been observed before in P. parvum or any other algae producing lytic compounds, and this discovery needs further scrutiny. Increased production and release of lytic toxins in dense populations will be favorable for the alga, because at low P. parvum cell concentrations there will be not benefits of such compounds, given their concentration dependent mode of action (Tillmann, 2003).

In the present study, prey mortality was significantly higher at low salinity, irrespectively of the nutrient conditions. This observation did not match with the measurements of prymnesin production. The experimental set-up and results do, however, not allow us to distinguish between potential higher toxicity due to increased release of prymnesins in the medium at low salinity, or higher prey sensitivity due to low salinity effects. Additionally, the toxin profile was also influenced by salinity and P, as the relative composition of prymnesin analogues was affected. Consequently, we cannot rule out, that the different analogues could perform various levels of lytic activity or that other unknown compounds may be involved in prey cell lysis.

The highest feeding rates were observed in P-deplete treatments, while the lowest feeding rates were observed at high salinity and P-replete treatments. At low salinity, the rates of prey mortality were initially very high. This could have led to low feeding rates because the prey cells had lysed before P. parvum had the opportunity to ingest them, but they may have benefited from osmotrophy of newly released or produced DOM of the lysed cells. Nevertheless, neither phagotrophy nor uptake of re-leased DOM did seem to boost the growth rate of P. parvum . However, the prey was lysed, before the P. parvum got nutrient limited in any of the treatments in the present study, and our data supports the idea that phagotrophy in P. parvum may be a strategy to compensate for nutrient limitation. Our findings indicate that low salinity is the main factor that induces higher degree of prey lysis, while P starvation seem to induce higher phagotrophy rates inP. parvum .

Phagotrophy involves the internalization of food via endocytosis and the subsequent formation of a phagosome, which fuses with lysosomes for further digestion. In mixotrophic Alexandrium , expression of endocytosis genes was more enhanced in a lytic strain as compared to a non-lytic strain, highlighting the potential coupling be-tween the processes of lytic toxin production and phagotrophy (Wohlrab et al. , 2016). In the current study, gene expression associated with phagotrophy-related processes were induced by both salinity and P depletion, indicating that mixotrophy is induced in P. parvumunder stress conditions. Genes involved in endocytosis were consistently up-regulated under low salinity and P starvation, with both treatments inducing strong prey lysis/mortality and feeding in P. parvum . Expression of phagotrophy-related genes has been previously related to ingestion rates in mixotrophic flagellates (McKie-Krisberg et al. , 2018). Moreover, factors that influence the expression of such genes include prey availability and light (Lie et al. , 2017). Our incubation experiments showed higher phagotrophy under P depletion independent of salinity. However, potential phagotrophy and/or DOM uptake might have been underestimated, due to higher excretion of toxins at low salinity, which could cause prey lyses.

Polyketide synthases

Polyketide synthase genes are secondary metabolites that have been extensive studied in marine protists (John et al. , 2008; Monroe and Dolah, 2008; Kohli et al. , 2016). Many phycotoxins are polyketides suggest that PKSs are involved in toxin biosynthesis. Prymnesins can have carbon backbone lengths, from 83 to 91 C (Binzeret al. , 2019). The considerable carbon length of toxins highlights the complexity of the underlying molecular biosynthetic mechanisms. Their biosynthesis could potentially involve both PKSs and perhaps partly fatty acid synthases, as has previously been suggested for the biosynthesis of antibiotics from bacteria (Masschelein et al. , 2013). Light-dark dependent expression of PKS has been previously studied in synchronized cultures of the ichthyotoxic haptophytePrymnesium polylepis , which led to the identification of 13 potential PKS sequences (John et al. , 2010). In P. parvum , physiological treatments such as high irradiance and low salinity shock have been shown to induce higher copy numbers of PKS (Freitag et al. , 2011). Polyketide synthase genes have previously been reported inP. parvum strains that produce different types of prymnesins (Anestis et al. , 2021), however, the development of efficient genetic manipulation techniques would be required for better understanding of the connection between gene and product. In the current study, the expression of PKSs was strongly influenced by cell density, coinciding with higher cellular content of prymnesins. The expression patterns highlight the presence of PKS gene clusters, which were, interestingly, higher expressed under P depletion. Under these treatments, phagotrophy was also higher, thus argue for the potential involvement of their product in prey capture or lysis. These finding certainly need further investigation.

Implications for bloom maintenance

Maintaining high cell concentrations over weeks during algal blooms challenges cellular and metabolic performances of algal cells as nutrients become limiting. Our findings suggest rearrangements of carbon metabolism and down-regulation of energy anabolism, with cellular transcription and translation, being considerably inhibited during P depletion. However, even under P depletion the expression of polyketide synthases was high and increased at high cell densities and under P depletion. Assuming that polyketide synthases are indeed the encoding genes for prymnesins and given the high amount of intracellular prymnesin, this would imply that cells under P depletion invest substantial amounts of carbon units and energy in the production of these energetically expensive secondary metabolites. The benefit of such an investment could be higher competitive abilities via prey capture, lyses and/or increased feeding. Such a response would be evolutionary advantageous when prey is available, and the production and release of lytic toxins could contribute to increased nutrient availability and/or direct prey uptake by P. parvum . The additional up-regulation of endocytic genes under P depletion and low salinity could suggest the preparation of the molecular toolkit for phagocytosis when nutrients become the limiting factor. Such a response would allow the further domination of P. parvum under bloom condition and the engagement in alternative nutrition modes for the maintenance of its population.

Conclusions

Salinity, phosphorus availability and cell density had important impact on different aspects of P. parvum physiology. Prymnesin contents were significantly higher in high cell density cultures, while salinity and P availability influenced the composition of prymnesin analogues. Gene expression of metabolic genes was in accordance with the observed growth rates, with up-regulation of catabolic genes under high growth rates. Prey lysis was faster at low salinity, while P starvation induced higher phago-trophic rates. Phagotrophy did not increase cell growth, supporting the idea that mixotrophy in P. parvum may be a survival strategy to cope with nutrient limitation, particularly expected during bloom conditions. The gene expression data supports the increased potential for mixotrophy under low salinity and P-starvation conditions, with up-regulation of endocytosis related genes. The expression of PKSs was generally higher at high cellular toxin content condition, and especially at P depleted cells, which were characterized by increased phagotrophy.

AUTHOR CONTRIBUTIONS

KA and UJ designed the study, KA performed the experiments, KA, SW, EV and UJ analyzed and interpreted the data, UJ supervised the study. All authors contributed to drafting the manuscript.

ACKNOWLEDGEMENTS

This work has been funded by European Union’s Horizon 2020 research and innovation programme under grant agreement No. 766327. UJ and SW were financially supported under the POF IV Topic 6 ST 6.2 of the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Germany and the Helmholtz Institute for Functional Marine Biodiversity. EV received funding of the Austrian Science Fund (FWF) through an Erwin-Schrödinger fellowship [J3895-N28] and by the University of Vienna.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The raw read sequences and transcriptome have been deposited at the National Center for Biotechnology Information under the BioProject PRJNA718746. The data for the physiological parameters is provided at the Appendix.

TABLE AND FIGURE LEGENDS

Table 1 | Regulation of the expression of catabolic and anabolic genes in Prymnesium parvum as fold-change in the comparison between low salinity and high salinity.

Table 2 | Contigs encoding for highly expressed polyketide synthases in Prymnesium parvum (clusters D and I of Figure 6), their length in base pairs and their corresponding functional domain organization.

Figure 1 | Overview of sample parametersPrymnesium parvum for different cell densities and treatments (HS = high salinity, LS = low salinity, HP = high phosphorus and LP = low phosphorus). Cell growth (a) and indicated with arrows are the sampling time points for low cell density (LD) and high cell density (HD). Further parameters include phosphate concentration (b); particulate organic carbon (c); particulate organic nitrogen (d); carbon:nitrogen ratio (e); chlorophyll-a concentration (f).

Figure 2 | Total cellular prymnesin content inPrymnesium parvum (attomol cell^-1) (a) and percentage prymnesin analog composition (b) for all treatments (HS = high salinity, LS = low salinity, HP = high phosphorus and LP = low phosphorus) and cell densities (LD = low density and HD = high density). Concentration of Teleaulax acuta (cells mL^-1), when incubated with Prymnesium parvum (c), and percentage ofP. parvum cells containing a food vacuole (d).

Figure 3 | Overview of statistics for the generated transcriptome of Prymnesium parvum (a), Principal Components Analysis (PCA) of the transcriptomes for the samples and their corresponding cell densities and conditions (b), and number of transcripts for which salinity (Sal), phosphorus availability (P), cell density (CD) and all combinatory interactions had a statistically significant (adjusted p < 0.01) effect on their expression (c).

Figure 4 | Number of up- and down-regulated genes inPrymnesium parvum , assigned per KEGG category for cellular processes (orange), environmental information processing (green), genetic information processing (blue) and metabolism (brown). The comparisons include (a) low phosphorus (LP) vs high phosphorus (LP), (b) low salinity (LS) vs high salinity (HS) and (c) low cell density (LD) vs high cell density (HD).

Figure 5 | Expression level of differentially expressed genes involved in transport and catabolism in Prymnesium parvum . The categories include genes involved in endocytosis, phagosome, lysosome and peroxisome. The heatmap shows the expression levels across different cell densities (low and high) and for all treatments (HS = high salinity, LS = low salinity, +P = high phosphorus and -P = low phosphorus).

Figure 6 | Expression level of polyketide synthase genes in Prymnesium parvum . The heatmap shows the expression value as transcripts per million (TPM) across different cell densities (low and high) and for all treatments (HS = high salinity, LS = low salinity, +P = high phosphorus and -P = low phosphorus). Hierarchical clustering was performed based on the expression profile and lead to the formation of 10 clusters (A-J).