Non-axenic cultures of A. minutum were obtained from the Collection Center of Marine Bacteria and Algae, Xiamen University, China. Cultures were maintained in F/2 medium (without silicate) (Guillard, 1975). Natural seawater (33 ppt salinity) was used to prepare the F/2 medium; the seawater was obtained from Dongchong (Shen Zhen), filtered through 0.22-μm filter, and autoclaved. Cells of A. minutum were grown in F/2 medium at 23 ± 1°C under 12 h light/12 h dark photoperiod and 60–70 μmol photons m⁻² s⁻¹ light intensity.

Alexandrium minutum cells in the logarithmic phase were collected as the inoculum and transferred to 120 ml sterile dialysis bags. All dialysis bags were placed in 10 L beakers filled with artificial seawater (the method of Schleicher and Schal) without N and P. The artificial seawater was exchanged at 2-h intervals for the first 12 h. The A. minutum cells were harvested using a 30-μm mesh at 0, 6, and 72 h, frozen in liquid nitrogen, and stored at −80°C until needed for cDNA library construction. Samples harvested at 0, 6, and 72 h are hereafter referred to as R, S6, and S72, respectively. This experiment was performed in triplicate at each time point.

Total RNA was extracted from nine samples (three treatments × three replicates) using the E.Z.N.A.® Plant RNA Kit (Omega, USA), according to the manufacturer's instructions, and treated with DNase I (Omega, USA) to remove traces of contaminating genomic DNA. RNA purity and concentration were determined using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA quality and integrity was assessed by gel electrophoresis on 1% agarose gels and by Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Only high quality RNA samples were used in subsequent experiments. Each of 1.5 μg RNA per sample was used for the RNA sample preparations. Meanwhile RNA-Seq libraries were constructed using the Illumina NEBNext® Ultra™ RNA Library Prep Kit (NEB, USA). Briefly, mRNA was purified and enriched from total RNA using poly-T oligo-attached magnetic beads. Under elevated temperature, fragmentation was carried out using NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized from the mRNA using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H⁻). Subsequently, second-strand cDNA was synthesized using the first-stand cDNA template, DNA Polymerase I, dNTPs, and RNase H. To each sample, index codes were added to attribute sequences with size-selected and adaptor-ligated cDNA at 37°C for 15 min followed by 95°C for 5 min. Subsequently, PCR amplification was performed using universal PCR primers, index (X) primer, and Phusion High-Fidelity DNA polymerase. Finally, both the purified PCR products and library quality were assessed on the Agilent Bioanalyzer 2100 system. The RNA-Seq libraries were sequenced on Illumina PE150 (Hiseq X Ten) platform at the BGI-Write Genomic Company. RNA-Seq data were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under the accession number SRP154845.

Prior to transcriptome assembly, reads were filtered according to strict criteria. Low quality reads and reads containing poly-N stretches and adapter sequences were removed from raw reads. After filtering, 39–62 million clean reads were obtained per sample (Table S1), which were used for downstream analyses. Additional parameters of clean reads, including GC content, Q20 and Q30 values, and sequence duplication level, were also calculated.

To analyze RNA-Seq data without a reference genome, left read1 files from the nine libraries were pooled into one left.fq file. Similarly, right read2 files from the nine libraries were pooled into the right.fq file. The left.fq and right.fq files were used to de novo assemble the reference transcriptome of A. minutum using Trinity (Grabherr et al., 2011) with min_kmer_cov set to 2.

