Variation in prostaglandin metabolism during growth of the diatom Thalassiosira rotula

Prostaglandins (PGs) are hormone-like mediators in many physiological and pathological processes that are present in all vertebrates, in some terrestrial and aquatic invertebrates, and have also been identified in some macroalgae. They have recently been reported also in marine microalgae but their role as chemical mediators is largely unknown. Here we studied the expression pattern of the PG biosynthetic pathway during different growth phases of the centric diatom Thalassiosira rotula and assessed the release of PGs in the surrounding environment for the first time. We show that enzymes responsible for PGs formation such as cyclooxygenase, prostaglandin E synthase 2-like and prostaglandin F synthase are mainly expressed at the end of the exponential phase and that PGs are released especially during the stationary and senescent phases, suggesting a possible signaling function for these compounds. Phylogenetic analysis of the limiting enzyme, COX, indicate the presence in diatoms of more than one enzyme related to the oxidative metabolism of fatty acids belonging to the peroxidase-cyclooxygenase superfamily. These findings suggest a more complex evolution and diversity of metabolic pathways leading to the synthesis of lipid mediators in diatoms.

Arachidonic acid (ARA), eicosapentaenoic acid (EPA), eicosatrienoic acid (ETrA) and docosahexaenoic acid (DHA) are polyunsaturated fatty acids that are physiologically important for animals at all taxonomic levels, including humans. EPA and DHA contribute to the healthy functioning of the cardiovascular system 1 and are precursors of important classes of fatty acid-derivatives playing multiple signaling roles in inflammation and immune responses, platelet aggregation and tumor growth 2,3 . Among these, the inflammation process is one of the most important mechanisms adopted by organisms in response to various external stimuli 3 . Inflammation involves a complex interplay of signaling molecules whose final goal is to restore the healthy status of a cell or tissue. Consequences of sustained inflammation are indeed the development of serious diseases such as cancer and autoimmune disorders 4 .
Included in the eicosanoids are prostaglandins (PGs) synthesized principally from arachidonic acid (ARA) in animals, but also from EPA and ETrA, through the enzymatic route initiated by cyclooxygenase (COXs) enzymes 5 . PGs are molecules with a hormone-like behavior playing a prominent role in many physiological processes that have been principally studied in animals 3 . The expression of the COXs enzymes is mandatory for their synthesis. COXs exist in two isoforms that differ for their subcellular localization and for their expression timing. COX-1 is located in the endoplasmic reticulum and is constitutively expressed at constant levels in many tissues unless external cues, such as tumor promoting factor, cytokine and growth factor, induce an increase in its expression level. COX-2 is the inducible form, which is not detectable unless a trigger similar to those that stimulate COX-1 expression occurs. COX-2 is located in the nuclear envelope and appears to be a target for cancer therapy 6 .
PGs synthesis, in mammals, is initiated by phospholipases (PLAs), a family of enzymes that hydrolyze membrane phospholipids liberating the precursors, ARA, EPA, and ETrA. These are rapidly converted, through cyclization and inclusion of molecular oxygen, into the unstable metabolite PGG 2 by the action of COXs enzymes that subsequently reduce it into PGH 2 . The PGH 2 is then transformed into the ultimate prostaglandin E 2 , D 2 , F 2α , prostacyclin or tromboxanes by the successive action of PGE, PGD, PGF, prostacyclin and tromboxanes A

