Identification of differentially represented GO terms

The transcriptomes of L. danicus (strain SZN-B650) and L. aporus (strain SZN-B651) were obtained by sequencing samples collected at the exponential phase, under identical growth conditions for the two species. The transcriptomes of the two species were assembled, annotated and mapped as described in the supplementary material (Additional file 1: Sequence assembly and annotation, and Table S1–Table S4; Additional file 2: Figure S1; Additional file 3: Figure S2; Additional file 4: Figure S3; Additional file 5: Figure S4; Additional file 6: Figure S5). Comparative analysis, following Fisher’s exact test (FDR <0.05), led the identification of 10 GO terms that were identified as differentially represented between the two species (Additional file 7: Table S5A). These classes contained 3917 and 2142 transcripts from L. danicus and L. aporus, respectively. Nine of the ten GO terms were over represented in L. danicus, while in L aporus only ‘protein tyrosine kinase activity’ was over represented. We investigated further to obtain lower level GO terms that were differentially represented in the two species transcriptomes (Additional file 7: Table S5B). Among other differences, Leptocylindrus danicus showed a higher representation of genes associated with sexual reproduction such as the GO terms ‘cilium’, ‘cilium morphogenesis’, ‘female meiosis’, ‘female gonad development’, ‘female sex differentiation’, ‘development of primary female sexual characteristics’ and ‘sex differentiation’ (Additional file 7: Table S5B), although no specific conditions that could stimulate sexualization had been applied to the culture used for RNA sequencing. Interestingly, in L. aporus often the genes associated with these GO terms were absent, whereas genes associated with the GO terms ‘androecium development’, ‘stamen development’, ‘stamen formation’ and ‘stamen morphogenesis’ were found, which were almost absent in L. danicus. Genes involved in signaling were also over represented in L. danicus. These included the GO terms ‘notch signaling pathway’ and ‘smoothened signaling pathway’. The observed differences in the transcriptomes of the two species, in terms of genes associated with the flagellar development necessarily linked to the male gamete formation, paved the way for a deeper investigation into the nature of the annotated transcripts and their expression during sexual reproduction.

Sexual reproduction and flagellar genes expression in L. danicus

Sexual reproduction was induced in L. danicus by using a nutrient-depleted, modified F/2 growth medium (T medium) and incubating cells at a slightly lower temperature, 16 °C [26]. In a typical, successful experiment, 36 h from the first subculture in T medium cells started to undergo meiosis, which was recognized by the presence of ‘tetrad cells’ as opposed to cylindrical vegetative cells (Fig. 1a and b). By 48 h most of the cells in the culture were undergoing meiotic divisions, producing gametes within gametangia (Fig. 1c). Gametes were released and auxospores became visible by 60 h (24 h post meiosis initiation), while fully developed resting spores that originate from auxospores were observed after 72 h (36 h post meiosis initiation) (Fig. 1d–f). However, the time taken for the different processes and the proportion of cells undergoing sexual reproduction were variable among induction experiments conducted at different times.

In addition to what occurred upon nutrient depletion and exposure at lower temperature, sexual stages were also observed in cultures grown under normal conditions. With respect to cell size, we did not observe signals of sexual reproduction in cultures with cells having the maximum diameter (>10 μm). Additionally, in maintenance cultures, we observed that cells would automatically undergo cell size expansion after reaching minimal size (<4 μm in diameter) but whether it occurred through sexual processes or vegetative expansion was not ascertained.

To test specific gene expression during sexual reproduction, we selected six genes that have a core function in the development or the functioning of the flagellum among the 22 transcripts that were associated with the GO terms ‘flagellum’, ‘cilium’, ‘cilium assembly’, ‘motile cilium’ and/or ‘cilium morphogenesis’. The genes selected included: ‘cytoplasmic dynein 2 heavy chain’, ‘dynein heavy chain axonemal’, ‘intraflagellar transport protein 172’, ‘tubulin tyrosine ligase-like 1 protein’, ‘tetratricopeptide repeat protein 30a’, and ‘B9 domain containing protein 1’. The expression levels of these transcripts were assessed through quantitative real-time PCR analysis at four stages of sexual reproduction: early meiosis (12 h before tetrad cell formation), meiosis (tetrad cells), auxospore formation and auxospore-resting spore formation and compared to vegetative growth stage. All the selected genes were significantly up-regulated at meiosis (>2-fold change), while at other stages of the sexual reproduction process the expression level of the six genes generally fluctuated within a smaller range (within +2 to −2) (Fig. 2).

Flagellar gene distribution and phylogeny