Gene function annotation was based on the following databases: NCBI non-redundant protein sequences (NR), NCBI non-redundant nucleotide sequences (NT), Protein Family (Pfam), euKaryotic Ortholog Groups (KOG), Swiss-Prot (a manually annotated and reviewed protein sequence database), Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog database (KO), and Gene Ontology (GO). Based on the NR and Pfam databases, Blast2GO v2.5 was used to determine the GO classifications and terms (Götz et al., 2008). KEGG Automatic Annotation Server (KAAS) annotation (http://www.genome.jp/tools/kaas/) was used to obtain KO codes for each unigene. The open reading frames (ORFs) of unigenes were defined according to the results of BLASTX against the NR and Swiss-Prot databases. The coding sequences (CDSs) of unigenes were translated into amino acid sequences, based on genetic codes. Potential ORFs of unigenes with no hits in NR and Swiss-Prot databases were predicted using ESTScan (version 3.0.3). Amino acid sequences were also functionally annotated using the Pfam database (Finn et al., 2008), and functional domains were identified using HMMER 3.0 (http://hmmer.janelia.org/) (Zhang and Wood, 2003).

Transcripts from the assembled reference transcriptome were regarded as reference sequences. Using RSEM (Li and Dewey, 2011), clean reads of all nine libraries were individually mapped onto the reference sequences, and read counts for each gene were calculated. Considering the effect of sequencing depth and gene length on read count, the FPKM parameter, which measures the expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequenced, was calculated from read-count numbers and used for the quantification of gene expression (Trapnell et al., 2010).

DESeq2 (Love et al., 2014) was used to analyze the differential expression of genes between control and nutrient deficient conditions. The resulting p-values were adjusted (p-adj) to control the false discovery rate. Genes with p-adj <0.005 and |log₂fold-change| > 1 were considered as significantly differentially expressed.

To better understand the function of differentially expressed genes (DEGs), GO enrichment analysis and KEGG pathway analysis were performed for DEGs identified in all comparisons, i.e., S6 vs. R, S72 vs. R, and S72 vs. S6). The GOseq R package (Young et al., 2010) based Wallenius non-central hyper-geometric distribution was implemented to test DEGs in GO term functional analysis. KEGG significant enrichment analysis, which provides information on DEGs involved in biochemical metabolic pathways or signal transduction pathways, was performed using KOBAS (Mao et al., 2005).

The results of GO classification revealed genes related to sexual reproduction (GO:0019953, Lev3), developmental process involved in reproduction (GO:0003006, Lev3), and single organism reproductive process (GO:0044702, Lev3). In this study, the description of GO terms (Lev4), including sex determination, sperm-egg recognition, sex differentiation, mating, and fertilization were used for differential gene expression analysis. The FPKM values of all DEGs identified in all three comparisons were used for hierarchical cluster analysis. The heatmap was created using the “pheatmap” package of R.

Five different expression sex-related unigenes (Cluster-20764.3338, Cluster-20764.90078, Cluster-19031.0, Cluster-9927.0 and Cluster-20764.31914) were performed by qRT-PCR to validate the expression of high-throughput sequencing, and primers of these unigenes have been displayed in Table 1. The β-actin gene (internal control) was amplified using β-actin-F and β-actin-R primers (Table 1). The qRT-PCR was conducted on ABI PRISM® 7300 Real Time PCR System using SYBR Premix Ex Taq II (Tli RNaseH Plus) Kit. The PCR was carried out in a total volume of 25 μl, containing 12.5 μl of 2X SYBR® Premix Ex TaqTM II, 2 μl cDNA template, 1 μl each of forward and reverse primer (10 μM each), 1 μl ROX Reference Dye I (50X), and 8.5 μl double distilled water. Cycling conditions were as follows: initial denaturation at 95°C for 30 s, followed by 30 cycles of 95°C for 5 s and 60°C for 31 s. The expression of the five unigenes were normalized relative to that of β-actin. Each sample was performed in triplicate. Dissociation curve analysis was performed after each assay to determine target specificity. The baseline was set automatically by the software, and the expression of sxtA was determined using the 2^−ΔΔCT method (Livak and Schmittgen, 2001).

