Mating type specific transcriptomic response to sex inducing pheromone in the pennate diatom Seminavis robusta

In this study, we generated a new time-series dataset to investigate the response of dark-synchronized MT+ cultures to SIP− filtrate (6 time points, 0–9 h, Fig. 2a) and to compare expression patterns with two existing datasets on the response of MT− to SIP+ [16, 18] (4 time points, 0–10 h, Fig. 2b). Here, we first report the results of a conventional DE analysis for each dataset separately. Next, we use these results to discover genes or biological processes related to sexual reproduction, either based on their predicted functional annotation or the current literature. Finally, we integrate the results of both the novel and publicly available datasets in an integrative analysis that aims at identifying key genes involved in either a single or both mating types.

General transcriptional response and identification of key SIP responsive genes

Multidimensional scaling plots of the RNA-seq data (Fig. 2a, b) showed that in both mating types the dark-to-light transition and time since illumination were the major drivers of gene expression change throughout the experiments. However, as time progresses, the effect of SIP becomes more pronounced. This is supported by DE analysis showing that the number of significant genes increased markedly at later time points (Fig. 2c). Overall, most genes are significantly DE in only a single time point on a 5% overall FDR (OFDR) level (Supplementary Fig. 1B). Combined, on a 5% OFDR level, 4037 genes are DE in response to SIP treatment for MT+, while 5486 genes are found to be DE in MT− in the first 3 h and 6079 genes after 10 h. The stronger response during the first 3 time points in MT− versus MT+ may be the result of the use of a chromatographic fraction of SIP+ for MT− while MT+ cultures were treated with a SIP− containing filtrate. Gene sets for many biological processes were significantly enriched in the conventional lists of DE genes: a total of 1081 and 740 enriched biological process terms were discovered for the response to sex pheromones in MT+ and MT−, respectively (Fig. 2d, Supplemenatry Fig. 2).

We discovered key genes by developing a statistical integrative analysis workflow that is capable of testing for equivalent, i.e., non-DE, expression between conditions. Coupling equivalence testing in one mating type with DE calls in the other allowed for the discovery of key genes exhibiting mating type specific responses to SIP, while DE calls in both datasets found key genes responsive in both mating types. This workflow revealed 52 key genes responding to SIP in both mating types (SRBs), 12 genes uniquely responding in MT+ (SRPs) and 70 genes uniquely responding in MT− (SRMs) (Fig. 3a, Supplementary Figs. 3, 4, 5). Similar to the conventional DE analysis, the response of MT− was more pronounced compared to MT+, likely due to technical differences such as different protocols for pheromone administration. Remarkably, while in MT− we discovered a comparable number of down- and upregulated SRMs, we only found upregulated SRPs and SRBs, possibly indicating that sexual processes induced by SIPs are mainly driven by the induction of key genes rather than the downregulation of inhibitory genes (Fig. 3a).

In what follows we will discuss in more detail the genes and pathways that are responding to SIP in both mating types or uniquely in one mating type and link these changes to physiological events in the mating process. In each section, we first discuss key genes highlighted by the integrative analysis (Fig. 3), after which we discuss results from the conventional DE analysis, focussing on selected biological processes (Fig. 4).

Responses to SIP conserved in both mating types

Integrative analysis reveals key genes responsive in both mating types