Results
In-silico structural reconstruction of TrotCOX. Since COX is the most relevant enzyme in the pathway, essential for the initiation of PGs synthesis, we used Phyre 2 (Protein Homology/analogY Recognition Engine V 2.0) 24 to predict the possible structure of TrotCOX protein and to confirm that it belongs to peroxidase-cyclooxygenase superfamily. The multiple alignments developed by the Phyre program retrieved as best hits 13 structures of proteins belonging to the peroxidase-cyclooxygenase family ( Table 1, Supplementary Fig. 1). All of them shared 100% confidence and 94-99% coverage ( Table 1). The first structure in the list shared a 19% identity and corresponded to a Bos taurus lactoperoxidase (BtLPO) already reported as representative of the peroxidase-cyclooxygenase superfamily 12 (Table 1, Supplementary Fig. 2). Interestingly the second to the fifth structures corresponded to an Ovis aries prostaglandin G/H synthase (OaCOX1) all sharing a 28% identity, confirming the cyclooxygenase structure of the TrotCOX protein (Table 1, Supplementary Fig. 3). Structures from line 6 to 13, sharing lower identity (26% to 19%, Table 1), were however representative of cyclooxygenases, myeloperoxidase or oxidoreductase protein structures from Homo sapiens, Mus musculus, Arabidopsis thaliana and Oryza sativa, all of which are representative of the peroxidase-cyclooxygenase superfamily. Figure 1b,c show the results of the comparison among TrotCOX (Fig. 1a) predicted structure and known models of the peroxidase-cyclooxygenase family, alias the BtLPO (3BXI.pdb, 10.2210/pdb3bxi/pdb) 12,25 and OaCOX1 structures. The main alpha helices that characterize the cyclooxygenases and accommodate the heme prosthetic group are present and almost perfectly overlapped with the references BtLPO and OaCOX1.
The in silico predicted model constructed on the BtLPO showed conserved domain sites in the diatom sequences ( Fig. 1d and Fig. 2): catalytic arginine was found (R275) in the distal heme side, embedded in the characteristic motif XRXXEX, while glutamate that is usually conserved, was not in TrotCOX but was mutated to leucine (L278). The latter seems to be a characteristic of diatoms as none of the diatom sequences ( Supplementary  Fig. 4) here analyzed showed a glutamate in that position. The proximal heme sides are relatively conserved in TrotCOX with the highly conserved arginine (R371), histidine (H373), isoleucine (I376) and glutamic acid (E454) unvaried in TrotCOX and in most of the other sequences. Surprisingly, the conserved asparagine and arginine that are normally bonded via a hydrogen bond were substituted with a threonine (T453) and a cysteine (C456), respectively. The calcium binding site in TrotCOX shows the relatively conserved Isoleucine (I184) and LYG motif (L198, Y199, G200) (Figs. 1d and Fig. 2). The sequence representation by Logo (Fig. 2) confirmed this residue conservation also among all the diatom sequences considered for the phylogenetic analysis.
Using a previous alignment 15 , populated with other sequences retrieved from NCBI via a pBLAST, we performed a refined phylogenetic analysis comparing the COX sequence from T. rotula with that of other organisms (Supplementary Tables 1 and 2). Most of the diatom sequences clustered in a large well-supported clade (Maximum Likelihood (ML) bootstrap value 100) closely related to the animal COX clade (family 4, cyclooxygenases), although the phylogenetic relationship was not well resolved (ML 35) (Fig. 3).
This large diatom clade was divided into a main well-supported clade including the vast majority of diatom sequences, among which the TrotCOX sequence (Fig. 3), and in an end-clade (ML 94) grouping both prokaryotic (Nostoc sp., a cyanobacterium, the actinomycetales Herbidospora mongoliensis and Rhodococcus gordoniae) and diatom (Grammatophora oceanica, Staurosira sp., and Chaetoceros cf. neogracile) sequences. This diatom/ bacterial clade can be hypothesized as being the 5a family within the peroxidase-cyclooxygenase super family 12 . Noteworthy, the pennate biraphid diatom sequence from Nitzschia sp. clusters basally to the animal COX clade (ML 56). In addition to that large clade, a second very well resolved (ML 100) super clade was composed of two clades. One of these forks includes only diatom sequences (Fistulifera solaris and S. marinoi 31509 and 1511); whereas the other bifurcates in two sister clades. One of these sister clades groups only two prokaryotic sequences, and the other bifurcates in two branches. The first branch groups plant α-DOX proteins (family 4 cyclooxygenase 12 ) and the other groups diatoms (Skeletonema gretae, Pseudo-nitzschia pungens and Fragilariopsis kerguelensis) (Fig. 3). Figure 4 illustrates differential expression analysis by qPCR of the three genes involved in the PG pathway during T. rotula growth (Fig. 4d). One-way ANOVA analysis shows, for each gene, a statistically significant difference (0.0001 < p-value < 0.0002) among the six time points considered along the 10-day T. rotula growth curve (Fig. 4, Table 2).