Although several sxt genes have been identified in A. minutum 454 libraries (Stüken et al., 2011), sequences of the remaining genes related to STX biosynthesis remain inconclusive. The known STX biosynthesis genes of Cylindrospermopsis raciborskii T3 (Kellmann et al., 2008) were used as queries to perform BLAST searches to identify homologs and additional sxt genes in A. minutum and to understand STX biosynthesis (Hackett et al., 2013). Subsequently, the BLAST search strategy was used to identify homologous sxtA genes in the transcriptome of A. minutum based on nucleotide sequences of Alexandrium fundyense (strain CCMP1719) sxtA long (GenBank accession No.JF343239.1) and short isoforms (GenBank accession No.JF343238.1) with an e-value <0.1 (BLASTN). TBLASTN was performed using the SxtA4 protein sequence of A. minutum (GenBank: AIN34681.1) as a query to identify homologous transcripts. Nested primers, including SxtA-longORFF1, SxtA-longORFR1, SxtA-longORFF2, SxtA-longORFR2, SxtA-shortORFF1, and SxtA-shortORFR1, were the used to amplify the full-length SxtA long and short isoform genes. Meanwhile, ten different treatment timepoints were selected to determine the expression profiles of sxtA genes by qRT-PCR. The long and short isoforms of sxtA were amplified using SxtA-LRT-F/SxtA-LRT-R and SxtA-SRT-F/SxtA-SRT-R primer pairs, respectively.

The basic information from RNA-Seq analysis of nine A. minutum samples using Illumina Hiseq X-Ten is summarized in Tables S1–S4. According to GO annotation, 102,495 genes were classified into three functional categories (Figure S1). KOG annotation classified 36,336 unigenes into 25 categories (Figure S2). Additionally, KEGG classification annotated 25,056 unigenes (13.1%) as involved in at least one of the 128 pathways in five major metabolic pathways (human diseases were not included) (Figure S3).

Analysis of gene expression patterns between R and S6 groups using DESeq2 indicated that 41,211 transcripts were differentially expressed (Figure 1A). Of the 41,211 DEGs, 38,719 were significantly up-regulated, and 2,492 were down-regulated in the S6 group compared with R. Comparison between R and S72 groups revealed 20, 280 DEGs, of which 19,473 were significantly up-regulated, and 807 were down-regulated in the S72 group compared with R (Figure 1B). Transcripts were ranked by log₂fold-change values; log₂fold-change > 1 indicated up-regulated transcripts, and log₂fold-change < 0 indicated down-regulated transcripts relative to the control.

Transcriptome analysis of A. minutum revealed that nutrient deficiency for 6 h significantly up-regulated genes involved in toxic substance binding, cell wall macromolecule process, sexual reproduction, germ cell development, P metabolic process, and N compound metabolic process (Figure 2A) and significantly down-regulated genes involved in transferase activity, protein kinase activity, and photosynthetic electron transport chain (Figure 2B). Nutrient deficiency for 72 h led to significantly enriched genes involved in N compound metabolic process, single organism cellular process, chlorophyll binding, and photosynthetic electron transport chain (Figures 2C,D).

Genes encoding protein subunits of photosystem I (PSI) and PSII and those encoding cytochrome (cyt) b6, cyt b6-f complex subunit 4, ferredoxin, and F-type H⁺-transporting ATPase subunits β and ε were significantly down-regulated in the S6 group compared with R. By contrast, psaL (encoding PSI subunit XI), psbO (encoding PSII oxygen-evolving enhancer protein 1), petA (encoding apocytochrome f), and atpE (encoding F-type H⁺-transporting ATPase subunit c) were up-regulated in the S6 group compared with R. Unigenes related to PsaC (PSI subunit VII) and F-type H⁺-transporting ATPase subunit α were both up- and down-regulated. In the S72 group, genes encoding proteins of PSI (PsaA and PsaB), PSII (PsbA, PsbD, PsbC, PsbB, and PsbF), cyt b6-f complex (PetB and PetD), and F-type ATPases (ATPF1B, ATPF1A, and ATPF1E) were significantly down-regulated compared with R (Figure 3).