A large fraction (22/52) of SRBs, i.e., key genes with a strong response to SIP in both mating types, lack any functional annotation and homology to sequenced genomes of other diatoms (Supplementary Table 2), suggesting that the molecular mechanisms underlying early mating are highly species-specific. The remaining 30 SRBs can be linked to energy metabolism, ROS signaling and meiosis, amongst others. Pyruvate kinase (Sro373_g129070) and isocitrate dehydrogenase (Sro492_g153950), respectively involved in glycolysis and the citric acid cycle, are strongly upregulated in both mating types (Fig. 3b), suggesting an increased energy demand. Interestingly, pyruvate kinase is also upregulated during gametogenesis in the brown alga Saccharina latissima [26] and the parasite Plasmodium berghei [27]. In addition, two enzymes from the pentose phosphate pathway (PPP) are among the SRBs: transketolase (Sro524_g159900) and transaldolase (Sro196_g083630) (Supplementary Fig. 3). The PPP generates NADPH, a reductive compound needed in various metabolic reactions and involved in detoxification of ROS by regenerating glutathione [28, 29]. Furthermore, one SRB encoding a heme peroxidase (Sro1252_g256250) exhibited strong upregulation upon SIP treatment (Fig. 3b). Upregulation of heme peroxidases was also reported during sexual reproduction in other eukaryotes, e.g., mosquitoes (Anopheles gambiae) [30] and fungi [31, 32]. Heme peroxidases promote substrate oxidation in various metabolic pathways and are essential for the detoxification of ROS [33], suggesting that ROS signaling plays a role in the response to SIP, as seen in the green algae Volvox carteri, where high ROS levels trigger sex [34]. Finally, a highly expressed SRB encodes a transmembrane protein containing an Epidermal Growth Factor-like (EGF-like) domain (Sro65_g036830, Fig. 3b), with potential orthologs encoded in pennate and centric diatoms including P. tricornutum and T. pseudonana (BLASTp, E < 1e-10). EGF-like domains are generally extracellular protein modules that play a role in receptor/ligand interactions, intracellular signaling, and adhesion [35]. Membrane bound proteins containing EGF-like domains are required for gamete fusion in green algae [36] and oocyte binding by animal sperm cells [37,38,39]. Furthermore, multiple EGF-like repeats were discovered in sexually induced genes Sig1–3 of the centric diatom Thalassiosira weissflogii [40]. Sig1 orthologs were later shown to be located on the mastigonemes of stramenopile flagella [41], suggesting their upregulation is related to the differentiation of flagellated male gametes in centric diatoms. We could not detect orthologs of Sig1–3 in the genome sequence of S. robusta (BLASTp, E < 0.001), in line with the lack of flagellated stages in pennate diatoms. Nevertheless, the presence of EGF-like domains in sexually induced genes in pennate and centric diatoms suggests that genes containing EGF-like domains play a role in diatom cell–cell communication or gamete fusion.

Conventional DE analysis uncovers mating-related processes differentially regulated in both mating types

We used flow cytometry to confirm a sex pheromone induced G1 phase arrest, which was proposed for S. robusta and Pseudo-nitzschia multistriata [16, 17]. Treatment with SIP significantly decreased the proportion of G2/M phase cells in MT+ after 3 h (control 12.4% versus SIP 3.5%, p = 0.0023) and 9 h (control 20.9% versus SIP 9.4%, p = 0.0123) and after 9 h in MT− (control 10.5% versus SIP 3.3%, p < 0.0001) (Fig. 5a, b). A temporary arrest of the cell cycle appears to be a prerequisite for a switch to meiosis in many diatoms [16, 17, 42, 43], although in the centric diatom Skeletonema marinoi no growth arrest was observed during sexual reproduction [44]. The cell cycle arrest is reflected in the transcriptomic data as a sequential downregulation of cell cycle genes in the conventional DE results (Fig. 4a), causing an enrichment in cell cycle related GO terms (Supplementary Fig. 2). In eukaryotes, S phase progression is controlled by E2F transcription factors forming a heterodimer with Dimerization Partner (DP) transcription factors [45]. In accordance with plants and animals [46, 47], we observe an increase in expression of E2F transcription factors in control conditions as cells go through the cell cycle. SIP treatment significantly repressed expression of two E2Fs (Sro2696_g334910 and Sro1798_g298290) in both mating types and MT−, respectively (Supplementary Fig. 6A). The transcriptional repression of E2Fs is likely caused by the inability of cells to enter S phase as a result of the G1 phase arrest, although the activity of E2Fs is generally also regulated by other mechanisms such as Retinoblastoma-related (Rb) protein binding and phosphorylation which we did not assess [46, 47]. Furthermore, two DP genes were DE in MT− in response to SIP+; one (Sro905_g218540) was repressed, while the other (Sro905_g218570) was induced, suggesting they play a contrasting role in the cell cycle arrest (Supplementary Fig. 6A).