Gene expression analysis.
Each gene showed a statistically significant expression peak around the late exponential/early stationary phase, i.e. day 5 for COX and at day 4 for mPTGES and PTGFS ( Fig. 4a and Table 2). In addition, COX-day 6 expression persisted with significant difference versus COX-day 7 and COX-day 10 expressions. mPTGES-day 3 expression was statistically significantly higher than mPTGES-days 6, 7 and 10 expressions ( Fig. 4b and Table 2) but had levels lower than mPTGES-day 4 expression. Similarly, PTGFS-day 3 expression was statistically significantly lower than day 4 but higher than PTGFS-days 7 and 10 expressions. Analogously, PTGFS-day 5 expression was still statistically significantly higher than PTGFS-day 7 and 10 expressions ( Fig. 4c and Table 2). The expression peaks were only transient because a quick decrease in the expression levels occurred soon afterwards, reaching very low levels at the end of the stationary (day 7) and senescent (day 10) phases. Moreover, the initial basal expression levels of the three enzymes were completely different. COX expression was much higher compared to mPTGES and PTGFS, and mPTGES was more expressed compared to PTGFS (Fig. 4, Day3). However, from day 7 to 10, the expression of all three genes decreased, reaching similar levels. The reference genes 23 utilized were stable, indicating that the down regulation of the analyzed genes was not due to senescence but rather to a specific pathway response.