Compared with the control group (R), the expression of psbF, which is involved in the regulation of KO-K02708 and encodes PSII cyt b559 subunit β, was down-regulated by ~1.1-fold in the S6 group and 7.5-fold in the S72 group. The petD gene, which is involved in the regulation of KO-K02637 and encodes cyt b6f complex subunit 4, was down-regulated by ~10-fold in the S6 group and 3.3-fold in the S72 group compared with R.

Approximately 1,071 unigenes, related to 50 types of proteins, were identified in the endocytosis pathway (Table 2). Among these, unigenes related to clathrin and adaptor protein AP2 were successfully annotated. The unigenes identified in this study represented the majority of genes involved in the endocytosis pathway (Figure 4), especially those encoding clathrin, AP-2, and Hsc70 (Mousavi et al., 2004). Most of the unigenes implicated in clathrin-dependent endocytosis were identified as DEGs under different experimental conditions (Figure S4). A total of 481 DEGs were involved in the endocytosis pathway, of which 470 unigenes, related to 47 types of proteins, were significantly up-regulated in the S6 group compared with R. The remaining three types of proteins (ACP33, VPS25, and CHMP6) were not encoded by DEGs. Compared with R, 5 clathrin-encoding genes, 4 AP-2-encoding genes, and 25 Hsc70 protein-encoding genes were significantly up-regulated, whereas 3 Hsc70 protein-encoding genes were significantly down-regulated in the S6 group. In the S72 vs. R comparison, 294 DEGs involved in endocytosis were identified, of which 288 unigenes were significantly up-regulated in the S72 group compared with R. Detailed information on other related functional genes is summarized in Table S5.

A total of 24 DEGs were annotated to play a role in sex determination (GO:0007530), of which 10 were significantly differentially expressed in the S6 group, and only 3 were significantly up-regulated in the S72 group compared with R (Table 3). In S6 and S72 groups, we selected two unigenes showing significant differential expression compared with the R group; one of these unigenes showed the lowest p-adj value, while the other unigene showed the highest log₂fold-change value. Gene ID, p-adj value, log₂fold-change value, and putative function of both unigenes are listed in Table S6. Sixteen unigenes were annotated as involved in sperm-egg recognition (GO:0035036), of which five unigenes were significantly differential expressed in the S6 group, and only two unigenes were significantly up-regulated in the S72 group compared with R. A total of 126 genes were annotated to play a role in sex differentiation (GO:0007548), of which 10 unigenes were significantly differentially expressed in the S6 group, and only 2 unigenes were significantly up-regulated in the S72 group compared with R. Of the 50 annotated genes in the mating group (GO:0007618), one unigene was significantly down-regulated and 14 unigenes were up-regulated in the S6 group compared with R, whereas only 4 unigenes were up-regulated in the S72 group compared with R. A total of 71 genes were annotated as involved in fertilization (GO:0009566), of which 33 were significantly differentially expressed in the S6 group compared with R. In the S72 vs. S6 comparison, one unigene showed significant down-regulation (Table 3). Cluster analysis was performed to examine the expression pattern of DEGs related to sexual reproduction under different experimental conditions (Figure 5).

Five representative sex-related unigenes were selected and carried out the qRT-PCR analyses Compared with the control groups, the expression levels of all the selected unigenes under treated conditions were significantly up- or down-regulated (P < 0.01; Figure 6), which confirmed the authenticity of high-throughput sequencing results.

To identify DEGs involved in cell wall biogenesis, we analyzed unigenes generated from R, S6, and S72 groups based on the comparison of |log₂fold-change| values. Among the 218 unigenes (Table 3) mentioned above, we selected unigenes with p-adj <0.05 from R vs. S6, R vs. S72, and S6 vs. S72 comparisons. In R vs. S6 comparison, only one unigene was significantly down-regulated and 53 unigenes were up-regulated in the S6 group. In R vs. S72 comparison, 31 unigenes were up-regulated in the S72 group. In S6 vs. S72 comparison, only one unigene was up-regulated and one gene was down-regulated (Table 3). These data demonstrate that ~24.3% of the unigenes were up-regulated and 0.46% were down-regulated upon stress treatment for 6 h, and 14.2% of the unigenes were up-regulated upon stimulation with the deficiency of N and P for 72 h.