During the S and G2 phase of the mitotic cell cycle, fission of both chloroplasts of S. robusta results in four daughter chloroplasts by the start of the M phase [48]. FtsZ and dynamin-related protein 5B (DRP5B) are key factors in the formation of a multi-ring structure that constricts the chloroplast during fission [49, 50]. In control conditions, we observe a expression peak of two FtsZ homologs (Sro1409_g270150 and Sro931_g221490) and DRP5B (Sro814_g206290), potentially coinciding with the timing of chloroplast division in S. robusta [48] (Supplementary Fig. 6B). After exposure to SIP, expression of these genes was significantly repressed in both mating types (Supplementary Fig. 6B), indicating that the cell cycle arrest results in an inhibition of chloroplast fission. Indeed, chloroplast division ceases during sexual reproduction of S. robusta so that each gamete will contain one chloroplast and the auxospore inherits two chloroplasts [51]. Noteworthy, while DRP5B is restricted to the chloroplast [52], stramenopiles can encode a mitochondrium-targeted FtsZ that is involved in mitochondrial division [53]. Thus, it is possible that the downregulation of FtsZ genes is linked with mitochondrial instead of chloroplastic division.

In the transcriptomic data, seven cyclins were downregulated in both mating types after treatment with SIP (Fig. 4a). Six out of seven appear in a cluster of genes, which peak late in the time series in control conditions (6 h, 9 h and 10 h; Supplementary Fig. 7), suggesting that they are mitotic cyclins involved in G2/M transition [54]. To determine which cyclin family they represent, a maximum likelihood phylogenetic tree of S. robusta cyclins was constructed (Supplementary Fig. 8). The repressed cyclins consist of three diatom-specific cyclins (dsCyc), three A/B type cyclins and one cyclin D (Supplementary Table 3).

Several genes involved in mitosis were downregulated in both mating types, including MAD2 (Sro3109_g343920), MAD3 (Sro690_g187690) and CDC20 (Sro589_g171680), which form the Mitotic Checkpoint Complex. Furthermore, we observe downregulation of five genes coding for subunits of the condensin complex, which play a central role in chromosome organization during mitosis and meiosis [55] (Fig. 4a).

Finally, at the end of each mitotic cell cycle, diatom cells produce a new valve in a silica deposition vesicle (SDV) prior to cytokinesis. As it was previously shown [16], treating cultures with SIP of the compatible mating type reduces the fraction of cytokinetic cells after 14 h (Supplementary Fig. 9A, p < 0.0001 for both MT+ and MT−). Accordingly, we observed downregulation of genes known to be important for silica cell wall formation after treatment with SIP (Fig. 4a), including two out of five S. robusta silicic acid transporter homologs, involved in the uptake of silicic acid from the environment [56]. In addition, we find two (out of a total of three) silicanin1-like genes as well as four (out of thirteen) frustulins, respectively proteins embedded in the SDV membrane and proteins which are found in the organic casing surrounding the cell wall [57, 58]. Moreover, four genes making up various subunits of V-type ATPase complexes were significantly downregulated in both species upon treatment with SIP. These proton pumps play a role in biomineralization of silica by acidifying the SDV [59].

While mitotic cell cycle genes were generally downregulated by SIP treatment, we observed an increase in meiotic gene expression, accompanied by an enrichment in the “meiotic chromosome condensation” GO term in both mating types (Supplementary Fig. 2). Surprisingly, one cyclin which is an SRB (Sro299_g111470, Fig. 4a) is upregulated rather than downregulated in response to the pheromone. Its expression pattern suggests it plays a role either in SIP induced physiological responses such as mate finding or in preparing the cell for meiotic cell cycle progression. The gene has presumably evolved independently from the sexually induced cyclins of other major eukaryotic clades, as phylogenetic analysis places the gene among diatom specific cyclins (dsCyc) (Supplementary Fig. 8). The expansion of the cyclin family in diatoms compared to other members of the SAR clade [60] might have been instrumental to allow the diversification of this sexually induced cyclin [61,62,63,64].