To understand the possible evolutionary origin and relationships in the genes encoding for the flagella, we investigated the relationships among the algal species utilising sequences retrieved from their transcriptomes for the IFT172 and IFT88 protein sequences. Maximum likelihood (ML) trees constructed using RAxML showed that the amino acid sequences of the IFT172 protein generally formed clades or grades reproducing known phylogenetic relationships in algal groups (Fig. 3, Additional file 10: Figure S6). Similar results were also obtained for the IFT88 protein, but with fewer taxa (Additional file 11: Figure S7).

Flagellar gene in L. danicus during sexual reproduction

Although diatom transcriptome data are now available for many centric species, including those known to have sexual reproduction, the information on genes activated during sexual reproduction has so far been restricted to a few species for which specific experiments have been performed. In this respect, L. danicus is a special case of a species spontaneously undergoing sexual reproduction even in normal cultivation procedures, while other species may need specific induction procedures. A few flagellar genes had already been found in the genome of another centric diatom, Thalassiosira pseudonana [16], yet this first transcriptome study clearly shows the over representation of several flagellar genes in L. danicus compared to the congeneric species L. aporus not known to undergo sexual reproduction. Although the molecular functions of some of these genes and other associated genetic networks in the flagella development and sexual reproduction have not been fully elucidated, the present analysis suggests they are all involved in the flagellum development in male gametes.

In general, sex inducing factors in diatoms include decrease in cell size, environmental cues (for example temperature) and possible endogenous signaling molecules (sex pheromones) [6]. In L. danicus, as already reported [26], depleted nutrients (low N and Si) and lower growth temperature (16 °C) were found to be important factors that can stimulate sexualization. Nonetheless, like in other diatoms, getting laboratory cultures to undergo sexual reproduction may produce variable results, in terms of time taken for the process to occur and the magnitude at which it occurred, which probably depends on differences in the physiological status of the cells. In addition hallmarks of sexual reproduction in cultures grown under normal conditions, including the one used to obtain RNA for sequencing within the MMETSP project, indicate that sexual reproduction can occur in culture even uninduced.

The six genes chosen to explore expression in relation to sexual reproduction gene have been associated with specific processes related to flagellum development and functions. Cytoplasmic dynein 2 heavy chain, dynein heavy chain axonemal and intraflagellar transport protein 172 are three proteins involved in the bidirectional transport along axonemal microtubules, a process that is essential for the formation (ciliogenesis) and maintenance of most eukaryotic cilia and flagella [27, 28]. Tubulin tyrosine ligase-like 1 protein and tetratricopeptide repeat protein 30a are involved in the post-translational polyglutamylation of tubulin in axonemal microtubules [29, 30]. Finally, ‘B9 domain containing protein 1’ (which has a cilia-specific role) prevents diffusion of transmembrane proteins between the flagella and plasma membranes [31]. These genes were all found to be expressed in the early phases of the process corresponding to meiosis, and this is consistent with the notion that flagella develop in male gametes (sperms) of centric diatoms and specifically in Leptocylindrus [20, 21]. Additionally, our results indicate that there is only a narrow time window during sexual reproduction where flagellar genes are expressed.

Flagellar genes distribution in diatoms and other protists

Since the transcriptomes of diatoms were not purposely obtained from cultures undergoing sexual reproduction, the absence of flagellum related-transcripts does not necessarily imply their absence in the genome of the species sequenced within the MMETSP. On the other hand, transcripts related to flagellar development in the transcriptome of several species for which sexual reproduction has not been described as yet (e.g., Proboscia alata and Rhizosolenia setigera) suggest that this process is widespread in diatoms. Their presence in so many cases might be due to background levels of transcription or, alternatively, to sexual reproduction occurring in the cultures undetected.

It is hypothesised that diatoms have lost both BBSome and IFT-A components during evolution, retaining only IFT-B components. In the current study, BBSome complexes were not found in any diatoms, confirming the hypothesis of their complete loss in this group. However, in contrast with the current hypothesis [16], the IFT-A subcomplex has been retained by centric diatoms, or at least by a number of them. In fact, L. danicus, L. hargravesii, Proboscia alata, Rhizosolenia setigera, and Thalassiosira weissflogii showed both IFT122 and IFT140, while Aulacoseira subarctica, Odontella aurita, Stephanopyxis turris and Ditylum brightwellii had IFT122 (Table 1 and Additional file 8: Tables S6 and Additional file 9: Table S7). Recently, IFT-A (IFT-140) was also reported in the brown algae Colpomenia bullosa and Ectocarpus siliculosus, which form flagellated gametes [32], corroborating the idea that IFT-A subcomplex is more widespread than thought so far. Within the IFT-B complex, many of the centric species were found to have IFT52, IFT80, IFT88 and IFT172 including L. hargravesii for which sexual reproduction has been reported [17].