To identify genes potentially involved in STX biosynthesis, 34 annotated Sxt peptides from the Cylindrospermopsis raciborskii T3 STX biosynthesis gene cluster were queried against A. minutum transcriptome assembled in this study using tblastn. A total of 489 homologs of 21 proteins were obtained (Table 4, and Table S7), which represented 21 of the 34 genes in C. raciborskii. Some of the sxt genes directly involved in STX biosynthesis were identified, including sxtA, sxtD, sxtG, and sxtS. Approximately 50% of the 21 genes were significantly up-regulated in the S6 group compared with R, while only 4 sxt genes were significantly up-regulated in the S72 group (Table S7). Furthermore, three putative homologs of sxtV, related to dioxygenase reductase, were first identified in dinoflagellates (Kellmann et al., 2008, 2010; Hackett et al., 2013). Based on the proposed pathway of STX biosynthesis, as previously reported, the putative STX biosynthesis gene homologs were presented in a revised STX biosynthesis pathway (Figure 7) (Stüken et al., 2011).

Next, we used two types of sxtA genes of A. fundyense (encoding the long and short isoforms) to identify sxtA gene homologs in A. minutum. Unigenes identified using blastn were subsequently searched in the NT and NR databases for gene annotation. Nine transcripts were annotated as sxtA long isoform, and one transcript was annotated as sxtA short isoform (Table S8). Additionally, three homologous transcripts were annotated as sxtA4 (Table S8). To examine the temporal expression profile of sxtA genes under nutrient deficient conditions, we performed qRT-PCR (Figure 8). The expression of sxtA short isoform and sxtA long isoform genes was down-regulated at 30 min after exposure to nutrient deficient conditions and decreased further over the next few hours. The expression of sxtA genes was the lowest at 6 h post treatment and then continued to increase until 24 h. Subsequently, the expression of sxtA genes declined until 48 h and then was slightly up-regulated at 72 h.

N is a major limiting factor affecting cellular biosynthesis and growth in dinoflagellates. Change in N concentration can greatly affect primary and secondary metabolism in dinoflagellates (Dagenais-Bellefeuille and Morse, 2013). N limitation causes a reduction in protein and chlorophyll levels (Lei and Song, 2011; Hockin et al., 2012; Pleissner and Eriksen, 2012; Zhang et al., 2019), resulting in the down-regulation of photosynthetic pathways, which may contribute to the reduced N demand of cells. P is one of the main limiting factors affecting algae growth. P is a component of a large number of molecules necessary for phytoplankton growth (such as DNA, RNA, and lipids), intracellular energy (e.g., ATP), and signaling pathways (Lomas et al., 2010).

In this study, genes involved in photosynthesis, particularly PSI and PSII, were significantly down-regulated under N and P deficient conditions. In conjunction with photosystem components, genes involved in light-harvesting complex II (LHCII; lhcb1 and lhcb2) and ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) were down-regulated under nutrient deficient conditions. This is consistent with the response of other algae to N deficiency. For example, the lhc genes showed 4-fold reduction in expression in Emiliana huxleyi under N-limitation compared with the control (Dyhrman and Anderson, 2003). In Prymnesium parvum, a toxin-producing microalga, N and P starvation down-regulates the expression of cyt genes (Beszteri et al., 2012). Transcriptome analysis of A. catenella (ACHK-T) and its non-toxic mutant form (ACHK-NT) shows the same results as our study, with up-regulated processes such as photosynthesis in ACHK-NT (Zhang et al., 2014b). Transcriptome analysis of S. trochoidea under N limiting condition revealed significantly inhibited photosynthesis (Cooper et al., 2016), which is consistent with our study. Reducing photosynthesis is perhaps a way to compensate for increases in cellular carbon: nitrogen (C:N) ratios; for example, in S. trochoidea and A. fundyense, the ratio of particulate organic C (POC) to PON increases under N limitation (Eberlein et al., 2016). Apoptosis of chloroplasts is a feature of resting cysts, which corresponds to the down-regulation of photosynthesis.