To investigate the dynamics of other meiotic genes throughout the experiment, we explored the expression of a set of 42 known meiotic and bifunctional mitotic/meiotic diatom genes [5], some of which were induced in response to sex pheromones in pennates [16, 17] and during sexual reproduction in the centric diatom S. marinoi [44]. Here, we found ten meiotic markers significantly responding to SIP in at least one time point in both mating types (Fig. 4a). These include genes encoding DNA replication licensing factors MCM4, MCM6, MCM7, which are involved in the initiation of replication during the mitotic and meiotic S phase [5]. Two homologues of SPO11 (SPO11–2 and SPO11–3) are upregulated in both mating types. SPO11 plays a role in the formation of double strand DNA breaks during homologous recombination. In diatoms and plants, the SPO11–2 homologue was shown to be meiotic while SPO11–3 is involved in vegetative growth [5, 17, 65]. However, the observed upregulation of SPO11–3 in response to SIP in S. robusta and during sexual reproduction in the centric S. marinoi [44] suggest that in addition to SPO11–2, SPO11–3 might also be involved in meiotic homologous recombination in diatoms. Three genes involved in DNA repair after the induction of double stranded breaks by SPO11 were upregulated in both mating types: MRE11, RAD50, and RAD51D [5]. Finally, we observed the upregulation of two genes not yet described in the meiotic toolkit of diatoms: Ataxia Telangiectasia Mutated (ATM, Sro154_g070060) which codes for a protein that controls double-strand break formation by SPO11 during meiosis [66] and Replication Protein A (RPA, Sro944_g222890), involved in the binding of ssDNA during replication and homologous recombination [67] (Fig. 4a). Interestingly, although meiotic genes are upregulated in response to SIP, the G1 phase arrest inhibits progression through the meiotic cell cycle and we did not observe gametogenesis in MT+ cells attracted to diproline loaded beads. Accordingly, chloroplast rearrangements and indications of the meiotic prophase were only observed microscopically after compatible cells form a mating pair [51]. We therefore hypothesize that a separate, local signal during cell pairing is required to break the G1 phase arrest and induce gametogenesis.

In response to SIP, 17 genes with a guanylate cyclase domain (GC) were significantly DE in both mating types, 8 of which form a bifunctional guanylate cyclase/phosphodiesterase (GC/PDE) fusion enzyme. Interestingly, the GC and PDE domains in these genes have contrasting functions, respectively synthesizing and breaking down the secondary metabolite cGMP. Although the genes with only a GC domain do not show a general direction of regulation, all GC/PDE genes were upregulated (Supplementary Fig. 10), including Sro991_g228730, which was previously described to be responding to SIP in S. robusta [16]. All significant GC/PDEs show the topology described by Moeys et al. (2016) [16] a small N-terminal intracellular domain, an extracellular stretch and a long C-terminal intracellular part containing the GC and PDE domains. Thus, cGMP signaling appears to be a conserved response to sex pheromones in pennate diatoms, as guanylate cyclases are also upregulated in P. multistriata [17].

Finally, orthologs of three uncharacterized genes induced during sexual reproduction in the centric diatom S. marinoi and the pennate diatom P. multistriata [44], were upregulated in response to SIP (Supplementary Fig. 6C). Sro587_g171310 (orthologous to STRINITY_DN12692_c0_g1_i1 and Pmu0010180) was significantly upregulated in both mating types. Protein domain analysis revealed the presence of a Homologous-pairing protein 2 (Hop2) domain in its S. robusta (IPR010776), P. multistriata (IPR010776) and S. marinoi (IPR040461) sequence. Although Hop2’s role in homologous recombination would be in accordance with the observed expression, more work is needed to confirm this gene as a Hop2 homologue, as Hop2 is presumed to be absent in diatoms [5, 44]. Another conserved sexual gene (Sro131_g062230, orthologous to Pc15065_g1_i1 and Pmu0009930), is upregulated in response to SIP only after 10 h of SIP treatment in MT−. Finally, Sro637_g179400 (ortholog of MTRINITY_DN9343_c0_g1_i2 and Pmu0061540) was upregulated after 3 h and 10 h in MT−. Interestingly, domain predictions show that the protein consists of an unknown N-terminal domain, followed by one transmembrane helix and a C-terminal beta-propeller domain in all three species (IPR013519, IPR011043 and IPR015943 in S. robusta, P. multistriata and S. marinoi respectively). This topology suggests a function as a receptor or in adhesion [68]. Since these three genes are upregulated in two pennate and one centric diatom species, they are interesting candidate marker genes for sexual reproduction in diatoms. Notably, we discovered five additional conserved sexual genes in the S. robusta genome which are not differentially expressed to SIP. If their function is conserved, we expect these genes to be upregulated during zygote or auxospore formation rather than during pheromone signaling, since the physiology of mate finding strongly differs between species.