LC-MS/MS analysis.
In addition to gene expression levels, we also studied the pathway activity by identifying the PGs released outside the cells by T. rotula and measured their concentration in the culture medium ( Fig. 5a-c). By LC-MS/MS in comparison with pure standards (Supplementary Figs. 5-8), we detected three PGs in the culture medium, namely PGE 2 , PGB 2 and 15-d-PGJ 2 . In addition, in the same medium, we also identified PGA 2 , PGD 1 , PGD 3 , PGE 1 , PGE 2 , PGE 3 , PGEM, PGFM, 15-d-PGD 2 and 2,3-dinor-11b-PGF 2 that we classified as "putative" since they have been identified through comparison of experimental mass fragmentation data obtained in this study with those reported in the literature (Supplementary Fig. 9, Supplementary Table 3)   www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 2. TrotCOX vs BtLPO amino acid sequence alignment at the distal and proximal heme sides and at the calcium binding sites 32 . The sequence Logo of the same portion of the alignment for the diatom specific clade is shown below. The sequence Logo is a resumed representation of an alignment among two or more sequences. The height of each letter is proportional to its frequency and the most common one is placed on the top 55 . Important amino acidic residues reported in the text are indicated on the T. rotula sequence. BtLPO = Bos taurus LPO; T. rot = Thalassiosira rotula. Figure 3. Phylogenetic analysis of the Thalassiosira rotula COX. Diatom sequences are scattered in the ML phylogenetic tree, either in diatom specific clades, or clustered together with sequences from other taxa. Legend: TrotCOX is indicated in red; cyan dot: diatom sequences not clustering the diatoms-specific clade; cyan-shaded clade: diatom specific clade; 5a: 5a-family in the diatom/bacterial clade; 4 COX: animal COX clade (family 4, cyclooxygenases); 4 α-DOX: plant α-DOX proteins (family 4, cyclooxygenases); 5b: 5b-family cyanobacteria. The nomenclature used for the family classification follows Zámocký et al. 12 .
Overall, a high variability among replicates was evident as indicated by the large error bars (Fig. 5a-c). In addition, we observed an up-and-down trend of each measured PG concentration from days 3 to 5. Conversely, measurements from day 6, corresponding to the full stationary phase, to day 10, showed a linear trend toward increasing concentration values, except for the PGE 2 , PGE 3 , 15-d-PGJ 2 , and PGFM (Fig. 5a-c). In particular, PGE 2 had a constant low level reaching zero on day 10 (Fig. 5a), while PGEM showed an exponential trend starting from day 4, i.e. late exponential phase of growth, reaching the highest value among all the other PGs identified (Fig. 5b).
One-way ANOVA analysis (Table 3) shows a statistically significant difference over the considered time points for PGA 2 , PGB 2 , PGD 1 , PGD 3 , PGE 1 , PGEM, 15-d-PGD 2 and 2,3-dinor-11b-PGF 2 . Comparison of day 4 versus day 10 values were statistically significant for the majority of the compounds, due to concentrations dropping to zero on day 4 (Fig. 5). PGEM concentrations, varied significantly from day 4 to 10, even if some differences were detected, i.e. Day5 vs Day6, Day5 vs Day7 and Day6 vs Day7 were not significant because of the large error bars (Fig. 5c, Table 3).
One-way ANOVA, calculated considering day 6 as reference point, shows the same significant differences over time points confirming that differences in concentrations at day 6 vs day 10 are statistically significant for all the PGs, except for PGE 2 , PGE 3 , 15-d-PGJ 2 , and PGFM (Table 4). www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Diatoms are a very important group of microalgae populating all aquatic niches and able to fix about 20% of global carbon production 28 . Their elevated adaptability is due to a rich set of metabolic pathways coded by their genome 29 . Their crucial ecological role 30,31 is coupled to an emerging biotechnological interest due to their ability to produce high added value molecules. One of the metabolisms recently identified in these microalgae is the pathway for the enzymatically-mediated synthesis of prostaglandins (PGs), principally studied in mammals but less in plants or microorganisms 15,23 .
The key step in PG biosynthesis is catalyzed by COX enzymes belonging to the heme-peroxidase protein superfamily. Evolutionary studies of this superfamily demonstrated an independent evolution of four superfamilies, including peroxidase-cyclooxygenase 12 , each possessing a peculiar folding of the heme peroxidase domain. In particular, the peroxidase-cyclooxygenase superfamily has a wide distribution in all living kingdoms 12 .
Phylogenetic analysis of the Thalassiosira rotula COX presented here widens the analysis reported previously in Di Dato et al. 15 with the addition of a few sequences that have better resolved some of the clades identified in that study. In particular, the sequences from Grammatophora oceanica 5553, Staurosira 23824, Chaetoceros cf. neogracile 1470 clustered with Nostoc WP015113127 thereby questioning the validity of the analysis. It was not clear enough whether the clade grouping the three diatom and Nostoc sequences could be identified as a diatom specific clade derived from cyanobacteria. In the present analysis, new bacterial sequences were added and clustered together with the above-mentioned diatom and Nostoc sequences showing that this clade is indeed a bacterial-like clade. Moreover, our phylogenetic analysis highlights the existence of more than one peroxidase-cyclooxygenase enzyme related to the oxidative metabolism of fatty acids in diatoms. The four sub-families in the peroxidase-cyclooxygenase superfamily (family 4: cyclooxygenase and α-DOX; 5a and 5b, bacterial peroxidase-cyclooxygenase families) did not evolved from one single protein, but rather appeared independently during evolution 32 . COX, were formerly known as 'animal heme-dependent peroxidases' . However, as demonstrated here and in our previous work 15,23 , this enzyme is present and active in diatoms as well, thus the denomination 'animal heme-dependent peroxidases' would deserve a revision in order to avoid confusion.
α-DOX belongs to the same peroxidase-cyclooxygenase superfamily as COX and some diatoms (e.g. S. marinoi) possess this protein along with COX. This questions the origin and evolution of this class of enzymes in diatoms and shows that this group of microalgae is capable of manipulating fatty acids in different ways: diatom α-DOX enzymes may act on medium chain saturated fatty acids, such as 16:0 33 , while COX on long chain-PUFAs (e.g. AA and EPA) 34 . Also, the close proximity of bacterial and diatom sequences and the topology of the tree we present, questions whether diatoms acquired the genes coding for COX via lateral gene transfer 35 , which occurred during diatom evolution [36][37][38] or whether these genes were acquired during endosymbiotic events, at least for species like S. marinoi possessing both α-DOX-like and COX-like enzymes. In the case of T. rotula, α-DOX-like is absent, while the pennate model diatom Phaeodactylum tricornutum completely lacks the enzymatic pathway leading to prostaglandins though it can produce isoprostanoids non-enzymatically 22 . Indeed, the ability to synthesize eicosanoids, including lipoxygenase-derived oxylipins and aldehydes, is not robustly conserved in diatoms 39,40 where genes involved in these metabolisms might have been lost during evolution.  www.nature.com/scientificreports www.nature.com/scientificreports/ In our in silico study of TrotCOX protein sequence and structure we found an overall homology and the conservation of motifs characterizing the heme-peroxidase protein superfamily, corroborating also the assignment of this sequence to the peroxidase-cyclooxygenase superfamily 12 . The catalytic sites, the motifs and their order along the protein sequence, are very conserved compared to the representative sequences of the peroxidase-cyclooxygenase superfamily, even if with some differences. The glutamate embedded in the characteristic motif XRXXEX, usually conserved, in TrotCOX is mutated into a leucine. This mutation seems to be characteristic of diatoms that not only show mutations in TrotCOX but also in other diatom sequences considered for the phylogenetic analysis. Considering the different cellular environments, the leucine residue may still play in diatoms the same role as glutamic acid in other organisms although the glutamic acid is a charged amino acid while leucine is a hydrophobic uncharged amino acid. In addition, in the proximal heme sides, the conserved asparagine and arginine normally linked via a hydrogen bond, playing crucial roles in alpha helix structures 41 , are substituted with a threonine and a cysteine, respectively. This change should not compromise the alpha helix stability present in this region.
Our experimental data demonstrate that TrotCOX is active and works like the animal COX, being able to synthesize the same PGs produced by animals. The presence of PGE 3 , a specific derivative of EPA, is in line with the abundance of this PUFA precursor, which represents a hallmark of diatoms 42 .
Interestingly, the general increase of PGs released outside the cells started from day 5 onwards, i.e., when cells were entering the stationary phase that coincides with the maximum expression of COX. This stage is generally associated with an increase in the naturally associated bacterial population. During the transition between the stationary and senescent phases, nutrients decrease while dying algal cells generate organic matter that can feed www.nature.com/scientificreports www.nature.com/scientificreports/ associated bacteria. Bacteria are able to assimilate phosphate better than algae, especially at low concentrations 43 , and compete with algae for inorganic nutrients 44 . This increased pressure may stimulate PGs synthesis in diatoms as already shown in animals, where COX expression is stimulated in the presence of bacteria 45 and in humans where PGE 2 have been correlated to viral load and infection severity in influenza [46][47][48][49] . In line with this hypothesis, similar compounds, such as hydroxylated eicosapentaenoic acid 15-HEPE, were shown to be up-regulated and released from cells when the diatom S. marinoi was exposed to pathogenic bacteria 50 . PGE 2 is the most abundant prostanoid in the human body, but it is also very unstable since it is rapidly converted into PGE-M 46-49 that is considered to mirror the systemic levels of PGE 2 formation. In our study, PGE 2 was present at low levels in the medium during the exponential growth phase and was almost absent in the late stationary phase (from day 7) while PGE-M levels increased exponentially up to the senescent phase (Fig. 4). The peak in the expression of the related enzyme, PTGES, at the onset of the stationary phase and the sustained release of PGs afterwards suggest a possible role of prostaglandins in communication and cell signaling. Their release outside the cells, in addition to their sustained presence in a saline environment, is quite striking and has never been observed before in diatoms. PGs are known to exert a wide range of effects in different organisms, including the induction of inflammatory processes, injury and pain in humans where they have been best studied. The fact that they have been identified in organisms ranging from unicellular diatoms, corals and jellyfish to arthropods, mollusks and mammals denotes that they regulate important physiological processes that have been conserved through evolution.    Table 4. One-way ANOVA analysis over the time points day 6 to day 10 of the ten day Thalassiosira rotula growth curve, considering Day 6 as reference time point.
In conclusion, this study confirms that diatoms possess a molecular toolbox generally believed to be unique to higher organisms such as mammals. The production and release of PGs by some diatoms and the variation in the expression levels of PG biosynthetic pathway during different growth stages, strongly suggest a relevant and possible signaling role for these molecules within the plankton community that needs to be further investigated. This characteristic may have contributed to render diatoms one of the most successful groups of organisms in the world's oceans.