Autotrophy, mixotrophy, and heterotrophy are the main trophic modes of marine dinoflagellates. Like the photosynthetic species, heterotrophic dinoflagellates also contribute significantly to the ecosystem (Tillmann and Hesse, 1998; Yang et al., 2004). Various feeding mechanisms are observed in mixotrophic or heterotrophic dinoflagellates, such as diffusion feeding and filter feeding (Jeong et al., 2005b, 2010). In our study, KEGG enrichment pathway showed that genes involved in the endocytosis pathway were significantly up-regulated. We identified 190 unigenes related to 45 types of proteins that participate in the endocytosis pathway (Table 2). An intracellular vesicle is formed when a substance enters a cell. With prolonged exposure to nutrient deficient conditions (72 h), unigenes encoding clathrin, AP-2, and Hsc70 proteins were significantly up-regulated. The results indicate that the entry of nutrients in A. minutum cells mainly via clathrin-dependent endocytosis. Similar results have been reported in the transcriptome analysis of A. catenella; 131 unigenes related to 23 types of proteins involved in clathrin-dependent endocytosis were observed (Zhang et al., 2014a). This implies that A. minutum has a tendency to exhibit mixotrophy under nutrient deficient conditions. Thus, it is possible that clathrin-dependent endocytosis is critical for the maintenance of heterotrophy in A. minutum.

In A. minutum, resting cysts are formed as a result of sexual reproduction involving the formation of gametes, followed by their fusion. Some dinoflagellate species exhibit facultative sex (Pfiester, 1989; Kremp, 2013; Figueroa et al., 2015). Cyst formation is a positive response to environmental stress such as genetic recombination, expansion of geographical distribution, and seedlings of the recurrent blooming (Dale, 1983). Cyst formation is one of the best ways for dinoflagellates to survive unfavorable conditions.

GO enrichment analysis of our transcriptome data revealed some GO terms related to sexual reproduction. During the first 6 h under N and P deficient conditions, the most significant change in expression was observed in genes encoding sex determination factors, which determine the sex of the gamete. Compared with the R group, Cluster-18154.1 showed the lowest p-adj value and the highest log₂fold-change in the S6 group. Based on the NR database, Cluster-18154.1 was annotated as ubiquitin carboxyl-terminal hydrolase 22. Ubiquitin carboxyl-terminal hydrolase enzymes hydrolyze thiol esters formed between small thiols (e.g., glutathione) and the ubiquitin carboxyl terminus (Rose and Warms, 1983). These enzymes also hydrolyze glycine methyl ester and spermidine (Pickart and Rose, 1985). Cluster-18154.1 may play an important role in G protein-coupled receptor signaling pathway. Vegetative cells differentiate into male and female gametes through sperm-egg recognition. The same phenomenon is observed in the life cycle of Chlamydomonas reinhardtii, where vegetative cells differentiate into two types of gametes (mt⁺ and mt⁻) during N starvation (Sekimoto, 2017).