While the centric diatoms contain both IFT-A and IFT-B subcomplexes, pennate diatoms do not seem to retain any of these proteins, which is consistent with the absence of flagellated gametes during sexual reproduction in these organisms. Indeed, the unexpected finding of IFT-A and IFT-B transcripts in the pennate diatom Pseudo-nitzschia fraudulenta WWA7 was eventually attributed to a contamination by a dinoflagellate in the culture submitted to MMETSP for transcriptome sequencing. Other transcripts related to flagellar functions were found in many pennate diatoms. However, because those proteins have additional roles, they do not represent a univocal indication of functional flagellum machinery.

Interestingly, the phylogenetic relationship for flagellar genes in protists displayed similar pattern as the ribosomal genes. For example, the diatoms species were found within the Stramenopiles and were sister to Bolidomonas spp., which matches what is generally seen in the phylogenetic trees based on ribosomal genes [17]. Thus it appears that the flagellar gene sequences are conserved within the groups and their evolution followed the general phylogenetic pattern expected from the tree of life.

Cell culture

Fresh cultures of L. danicus strain SZN-B650 and L. aporus strain SZN-B651 (the strain was submitted as L. danicus strain SZN-B651 as only later it has become clear that it was a separate species, [17, 33]) were inoculated at a starting cell density of 200 cells ml⁻¹ in 250 ml polystyrene culture flasks with 100 ml of Keller medium. Cultures were maintained at 20 °C with 100 μE m⁻² s⁻¹ light and a diurnal cycle of 12:12 h light: dark periods. The cultures were grown till they attained a cell density of 10⁴ cells ml⁻¹ and later diluted to 1 L to give a cell density of 10³ cells ml⁻¹. Then the cultures were further incubated under the same conditions for 3 more days. In this way the cells were always retained in the log phase while reducing the time of growth in larger volume. Once they attained the density of 10⁴ cells ml⁻¹, the cultures were harvested by filtering on 1.2 μm pore-size nitrocellulose filters (Millipore, Italy). The circular filter was cut into four parts and each quarter was dipped in 1 ml of TRIzol®-reagent (Life Technologies Europe, Italy). The samples were either immediately processed or frozen in liquid nitrogen and stored at −80 °C for later RNA extraction.

RNA extraction and sequencing

Frozen samples were thawed on ice and homogenised by vortexing for 2–3 min. The samples were then incubated on thermo shaker at maximum speed for 10 min at 60 °C, to help cell disruption and release of RNA. Subsequently, the samples were cooled to room temperature and 300 μl of chloroform was added under fume hood and was vigorously mixed by inverting the tubes. The tubes were incubated at room temperature for 15 min for complete precipitation. Thereafter, the samples were centrifuged at 10,000 g for 15 min at 4 °C. The aqueous supernatant (≈750 μL) above the white phase was carefully transferred to new RNase free pre-chilled 1.5 ml Eppendorf tube. Equal volume (750 μL) of isopropanol was added to the collected supernatant. The RNA in supernatant-isopropanol mix was allowed to precipitate under room temperature for 10 min. The mix was then centrifuged at 10,000 g for 10 min at 4 °C to pellet RNA. The obtained RNA pellet was washed once with 750 μL of 75% ethanol (in DEPC treated water), by centrifuging at 10,000 g for 10 min at 4 °C. Finally, the RNA pellet was resuspended in 50 μL of RNase free water and was stored in deep freezer.

Total RNA extracted using the above method was purified with an RNeasy kit (Qiagen), following the manufacturer’s protocol. Measures of RNA quality (purity and integrity) and abundance were determined for each sample using a NanoDrop Spectrophotometer (Thermo Fisher Scientific Inc., UK) to determine the A260/A280 ratio and RNA concentration, and the Agilent Bioanalyzer 2100 RNA Nano LabChip (Agilent, Palo Alto, CA) to generate an electropherogram and an RNA Integrity Number (RIN). A good quality RNA sample was used for total transcriptome sequencing.

Total RNA was submitted to the National Center for Genome Resources (NCGR, Santa Fe, NM) for sequencing under the project, Marine Microbial Eukaryote Transcriptome Project (MMETSP; Gordon and Betty Moore Foundation). RNA libraries were prepared from total RNA isolated using the standard Tru-Seq™ RNA protocol - poly-A+ selection (Illumina, USA) with an insert size of ~200 bp. Sequencing was done from both ends (paired-end reads 2 × 50 nt) on the Illumina Hi-Seq 2000 to obtain approximately 2.5 Gbp of data (L. danicus, MMETSP0321 and L. aporus, MMETSP0322). Reads were filtered for quality (Q15) and read length (0.5). The two transcriptomes can be downloaded from the iMicrobe Data Distribution Center [34], in which all MMETSP data have been deposited.