Responses to SIP specific for MT+

We identified 12 genes displaying a MT+ specific response to SIP (SRPs, Fig. 3a, Supplementary Table 2), among them a Myb-like transcription factor (Sro94_g048900), which is likely regulating downstream events associated with the perception of SIP− by MT+. One such downstream response unique for MT+ is responsiveness to the attraction pheromone diproline [16], which we confirmed using a bead attraction assay (Supplementary Fig. 9B). This might be caused by SIP− induced expression of a diproline receptor in MT+. While biosynthetic pathways of 2,5-diketopiperazine cyclodipeptides such as diproline have been elucidated, cyclodipeptide receptors are not well characterized [69]. Among the SRPs, we have identified two functionally unidentified genes which satisfy the requirements for a potential diproline receptor: Sro1719_g293450 and Sro1_g000260. Both encode transmembrane proteins and show high expression uniquely in response to SIP− in MT+ (Fig. 3b). Interestingly, the former gene is predicted to contain 7–8 transmembrane helices, reminiscent of G-protein coupled receptors (GPCR) which often associate with cyclodipeptides in humans [69] and, moreover, a GPCR was upregulated in response to sex pheromones in P. multistriata [17]. Sexual S. robusta MT− cells were shown to predominantly secrete the cyclo(L-Pro-L-Pro) enantiomer [14]. However, while pheromone receptors are usually stereosensitive, synthetic cyclo(D-Pro-D-Pro) is also biologically active, suggesting a highly unusual stereo-insensitive receptor or racemization of the pheromone [14, 70]. Further, the SRPs include 6 genes lacking functional annotation, among which is SRP12 (Sro2882_g339270) that belongs to a gene family containing five S. robusta genes of which three are located adjacently on the same genomic contig (Sro2882_g339270, Sro2882_g339280 and Sro2882_g339290). Four out of five show a MT+ specific upregulation to SIP− with very high expression ranging up to 1.5% of the total transcriptome library (Fig. 3c). Furthermore, all five genes exhibit a high baseline MT+ expression and negligible expression in control MT− conditions (Fig. 3c), suggesting they play a role in mating type differentiation, comparable to MRP genes in the pennate P. multistriata [71].

When exploring ROS signaling related genes, the conventional DE analysis uncovered an NADPH oxidase (NOX, Sro1572_g283390) with a pronounced MT+ specific response to SIP (Fig. 4b), catalyzing extracellular production of superoxide anions by moving an electron from NADPH through the plasma membrane to molecular oxygen. The observed NOX contains six transmembrane domains [72] and an EF-hand domain suggesting a potential link with calcium signaling. As superoxide is cell impermeable, extracellular superoxide is most likely dismutated to H₂O₂, which can enter neighboring cells for ROS signaling, e.g., through aquaporin channels [73, 74]. NOX activity during sexual reproduction is a common theme in eukaryotes: it is required for gametogenesis and fertilization in plants [75] and is required for the formation of sexual fruiting bodies in the fungus Aspergillus nidulans [32]. Furthermore, gametophytes of the kelp S. latissima show female specific expression of NOX, suggesting mating type specific expression of NOX during sexual reproduction may be conserved among stramenopiles [26]. Notably, enzymes from the NADPH producing PPP were upregulated in both mating types, potentially supplying NADPH to support NOX-mediated superoxide production.