Mating and fertilization are important processes during sexual reproduction. Planktonic zygotes are formed by the fusion and mating of haploid gametes. Subsequently, diploid zygotes swim for several days before transforming into resting cysts. In this study, at 72 h under N and P starvation, 71 unigenes were annotated to play a role in fertilization, and only 1 unigene (Cluster-20764.19572) was significantly down-regulated in S72 compared with S6. The most likely function of is the fusion of sperm to the egg plasma membrane. The annotation from Pfam database indicated that it was an acrosome formation-associated factor (Afaf). Acrosome is a specialized apical vesicle of spermatozoa derived from the Golgi apparatus. Acrosome overlies the nucleus of mature spermatozoon and is formed in response to fertilization. Afaf antibodies reduce the rate of in vitro fertilization by significantly inhibiting the penetration of eggs by sperm (Hu et al., 2010). Therefore, down-regulation of Cluster-20764.19572 possibly enhanced fertilization. Thus, expression profiles of the abovementioned unigenes were significantly different between the control (R) and S6 groups, suggesting that A. minutum may carried out the ecophysiological adaption under N and P limiting conditions. With the treatment time last for 72 h, we found the expression quantity of most unigene had returned to normal level. These results revealed that nutritional deficiency might not be a strong stimulating factor for a long time. Our results were similar with other studies, i.e., nutrient deficiency alone is not enough to induce sexuality either in nature (Garcés et al., 2004; Kremp et al., 2009; Brosnahan et al., 2017) or in culture (Kremp et al., 2009; Figueroa et al., 2011).

Additionally, a thicker protective cell wall was observed during the sexual reproduction process. This suggests that genes involved in cell wall biogenesis are up-regulated in A. minutum to resist the adverse environment. However, it should be noted that under prolonged exposure to nutrient deficient conditions for 72 h, the expression level of most unigenes returned to the level observed at 0 h. These results suggest that nutrient deficiency does not act as an inducer of sexual reproduction over a long term.

Comparison between toxin-producing and non-toxic dinoflagellates shows that sxtA4 is critical for toxin biosynthesis (Stüken et al., 2011). Additionally, phylogenetic analysis indicates that the C-terminal end of sxtA and sxtG plays important roles in STX biosynthesis (Hackett et al., 2013). In this study, we not only found the SXT genes related to toxin biosynthesis but also identified the sxtV gene. This is the first report of the identification of sxtV in the Alexandrium genus. The sxtV gene putatively encodes dioxygenase reductase, which belongs to the succinate dehydrogenase family. We also identified and characterized the long and short isoforms of sxtA gene. In the majority of toxin-producing dinoflagellates, sxtA produces two types of transcripts: the short transcript contains three domains (sxtA1–3), while the long transcript harbors four domains (sxtA1–4) (Stüken et al., 2015).

The long isoform of sxtA initiates toxin biosynthesis in cyanobacteria, indicating the long sxtA transcript in A. minutum may be directly involved in toxin biosynthesis (Kellmann and Neilan, 2007). Based on the FPKM values, homologous sxtA4 transcripts were expressed in all experiment groups (R, S6, and S72). This suggests that the strain of A. minutum used for RNA-Seq analysis and transcriptome assembly may be a toxin-producing strain. Previously, analysis of two kinds of Alexandrium ostenfeldii (PSP toxin-producing and non-producing) strains showed that sxtA4 and sxtA1 fragments were not present in the non-PSP-producing strains, whereas the PSP-producing strains contained both sxtA1 and sxtA4 fragments (Suikkanen et al., 2013). Additionally, neither PSP nor sxtA4 gene fragments were detected in the non-toxic A. minutum strain, VGO663 (Stüken et al., 2015). Compared with the non-PSP-producing S. trochoidea strain, neither sxtA1 nor sxtA4 gene fragments were detected in its transcriptome (Cooper et al., 2016). We speculate that the long isoform of sxtA performs the first step in toxin production. Both sxtA genes were down-regulated in S6 and S72 groups compared with the R group, suggesting that A. minutum produces less toxins under nutrient deficient conditions. These results are consistent with those of previous studies (Wang and Hsieh, 2002; Dyhrman and Anderson, 2003; Leong et al., 2004; Hu et al., 2006; Vanucci et al., 2012). The toxin contents of Alexandrium and Ostreopsis ovata decrease under N limiting conditions, possibly because algae divert more energy toward sexual reproduction under adverse environments, in order to propagate, than toward toxin production (Vanucci et al., 2012; Zhang et al., 2014b).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.