Assembly and annotation

Sequences were assembled using AbySS [35] to generate 20 assemblies across a k-mer sweep from 26 to 50 nt. The obtained assemblies were then subjected to gap closing with GapCloser v 1.10 (Short Oligonucleotide Analysis Package, SOAPdenovo) [36]. Gap closed assemblies from the k-mer sweep were reconciled into one assembly using miraEST [37]. The assembled transcripts were annotated using blastx search in Blast2GO [38] against the NCBI non-redundant protein database with 10⁻⁶ e-value cut-off. Sequences were mapped and functionally annotated on Blast2GO with the default settings. GO functional classification for transcripts was performed using WEGO (Web Gene Ontology Annotation Plot) [39]. KEGG metabolic pathway annotation and COG classification of transcripts were determined through blastx search against KEGG database and COG database, respectively [40, 41]. The difference in GO term abundance between the two species was analysed using Fisher’s exact test with the GOSSIP package [42] implemented in Blast2GO, with support from p-value and subsequently the false detection rate (FDR) at 0.05.

Inducing sexual reproduction

Meiosis was induced in cultures by reduction of the nutrient levels in the medium, designated as “T” medium (containing ammonia at 15 μM; phosphate at 7 μM; silicate at 50 μM; trace metals, vitamins, and Trizma as in f/2) [43], with at a slightly lower temperature of 16 °C, light of 80 μmol photons m⁻² s⁻¹ intensity and a 12:12 L:D cycle [26]. The strains submitted to these reduced nutrient conditions were observed regularly in the LM. On cultures where spore formation was observed, microscope pictures were taken and samples were collected at different stages for RNA extraction and real-time PCR analysis.

Real-time PCR assays

Six annotated transcripts that may relate to functionally-conserved genes for flagella development in centric diatoms were selected and analysed using real-time PCR (specific primers are listed in Additional file 12: Table S8). For real-time PCR analysis, sexual reproduction in Leptocylindrus danicus culture was inoculated into a 200 ml polystyrene flask with 150 ml T media [44] to a cell density of 2000 cells ml⁻¹ and incubated at 16 °C with 100 μE m⁻² s⁻¹ light and a photoperiod of 12:12 light:dark. Samples were collected at various stages of sexual reproduction and vegetative growth. Cultures were filtered and RNA was extracted following the same protocols described above. cDNA was synthesized from the extracted total RNA with SuperScript II RT (Invitrogen, USA), and quantitative real-time PCR was triply performed using the Applied Biosystems with a volume of 10 μL containing 1 μL of 1:5 cDNA diluted with RNase free water, 5 μL of 26SYBRGreen Master Mix, and 250 nM of each primer. The cycling parameters included 95 °C for 2 min, followed by 40 cycles of 95 °C for 1 s, 62 °C for 1 min, and 70 °C for 1 s. During each cycle fluorescence measurements were taken at 70 °C for 1 s. All the qPCR standardisations were performed using the cDNA prepared using cultures in vegetative growth. Expression levels of each gene for sexual reproduction were normalized to histone H4, transcription initiation factor subunit 12 and translation elongation factor 1 alpha. Statistical analysis was performed with the relative quantification method described by Muller et al. [45].

Homology based searches

tblastn searches were conducted using L. danicus protein sequences as queries against the nucleotide collection of MMETSP currently available. Only hits with an e-value smaller or equal to 1e-10 were considered. Blast results for diatoms only are shown in Additional file 7: Table S5. Peptide sequences for selected protein hits were retrieved from iMicrobe [34] and blasted in public databases (blastp in nr) to confirm protein identity.

Phylogenetic analysis

To study the evolutionary pattern of the flagellum coding genes in diatoms and other algae for which the sequence information is available, we applied maximum likelihood relationship approach. Sequences for the IFT172 and IFT88 genes involved in the flagellum development identified in L. danicus were downloaded from NCBI nucleotide database through homology searches using blastn. Sequences were also obtained by local BLAST searches from the transcript dataset that constituted many publically available transcriptomes of species sequenced within the Marine microbial eukaryotic transcriptome sequencing project (e-value smaller or equal to 1e-10). Obtained sequences were gathered into datasets and sequence alignment was generated using MUSCLE version 3.8.31 [46]. Fewer taxa were used for IFT88 than IFT172 analysis. The IFT172 dataset was comprised of 137 sequences which were aligned to give an alignment length of 2732 bp. The IFT88 dataset included 99 sequences that were aligned with an alignment length of 1324 bp. Any misalignments were corrected manually and phylogenies for each of the dataset was constructed following a maximum-likelihood approach [47] using the WAG amino acid substitution model [48] with a discrete gamma distribution [49] and nonparametric bootstrap of 1000 replicates.