To assess whether SIP− affects the motility of MT+, the displacement of control and SIP− filtrate treated cells was tracked for 30 sec over the course of 9 h (Supplementary Fig. 11). We observed mate-searching behavior starting 6 h–7 h after treatment with SIP−. Not only did the fraction of motile cells (i.e., cells showing any displacement) in the culture increase (Omnibus test p < 0.001, Fig. 5c, Supplementary Fig. 12A), the speed of the subset of motile cells also increased: after 9 h, the average speed of motile SIP− treated cells was 3.8 µm/s compared to 1.1 µm/s for untreated cells (Omnibus test p < 0.0001, Fig. 5d, Supplementary Fig. 12B). Thus, exposure to SIP does not only prime MT+ cells to become responsive to the attraction pheromone diproline, the increased motility of gametangia further increases the probability to encounter an immotile, diproline producing MT− cell. Furthermore, the pre-activation of cellular motility machinery by SIP− might explain the almost immediate attraction to a new diproline source [76]. The timing of the mating behavior in S. robusta appears to be remarkably synchronized between mating types: significant increases in motility occurred simultaneously with the first noticeable amount of diproline produced by MT− [14] and the onset of responsiveness of MT+ to diproline (Supplementary Fig. 9B). As cell motility in raphid diatoms is achieved through the secretion of adhesive molecules that attach the cell to the substratum and provide traction for their gliding movement [77], we performed BLAST searches to identify adhesive proteins containing a GDPH-domain, named after a conserved Gly-Asp-Pro-His amino acid motif [78]. Among 87 GDPH-containing proteins identified, four are upregulated in MT+ in response to SIP− treatment in the conventional DE analysis (Supplementary Fig. 6D). However, three out of four GDPH-domain-containing genes are also significantly upregulated in immotile SIP+ treated MT− cells and one is upregulated uniquely in MT−. Therefore, as previously suggested, the GDPH-domain containing proteins may play additional roles related to extracellular adhesion in addition to cell motility, such as mucilage pad, stalk and chain formation as well as cell–cell adhesion during mating [78].

Responses to SIP specific for MT−

A relatively high number (70) of SRMs, i.e., genes with a MT− specific response to SIP, were discovered (Fig. 3a) of which 16 lack any functional annotation (Supplementary Table 2). Molecular functions of SRMs are diverse, including receptors, membrane channels, guanylate cyclases and other signaling enzymes. Among the SRMs are two E3 ubiquitin ligase genes (Sro1305_g261220 and Sro25_g017160) representing the Ring-Between-Ring (RBR) and U-box family, respectively. Since ubiquitin ligases are important players in signaling by targeting downstream proteins for ubiquitination [79], these genes are potential key regulators of MT− specific responses such as the production of diproline. Ubiquitin ligases also play a role in meiosis across eukaryotes by targeting proteins for proteasomal degradation [80]. However, a function in meiosis for the identified ubiquitin ligases is hard to reconcile with their mating type specific response, as meiosis occurs in both partners. SRMs further include a gene with a DOMON domain, which is a presumed heme or sugar sensor domain (Sro7_g006210) [81] (Fig. 3b). Interestingly, two SRM genes belong to an unknown and S. robusta specific gene family containing a zinc finger domain (Supplementary Fig. 6E).

Our conventional DE analysis also shows that glutamate-to-proline conversion enzymes Δ1-pyrroline-5-carboxylate synthetase (Sro2012_g310890, P5CS), and Δ1-pyrroline-5-carboxylate reductase (Sro216_g089310, P5CR) were significantly upregulated after SIP treatment in MT− after 3 h and 10 h, while in MT+ both genes were not significantly responding (Fig. 4b). The unique increase of proline biosynthesis in MT− is further supported by the enrichment of the “proline biosynthetic process” GO term in the 3 h time point only in MT− (Supplementary Fig. 2). Their MT− specific response supports the hypothesis that upregulation of P5CS and P5CR increases the cellular proline pool as a precursor for diproline biosynthesis [16]. Furthermore, a proline-tRNA ligase (Sro27_g018090) exhibits a strong significant upregulation uniquely in MT−, with expression levels exceeding 4000 CPM (Fig. 4b). This enzyme attaches proline to transfer RNA (tRNA), which serves as a substrate for ribosomal protein synthesis, explaining its consistent expression in control conditions in this dataset (Fig. 4b). The cyclodipeptide ring of diketopiperazines such as diproline is typically synthesized by either nonribosomal peptide synthetases or cyclodipeptide synthases (CDPS) [82]. Interestingly, CDPS require aminoacyl-tRNA as a substrate for the reaction instead of a free amino acid [83]. Thus, the observed MT− specific upregulation of a proline-tRNA ligase likely caters to the increased need for proline-tRNA^pro driving CDPS dependent diproline biosynthesis. BLAST searches in the S. robusta genome revealed several candidate CDPS genes, with rather low levels of conservation (Supplementary Table 4), of which none show significant regulation at any time points in any mating type. Thus, either an unidentified, transcriptionally controlled CDPS exists or one of the identified CDPS is non-transcriptionally regulated to be active only in SIP+ treated MT− cells with a cell size below the SST.

Source: Ecology - nature.com