Methods
Strain and cell cultures. Thalassiosira rotula, strain FE80C, was isolated in 2011 in the Gulf of Naples (40°48.5'N, 14°15'E), Mediterranean Sea. Clonal cultures were established by isolating single cells or short chains from phytoplankton net samples collected from the surface layer of the water column. Cultures were grown in sterile filtered oligotrophic seawater at 36 psu amended with f/2 51 nutrients and maintained at a temperature of 20 °C, at 12:12 h light:dark cycle, and with a photon flux of 100 μmol photons m-2 s −1 .
10-liter cultures of T. rotula, in triplicate, were used to follow their growth from day 3 to day 10. Every day, 250 mL of each culture was harvested by filtration onto 1.2 μm pore size filters (RAWP04700 Millipore) and immediately frozen in liquid nitrogen. 100-200 mL of culture media recovered from the cell filtration was collected and stored at −80 °C until sample processing. Initial cell concentrations were about 5000 cells/mL upon inoculation. Culture growth was monitored daily from samples fixed with one drop of Lugol (final concentration of about 2%) and counted in a Bürker counting chamber under an Axioskop 2 microscope (20×) (Carl Zeiss GmbH, Jena, Germany).

RNA extraction and reverse transcription.
To proceed to total RNA extraction, filters were covered with 1.5 mL TRIsure (BIO-38033, Bioline) to which glass beads (G1277, Sigma-Aldrich) were added. Cells were disrupted on a thermo-shaker (Eppendorf) at 60 °C for 10 minutes at 1200 r.p.m. Filter and glass beads were discarded, and the extraction was continued according to the TRIsure instructions. DNase treatment was carried out using DNase I recombinant, RNase-free (Roche, Basel, Switzerland) according to manufacturer's protocol to eliminate potential genomic DNA contamination. The efficiency of the DNase digestion was checked by testing DNA primer capability to amplify an amplicon with traditional PCR. PCR reactions were carried out in 25 μL volume with 2,5 μL of 10× PCR reaction buffer (Roche, Basel, Switzerland), 2,5 μL of 2 mM dNTP, 0.3 μL of 5 U/ μL Taq (Roche, Basel, Switzerland), 1 μLof 20 pmol/μL for each oligo, 1 μL of RNA templates and nuclease-free water up to 25 μL. The PCR program consisted in a denaturation step at 95 °C for 5 min, 35 cycles at 95 °C for 30 s, 53 °C for 30 seconds and 72 °C for 45 s, and a final extension step at 72 °C for 7 min. Amplified PCR products were analyzed by agarose gel electrophoresis.
Total RNAs were purified and concentrated using RNeasy MinElute Cleanup Kit (Qiagen, Venlo, Netherlands) and eluted in 30 μL of RNase-free water. Concentrations of the resulting RNA samples were assessed by absorbance at 260 nm (ND-1000 Spectrophotometer; NanoDrop Technologies, Wilmington, DE, USA). The integrity of total RNA was checked by agarose gel electrophoresis. One microgram of each RNA sample was retro-transcribed in complementary DNA (cDNA) following the manufacturer's instructions (5X All-In_One RT Master Mix, abm, Applied Biological Materials Inc., Richmond, Canada) using the T100 Thermal cycler (Bio-Rad Laboratories, Hercules, CA, USA). Retro-transcribed samples were again checked for DNA contamination by traditional PCR, as above, using intron-spanning primer pairs. Phylogenetic analyses. For cyclooxygenase phylogeny, the alignment from Di Dato et al. 15 was retrieved and all the sequences were blasted against the NCBI database. The best hits from each pBLAST were added to the analysis. Sequences were visualized and aligned in BioEdit Sequence Alignment Editor 52 using the ClustalW alignment algorithm implemented in BioEdit. The alignment was performed using MUSCLE online software and compared to the ClustalW output. The two alignments were identical. For maximum likelihood (ML) phylogenetic analysis both MEGAX 53 Supplementary Table 1. The corresponding sequences are reported at the end of the supplementary file, both the aligned sequences (with gaps) and the ungapped sequences. The sequences from strain SkelB that were annotated as Skeletonema dohrnii in MMETSP database were changed to S. marinoi after ribosomal sequence identification. Sequence ID were left unchanged. The phylogenetic tree was built on scale in the branch length (scale bar reported on Fig. 3). Sequences are identified by the species name, followed by their MMETSP ID (deprived of the taxon).
In-silico TrotCOX protein structure modelling and comparison. We used Phyre 2 (Protein Homology/ analogY Recognition Engine V 2.0) 24 to predict the possible structure of TrotCOX protein. In order to compare the TrotCOX predicted structure to known models of the peroxidase-cyclooxygenase family, the BtLPO (3BXI. pdb, 10.2210/pdb3bxi/pdb) 12,25 and OaCOX1 structures were retrieved from the RCSB database (https://www. rcsb.org) and compared to TrotCOX using PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC.
Bioinformatic identification of the prostaglandin pathway. Pathway annotated as 'prostaglandin biosynthesis' was found among the second level pathways list generated within the Annocript pipeline annotation Primer design and real time quantitative PCR. Candidate reference genes and genes of interest were selected considering the annotation of the peptides reported in the annotated transcriptome of T. rotula FE80 (CCMP1647) 23 .
Oligo sequences utilized are listed below 23 : Each sequence was initially tested by standard PCR in a 25 μL final volume with 2,5 μL of 10× PCR reaction buffer (Roche, Basel, Switzerland), 2,5 μL of 10 × 2 mM dNTP, 0.3 μL of 5 U/μL Taq (Roche, Basel, Switzerland), 1 μL 10 μΜ of each oligo, 1 μL of cDNA templates and nuclease-free water up to 25 μL. Standard PCR amplification program was used, i.e. 95 °C for 3 min, 40 cycles at 95 °C for 30 s, 53 °C 30 s, 72 °C for 30 s, and a final extension step at 72 °C for 7 min. Amplified PCR products were analysed by agarose gel electrophoresis and verified by sanger sequencing.
Reverse transcription-quantitative PCR (rt-qPCR) experiments were performed in MicroAmp Optical 384-Well reaction plate (Applied Biosystems, Foster City, CA, USA) with Optical Adhesive Covers (Applied Biosystems, Foster City, CA, USA) in a Viia7 Real Time PCR System (Applied Biosystem, Foster City, CA, USA). Five serial dilutions of mixed cDNAs were used to determine primer reaction efficiency 23 using the formula: E = 10 −1/slope . The PCR volume for each sample was 10 μL, with 5 μl of SensiFAST TM SYBR ® Lo-ROX Kit (BIO_94020, Bioline), 1 μL of cDNA template (1 to 5 dilution each template) and 4 μL of 0.7 μM oligo mix (forward and reverse). Program reaction used was: 95 °C for 20 s, 40 cycles of 95 °C for 1 s and 60 °C for 20 s. The program was set to reveal the melting curve of each amplicon from 60 °C to 95 °C, and read every 0.5 °C. Single peaks for all genes confirmed gene-specific amplification and the absence of primer-dimers. All RT-qPCR reactions were carried out in triplicate to capture intra-assay variability. Each assay included three no-template negative controls for each primer pair.
The normalized expression levels of each gene of interest relative to the most stable reference genes 23 , actin and TBP, were calculated by using the Q-Gene tool 54 . Only TBP normalized values were reported in the main text and figure results. Relative expression ratios above two fold were considered significant 23 .
Prostaglandin extraction. One µg of Prostaglandin E2-d4 (CAYMAN CHEMICAL, Michigan, USA) was added as internal standard to the media (100-200 mL) recovered from the cell culture, from day 3 to day 10, through a filtration step on 1.2 µm nitrocellulose filters (Millipore RAWP04700). Culture media were loaded onto pre-packed CHROMABOND ® HR-X cartridges (500 mg/6 mL) previously activated with methanol (12 mL) and milliQ water (12 mL). After a preliminary desalting step with 12 mL of milliQ water, collection of the organic components was achieved by elution with 16 mL of methanol followed by 16 mL of methanol/dichloromethane (1:1). The two organic fractions were combined and dried under reduced pressure for LC-MS/MS analysis.

LC-MS/MS analysis.
Thalassiosira rotula water samples extracted as above, were re-suspended in 0.5 mL methanol and analyzed by tandem mass spectrometry in MRM (Multiple Reaction Monitoring) mode on a 4000 QTRAP ® LC/MS/MS System (Applied Biosystems, Toronto, Canada), working in negative ion mode and coupled to a 1100 nanoHPLC system (Agilent Technologies, Waldbronn, Germany). Prostaglandins were separated by using a micro C18 column (10 cm ×1.0 mm, 5 µm). The mobile phase was generated by mixing eluent A (water, 0.1% Acetic Acid) and eluent B (acetonitrile/isopropanol 50/50) and the flow rate was 30 nL/min. Elution started at 20% B up to 95% B in 15 minutes. Tandem mass spectrometry was performed using a turbo ion spray source operating in negative mode, and the multiple reaction monitoring (MRM) mode was used for the selected analytes. Mass parameters (4000 qtrap AB science) were as follows: curtain gas 20 psi, GS 1/2 50/50 psi, ion spray voltage −5500 V, DP, −60 V; Dwell Time 25 ms and Temperature of 550 °C. Quantitative analysis was performed by monitoring a unique production 26 arising from collision-induced fragmentation of the deprotonated selected parent compound after proper optimization of mass spectral parameters. After testing the scan mode including all the transitions, only the most intense transitions were chosen for each molecule.
Statistical analysis. One-way ANOVA (α = 0.05) with Tukey's post-hoc test was performed using GraphPad Prism6.0 (GraphPad Software Inc., San Diego, CA, USA). The analysis was performed to determine significant differences among the time points (days) considered during the T. rotula growth. Both qPCR amplification and LC-MS/MS data have been analysed.