Biecheleriopsis adriatica was isolated from plankton samples collected from surface waters off the coast of Tongyoung, southern Korea, using plankton samplers in August 2006, when the water temperature and salinity were 28.0°C and 31.0, respectively (Jang et al., 2015). The collected samples were filtered using a 154-μm Nitex mesh, and a clonal culture of B. adriatica BATY06 was established using two consecutive single-cell isolations (Jang et al., 2015). When B. adriatica concentration increased sufficiently, it was transferred to 32-, 270-, and 500-mL polycarbonate bottles containing f/2-Si medium (Guillard and Ryther, 1962). The bottles were placed on a shelf at 20°C for the cold stress experiment and 25°C for the optimal temperature and heat stress experiments under a 14 h:10 h light:dark cycle and illuminated with 50 μmol photons m^–2 s^–1 using cool white fluorescent light.

To determine the optimal, lowest, and highest water temperatures for survival, the growth rates of B. adriatica were measured at 14–34°C. The growth rates were almost zero at 14 and 32°C; thus, these two temperatures were used for cold and heat stress. The optimal temperature for maximum growth rate was 25°C.

For cold stress, a dense B. adriatica culture (approximately 100,000 cells mL^–1) grown autotrophically at 20°C was distributed to three 800-mL tissue culture flasks (Falcon; Heidelberg, Germany) containing fresh f/2-Si medium. For heat stress and optimal temperature, a dense B. adriatica culture (approximately 100,000 cells mL^–1) grown autotrophically at 25°C was distributed to six 800-mL tissue culture flasks containing fresh f/2-Si medium.

For RNA sequencing, the triplicate flasks for each target temperature were incubated in three temperature-controlled chambers set up at 14, 25, and 32°C for 11 days after gradual acclimation (Figure 1A). The temperature in each chamber was monitored using HOBO data loggers (HOBO MX2202 Pendant MX Temperature/Light Data Logger; Onset; Bourne, MA, United States).

To monitor B. adriatica cell density at each target temperature, a 5-mL aliquot was collected from each flask and fixed with 5% acidic Lugol’s solution daily. Aliquots were collected from the fixed samples, and ≥ 200 B. adriatica cells were enumerated using a Sedgewick-Rafter counting chamber to determine B. adriatica density. After incubating for 6 days at 25°C, B. adriatica reached the stationary phase and were transferred into new 800-mL flasks containing fresh f/2-Si medium.

To confirm that nutrient concentration did not affect B. adriatica growth, a 10-mL aliquot was collected from each flask on days 0, 6, and 11, gently filtered through GF/F filters (Whatman Inc., Floreham Park, NJ, United States), and stored at –20°C until nutrient concentration measurement. The concentrations of nitrate plus nitrite (hereafter NO₃), phosphate (PO₄), and silicate (SiO₂) were measured using a nutrient auto-analyzer system (QuAAtro; Seal Analytical GmbH; Nordestedt, Germany). The samples diluted to 1/10 were used for the measurement to consider the limit concentration of the auto-analyzer system, and it was considered when the final nutrient concentration was calculated after the measurement.

To harvest the similar cell number of B. adriatica (1 × 10⁷ cells) after incubating for 11 days, 600, 200, and 300 mL aliquots were collected from each flask at 14, 25, and 32°C, respectively. The cells in the aliquots were collected by centrifugation at 2,808× g for 5 min (Labogene 1696R; Gyrozen Co., Gimpo, Korea), immediately frozen in liquid nitrogen, and stored at –80°C until RNA extraction.

Total RNA was extracted individually from each sample using the RNeasy Plant Mini Kit (Cat. No. 74903; Qiagen; Hilden, Germany) and treated with RNase-Free DNase set (Cat. No. 79254; Qiagen) to eliminate residual genomic DNA contaminants. The extracted RNA samples were purified using a RNeasy Kit (Cat. No. 74104; Qiagen), and the total RNA amount was assessed using Quant-IT RiboGreen (Cat. No. R11490; Invitrogen; Carlsbad, CA, United States). The integrity of the total RNA was evaluated using the TapeStation RNA Screen Tape (Cat. No. 5067-5576; Agilent Technologies; Santa Clara, CA, United States), and only RNA preparations with RNA integrity number (RIN) > 7.0 were used for RNA library construction.

The nine cDNA libraries were separately constructed using an Illumina TruSeq Stranded mRNA LT Sample Prep Kit (San Diego, CA, United States), and paired-end (2 × 100 bp) RNA of the pooled libraries was sequenced on an Illumina NovaSeq 6000 at Macrogen (Seoul, Korea). In the first step of this process, only eukaryotic polyadenylated mRNA was filtered to eliminate potential bacterial contaminants from the total RNA of B. adriatica.

The quality, total bases, total reads, GC content (%), and basic statistics of the sequenced raw reads of the nine libraries were calculated, and the raw reads of low-quality, adapter, contaminant DNA, or PCR duplicates were removed using Trimmomatic program (Bolger et al., 2014). The quality of the trimmed reads (clean reads) was checked using FastQC v0.11.7² before assembly.

The clean reads for the nine libraries were assembled de novo using Trinity with the “–SS_lib_type FR” option, which is typically used for de novo transcriptome reconstruction (Grabherr et al., 2011). It combines the read sequences of a certain length of overlap to form long fragments without N gaps, called contigs. Only contigs longer than 200 bp were selected. The longest assembled contigs were filtered, and redundant contigs were removed (i.e., unigenes) using CD-HIT-EST, with a clustering identity of 0.98 (Li and Godzik, 2006).

For functional annotation, the unigenes were searched against NCBI non-redundant (NCBI nr) and Gene Ontology (GO) using BLASTx of DIAMOND software with an E-value default cutoff of 1.0E–5 (Buchfink et al., 2015).

Unigene expression levels were estimated using RSEM v1.3.1 (Li and Dewey, 2011). The clean reads from each sample were mapped back onto the assembled unigenes, and the read count value of each unigene was obtained from the mapping results. The unigene expression level (log₂ fold change) between the samples under cold stress and optimal temperature or under heat stress and optimal temperature was calculated as fragments per kilobase of unigene per million mapped reads (FPKM). If more than one FPKM value was 0, it was excluded from the analysis to avoid reducing the power of the statistical hypothesis testing. After the addition of + 1 to each of the FPKM values from filtered unigenes, they were transformed to log₂ scale and subjected to relative log expression (RLE) normalization. Statistically significant results of differential gene expression were obtained with |log₂ fold change| ≥ 1 and p < 0.05, using the nbinomWaldTest of DESeq2 R library (Anders and Huber, 2010). Hierarchical clustering analysis of differentially expressed genes (DEGs) was performed using the Euclidean method and complete linkage, representing the similarity of expression patterns between unigenes and samples.

Based on the GO annotation results, functional unigene classification and enrichment were analyzed using Web Gene Ontology Annotation Plot (WEGO)³. The Pearson’s chi-square test was automatically carried out using WEGO to test whether the GO categories were significantly differentiated between the upregulated and downregulated gene datasets (Ye et al., 2006). Differences in GO categories between the two datasets were considered significant at a p-value < 0.05. In addition, the number of DEGs obtained with |log₂ fold change| ≥ 1 and adjusted p-value (padj) < 0.05 was also determined (Benjamini and Hochberg, 1995).

The unigenes were classified into nine types based on the trends in expression under cold and heat stress against that under optimal temperature (Figure 2). Type 1 was composed of DEGs upregulated under both cold and heat stress; Type 2 consisted of unigenes upregulated under cold stress, but not significantly differentially expressed between optimal temperature and heat stress; Type 3 was composed of DEGs upregulated under cold stress but downregulated under heat stress; Type 4 was composed of unigenes not significantly differentially expressed between cold stress and optimal temperature but upregulated under heat stress; Type 5 consisted of unigenes not significantly differentially expressed between cold stress and optimal temperature and also between heat stress and optimal temperature; Type 6 consisted of unigenes not significantly differentially expressed between cold stress and optimal temperature but downregulated under heat stress; Type 7 was composed of DEGs downregulated under cold stress but upregulated under heat stress; Type 8 was composed of unigenes downregulated under cold stress but not significantly differentially expressed between optimal temperature and heat stress; Type 9 consisted of DEGs downregulated under both cold and heat stress. The number of unigenes in each type and the proportion of each type to all unigenes or DEGs was determined.

To investigate the functions of representative DEGs belonging to each type, the unigenes annotated against eukaryotic organisms in the NCBI nr database with E-value cutoff E–30 were ranked by their expression level, and the functions of the high-ranking unigenes were examined. The hypothetical or predicted proteins were excluded from the identification and further screened using an NCBI conserved domain search. Unigenes with functions revealed previously were obtained. The rank of the unigenes was determined as follows: In the first step, the fold change values of the unigene expression under cold stress were arranged in the descending or ascending orders if they were upregulated or downregulated, respectively. In the second step, the unigenes were ranked. In the third step, the first and second steps were repeated for the fold change values of the same unigenes under heat stress. In the fourth step, the two rank numbers assigned to each unigene under cold and heat stress were added, and the unigenes were arranged according to the ascending order of the added number. If there was no differential expression of unigenes under cold or heat stress depending on the type, they were not calculated. Therefore, Type 5 was excluded from this calculation and examination of the function of the unigenes.

Under cold stress, the average B. adriatica cell density (± standard error) decreased slowly from 24,500 (± 1,030) cells mL^–1 on day 0 to 13,330 (± 390) cells mL^–1 on day 11, when they were harvested for RNA sequencing (Figure 1B). Under heat stress, the average B. adriatica cell density increased from 21,500 (± 1,410) cells mL^–1 on day 0 to 54,400 (± 2,390) cells mL^–1 on day 3 and then decreased slowly to 28,900 (± 4,320) cells mL^–1 on day 11 (Figure 1B). Under optimal temperature, the average B. adriatica cell density increased from 20,500 (± 920) cells mL^–1 on day 0 to 138,700 (± 3,180) cells mL^–1 on day 6, entering the stationary phase. Thus, these cultures were transferred into new 800 mL flasks at an initial density of 4,200 (± 140) cells mL^–1, and the B. adriatica cell density increased to 64,700 (± 770) cells mL^–1 on day 11 (Figure 1B).

Under cold stress, the average (± standard error) NO₃ concentration at day 0 was 882.8 (± 3.8) μM and remained unchanged on days 6 and 11 (Figure 1C). Similarly, the average PO₄ and SiO₂ concentration at day 0 were 26.9 (± 0.4) and 60.9 (± 0.5) μM, respectively, and remained unchanged on days 6 and 11 (Figure 1D).

The average NO₃, PO₄, and SiO₂ concentrations on days 0, 6, and 11 under optimal temperature and heat stress were similar to those under cold stress (Figures 1C,D).

Illumina sequencing of B. adriatica generated 130–202 million paired-end 100-bp clean reads per library (nine libraries = three triplicates at each temperature condition), with a trimming process for removing adapter sequences and poor-quality bases (Table 1). The GC content of the clean reads was between 54.30 and 54.92% and the ratio of bases with Phred quality scores ≥ 20 and 30 (Q20 and Q30) of the clean reads exceeded 99.14 and 96.70%, respectively (Table 1). A total of 198,703 unigenes were generated by a combined de novo assembly with an average 850 bp length and a 1,562 bp N50 (Table 1). Among the 198,703 unigenes, 35,370 (17.8%) and 65,580 (33.0%) were functionally annotated against the GO and NCBI nr databases, respectively, with an E-value default cutoff of 1.0E – 5.

The number of DEGs in B. adriatica under cold and heat stress were 24,916 and 13,366, respectively, at the expression level of |log₂ fold change| ≥ 1 and p < 0.05 (Figure 3). Under cold stress, 13,159 and 11,757 unigenes were upregulated and downregulated, respectively. In comparison, 6,345 upregulated and 7,021 downregulated unigenes were identified under heat stress (Figure 3). At the expression level of |log₂ fold change| ≥ 1 and padj < 0.05, the number of DEGs in B. adriatica under cold and heat stress were 23,965 and 12,185, respectively (Supplementary Figure 1), and 12,812 and 11,153 genes were upregulated and downregulated, respectively, whereas 5,892 upregulated and 6,293 downregulated genes were identified under heat stress.

The heatmap for the DEGs depicting the result of hierarchical clustering analysis showed more similar patterns of DEGs between heat stress and optimal temperature than those between cold stress and optimal temperature (Figure 4).

The biological function of the identified DEGs of B. adriatica under cold stress was investigated using GO analysis. The results of GO functional classification and enrichment were visualized and compared using WEGO (Figure 5).

Under cold stress, 4,370 upregulated and 2,866 downregulated DEGs were annotated to the GO database (29.0% DEGs under cold stress), which were classified into 21, 18, and 37 GO subcategories in cellular component (CC), molecular function (MF), and biological process (BP) categories, respectively, in GO level 2 (Figure 5 and Supplementary Table 1).

In the CC category, “cell” and “cell part” subcategories possessed most DEGs, followed by “organelle” and “membrane” (Figure 5A and Supplementary Table 1). Six of the 21 subcategories had significantly different percentages of upregulated and downregulated DEGs (Pearson’s chi-squared test; p < 0.05) (Figure 5A, Table 2, and Supplementary Table 1). The percentage of upregulated DEGs was higher than that of downregulated DEGs in “membrane-enclosed lumen” and “organelle part” subcategories (Pearson’s chi-squared test; p = 0.001, p = 0.001), whereas the percentage of downregulated DEGs was higher than that of upregulated DEGs in “supramolecular complex,” “membrane,” “synapse,” and “cell junction” subcategories (Pearson’s chi-squared test; p < 0.001, p = 0.001, p = 0.013, p = 0.026; Figure 5A, Table 2, and Supplementary Table 1).

In the MF category, most DEGs belonged to “catalytic activity” and “binding” subcategories (Figure 5B and Supplementary Table 1). Four of the 18 subcategories had significantly different percentages of upregulated and downregulated DEGs (Pearson’s chi-squared test; p < 0.05; Figure 5B, Table 2, and Supplementary Table 1). The percentage of upregulated DEGs was higher than that of downregulated DEGs in the “structural molecule activity” subcategory (Pearson chi-squared test; p = 0.001), whereas the percentage of downregulated DEGs was higher than that of upregulated DEGs in “transporter activity,” “signal transducer activity,” and “binding” subcategories (Pearson’s chi-squared test; p < 0.001, p < 0.001, p < 0.001; Figure 5B, Table 2 and Supplementary Table 1).

In the BP category, “cellular process” subcategory had the highest percentage of DEGs, followed by “metabolic process” and “biological regulation” subcategories (Figure 5C and Supplementary Table 1). Eleven of the 37 subcategories had significantly different percentages of upregulated and downregulated DEGs (Pearson’s chi-squared test; p < 0.05; Figure 5C, Table 2, and Supplementary Table 1). Among them, the percentage of upregulated DEGs in “metabolic process” subcategory was higher than that of downregulated DEGs (Pearson’s chi-squared test; p = 0.001; Figure 5C, Table 2, and Supplementary Table 1).

At the expression level of |log₂ fold change| ≥ 1 and padj < 0.05, there was no change in the trend of the results of GO analysis at level 2 under cold stress (Supplementary Figure 2).

In GO level 3, 77 subcategories had significant differences between the percentages of the upregulated and downregulated genes under cold stress (Pearson’s chi-squared test; p < 0.05; Supplementary Table 2). Among them, 23, 13, and 41 subcategories belonged to the CC, MF, and BP categories, respectively (Supplementary Table 2). Furthermore, the percentage of upregulated DEGs of 17 subcategories was higher than that of downregulated DEGs, whereas the remaining 60 subcategories had a higher percentage of downregulated DEGs than that of upregulated DEGs (Supplementary Table 2). The top 5 subcategories in GO level 3 with the largest differences between the percentage of the upregulated and downregulated genes significantly in B. adriatica under cold stress, based on the p-value, were “methylation,” “ion binding,” “cell communication,” “transmembrane transporter activity,” and “catalytic activity, acting on RNA” in sequence (Table 3).

When the biological function of the identified B. adriatica DEGs under heat stress was investigated, 1,451 upregulated and 1,504 downregulated genes were annotated to the GO database (22.1% DEGs under heat stress), which were classified into 21, 12, and 33 GO subcategories in the CC, MF, and BP categories, respectively, at GO level 2 (Figure 6).

In the CC category, “cell,” “cell part,” “organelle,” and “membrane” subcategories were dominant (Figure 6A and Supplementary Table 3). There were no subcategories with significant differences between the percentages of upregulated and downregulated genes (Figure 6A, Table 2, and Supplementary Table 3).

In the MF category, “catalytic activity” subcategory had the highest percentage of DEGs, followed by “binding” subcategory (Figure 6B and Supplementary Table 3). The “catalytic activity” subcategory had a higher percentage of downregulated genes than that of upregulated genes (Pearson’s chi-squared test; p < 0.001), but the percentage of downregulated genes in other subcategories was not significantly different from that of the upregulated genes (Figure 6B, Table 2, and Supplementary Table 3).

In the BP category, most DEGs were involved in the “cellular process” and “metabolic process” (Figure 6C and Supplementary Table 3). Four of the 33 subcategories had significantly different percentages of upregulated and downregulated DEGs (Pearson’s chi-squared test; p < 0.05; Figure 6C, Table 2, and Supplementary Table 3). The percentage of the downregulated genes in “nitrogen utilization,” “response to stimulus,” “metabolic process,” and “cellular process” subcategories were significantly higher than that of the upregulated genes (Pearson’s chi-squared test; p = 0.005, p = 0.006, p = 0.012, p = 0.016), whereas there were no significant differences between the percentage of the upregulated and downregulated genes in other subcategories (Figure 6C, Table 2, and Supplementary Table 3).

At the expression level of |log₂ fold change| ≥ 1 and padj < 0.05, there was no change in the trend of the results of GO analysis under heat stress (Supplementary Figure 3).

In GO level 3, 29 subcategories had significant differences between the percentages of the upregulated and downregulated genes (Pearson’s chi-squared test; p < 0.05) under heat stress, with 9, 7, and 13 subcategories belonging to the CC, MF, and BP categories, respectively (Supplementary Table 4). Most subcategories had a higher percentage of downregulated genes than that of the upregulated genes (24 subcategories; Supplementary Table 4). The top 5 subcategories in GO level 3 with the largest differences between the percentage of the significantly upregulated and downregulated genes in B. adriatica under heat stress, based on the p-value, were “regulation of biological quality,” “small molecule metabolic process,” “response to abiotic stimulus,” “carbohydrate derivative binding,” and “ion binding” in sequence (Table 3).

When the unigenes were classified into nine types according to the expression trend, Type 5 had the greatest percentage of the number of the genes relative to that of all unigenes (83.4%; Figure 7A), followed by Type 2 (5.2%), Type 8 (4.8%), Type 6 (2.3%), Type 4 (1.8%), Type 1 (1.2%), Type 9 (1.0%), Type 3 (0.2%), and Type 7 (0.1%) (Figure 7A). Furthermore, Type 2 had the greatest percentage of the number of genes relative to that of the DEGs (31.0%; Figure 7B), followed by Type 8 (28.6%), Type 6 (13.6%), Type 4 (11.0%), Type 1 (7.4%), Type 9 (6.2%), Type 3 (1.4%), and Type 7 (0.8%) (Figure 7B). At the expression level of |log₂ fold change| ≥ 1 and padj < 0.05, there was no change in this trend (Supplementary Figure 4).

Among the unigenes ranked high by expression level, unigenes related to anti-stress, metabolic pathway, photosynthesis, and cell motility were detected in each type (Table 4). In Type 1, unigenes encoding “apolipoprotein d lipocalin” and “alpha-protein kinase vwkA” were detected; in Type 2, unigenes encoding “facilitated trehalose transporter Tret” and “glutathione S-transferase class-mu 28 kDa isozyme” were detected; in Type 3, unigenes encoding “fatty acid desaturase” and “fructose-bisphosphate aldolase” were detected; in Type 4, unigenes encoding “FAD-binding domain-containing protein” and “malate l-lactate dehydrogenase” were detected; in Type 6, unigenes encoding “proline/betaine transporter” and “aldo/keto reductase” were detected; in Type 7, unigenes encoding “chloroplast light harvesting complex protein, partial” and “intraflagellar transport protein 46-like” were detected; in Type 8, unigenes encoding “DnaJ-like subfamily B member 6-A” and “copia protein” were detected; in Type 9, unigenes encoding “troponin C, skeletal muscle” and “palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase, chloroplastic” were detected (Table 4).

Some studies determined the number of DEGs of dinoflagellates at the expression level of p < 0.05 (e.g., Guo et al., 2016; Jang et al., 2019), whereas others determined that at the expression level of padj < 0.05 (e.g., Barshis et al., 2014; Bellantuono et al., 2019; Lin et al., 2019). The number of DEGs in B. adriatica at the expression level of p < 0.05 and padj < 0.05, were 33,049 and 30,877, respectively, indicating a 6.6% difference. The results of GO analysis and expression type at the expression level of p < 0.05, and padj < 0.05, did not change significantly. Furthermore, there was no change in the trends of any of the results.

At the expression level of p < 0.05, the larger number of upregulated DEGs than downregulated DEGs under cold stress indicates that B. adriatica expresses more genes in response to cold stress. In contrast, the higher number of downregulated DEGs than upregulated DEGs under heat stress suggests that B. adriatica restricts gene expression to save energy or that genes are damaged by heat shock. At the geological scale, the sea surface temperature during the mid- and late-Pleistocene (∼1.2 Ma) showed a cooling trend, although a glacial–interglacial cycle with large amplitudes of air temperature change was repeated (McClymont et al., 2013). The period of global warming in the last century was significantly shorter than that in the ice age. Thus, dinoflagellates may have acclimated to cold stress better than to heat stress. Furthermore, dinoflagellates conduct diurnal vertical migration, and some dinoflagellates pass through the thermocline and remain in cold water, although the surface water temperature increases significantly (Jeong et al., 2015). However, vertical migration may not provide a refuge when the surface water temperature decreases significantly. Thus, dinoflagellates may upregulate or downregulate gene expression under cold stress than under heat stress.

The “structural molecule activity,” “membrane-enclosed lumen,” and “organelle part” subcategories in GO level 2 are related to cell structure and organelles. Cold stress may damage the protist cell structure. The morphology of the green alga Chlamydomonas reinhardtii changed significantly and its organelles were progressively disorganized and disrupted after decreasing the water temperature from 25 to 5°C (Valledor et al., 2013). After 24 h of water temperature decrease, the nucleolus of C. reinhardtii could not be distinguished; vacuoles and starch sheaths increased. Furthermore, the chloroplast changed its shape and flagella were shortened and disappeared after 120 h (Valledor et al., 2013). Similarly, under cold stress, B. adriatica might upregulate genes related to maintaining cell integrity instead of genes with other functions.

Under stress, microalgae tend to reserve extra metabolites by promoting metabolic processes or save energy by reducing metabolic processes (Levin et al., 2016; Sun et al., 2018). In this study, the “metabolic process” subcategory had a significantly higher percentage of upregulated DEGs than that of downregulated DEGs under cold stress. Therefore, B. adriatica may respond to cold stress by enhancing metabolic processes to reserve extra metabolites.

Under cold stress, the “signaling” subcategory had the greatest difference between the percentages of upregulated and downregulated DEGs, that is, the smallest p-value. This was followed by “response to stimulus” and “biological regulation.” The basic activities of B. adriatica, which perceive environmental change, transmit signals, and respond to changes by modulating biological processes or functions, seemed to be disturbed by cold stress.

Under cold stress, “methylation” had the highest difference between upregulated and downregulated genes in GO level 3. In organisms, methylation plays diverse roles, such as heavy metal modification, gene expression and protein function regulation, and RNA processing (Lee et al., 2010; Li et al., 2018). DNA methylation, which adds a methyl group to the DNA, can affect the gene expression of dinoflagellates, and temperature change can regulate DNA methylation levels (Bayer et al., 2012; Naydenov et al., 2015; Li et al., 2018). Therefore, cold stress may affect DNA methylation and gene expression in B. adriatica.

The number of subcategories showing a significant difference between the percentages of upregulated and downregulated DEGs under heat stress was lower than that under cold stress. The catalytic activity generally increases with temperature, but it decreases sharply due to enzyme denaturation if the temperature exceeds a certain threshold (Daniel and Danson, 2013). Therefore, enzyme denaturation in B. adriatica cells may occur at 32°C. Furthermore, nitrogen utilization efficiency decreases at non-optimal temperatures owing to the changes in cytoplasmic viscosity. Heat stress may also affect the cytoplasmic viscosity and nitrogen utilization efficiency in B. adriatica (Raven and Geider, 1988). Levin et al. (2016) reported that the downregulated DEGs of the thermo-tolerant dinoflagellate Cladocopium sp. under heat stress were distributed in more GO categories related to the metabolic and biosynthetic processes than those of thermo-sensitive Cladocopium sp. In some organisms, reduced metabolic and biosynthetic activity under stress is correlated with the increased survival time of organisms while saving energy (Hand and Hardewig, 1996). Therefore, B. adriatica may try to acclimate to high water temperatures with reduced energy consumption.

Under heat stress, only “regulation of biological quality” had a higher percentage of the upregulated DEGs than that of the downregulated DEGs among the top 5 subcategories. It includes many processes that modulate a qualitative or quantitative trait with a measurable attribute in an organism or part of an organism, such as size, mass, shape, and color (Gene Ontology). It was found that increasing the water temperature by 1°C caused an approximately 2.5% intra-specific decrease in protist cell volume (Atkinson et al., 2003). Furthermore, changes in cell morphology, such as cell shape, size, and color, due to heat stress have been observed in microalgae (Neustupa et al., 2008; Chokshi et al., 2015). The “responses to abiotic stimulus” subcategory had a higher percentage of downregulated genes than that of upregulated genes, suggesting that high temperature may not cause significant stress to B. adriatica or it adapts to high temperature rapidly.

The number of unigenes differentially expressed under only cold or heat stress from that under optimal temperature (Types 2, 4, 6, and 8) was significantly more than that detected under both cold and heat stress (Types 1, 3, 7, and 9). This suggests that the mechanism of molecular responses of B. adriatica to cold stress is different from that to heat stress. Interestingly, the number of unigenes belonging to Types 1 (upregulated under both cold and heat stress) and 9 (downregulated under both cold and heat stress), Types 2 (upregulated under cold stress and unchanged under heat stress) and 8 (downregulated under cold stress and unchanged under heat stress), Types 3 (upregulated under cold stress and downregulated under heat stress) and 7 (downregulated under cold stress and upregulated under heat stress), and Types 4 (unchanged under cold stress and upregulated under heat stress) and 6 (unchanged under cold stress and downregulated under heat stress) were similar to a mirror image. Thus, B. adriatica may balance upregulation and downregulation of gene expression in response to temperature stress.

The “apolipoprotein d lipocalin (ApoD)” identified in Type 1 is an evolutionarily conserved anti-stress protein that enhances tolerance to freezing and oxidative stress in Arabidopsis thaliana. It was upregulated in the branching stony coral Acropora millepora under cool and warm treatments and the phototrophic dinoflagellate Prorocentrum donghaiense under phosphorus-depleted conditions (Charron et al., 2008; Ganfornina et al., 2008; Muffat and Walker, 2010; Wuitchik, 2018; Zhang et al., 2018). Furthermore, the unigene encoding “alpha-protein kinase vwkA” belonged to Type 1. The alpha-kinase family is involved in diverse cellular processes such as protein translation, Mg²⁺ homeostasis, intracellular transport, cell migration, adhesion, and proliferation (Middelbeek et al., 2010). Betapudi et al. (2005) reported vwkA involvement in regulating myosin II abundance and assembly behavior in the amoeba Dictyostelium discoideum. Myosin II is an actin-dependent molecular motor protein important for cell motility and cytokinesis (Maciver, 1996; Foth et al., 2006). Therefore, temperature stress may affect B. adriatica cell movement and division processes.

The unigenes encoding “facilitated trehalose transporter Tret1” and “glutathione S-transferase class-mu 28 kDa isozyme” belonged to Type 2. Trehalose improves insect and plant tolerance to cold stress (Kikawada et al., 2007; Iordachescu and Imai, 2008; Benoit et al., 2009). Therefore, trehalose may play the same role in B. adriatica as in insects and plants. Furthermore, glutathione S-transferase (GST) is an antioxidant; thus, it is considered to be a potential oxidative stress biomarker (Guo et al., 2014). In the dinoflagellate Fugacium sp., most of the 10 GST transcripts detected were upregulated under heat stress (Gierz et al., 2017). However, in this study, the B. adriatica GST gene was highly expressed under cold stress, but not heat stress. Therefore, this gene may be expressed by temperature stress, but whether the temperature stress is cold or heat stress is dependent on dinoflagellate species.

The unigenes encoding “fatty acid desaturase” and “fructose-bisphosphate aldolase (FBA)” belonged to Type 3. Temperature change affects the membrane fluidity of poikilothermic organisms such as bacteria, fungi, protozoa, plants, and animals (Murata and Los, 1997). Membrane fluidity perceives changes in environmental conditions and transduces signals to acclimate to these conditions (Murata and Los, 1997; Los and Murata, 2004). The cell membrane becomes rigid when exposed to low temperatures and becomes fluid at high temperatures (Murata and Los, 1997). The reduced membrane fluidity at low temperatures can be improved via membrane lipid desaturation by fatty acid desaturases (Sakamoto and Murata, 2002). Conversely, the increased membrane fluidity at high temperatures can be ameliorated by integrating de novo synthesized saturated fatty acids into membrane lipids and the presence of membrane-stabilizing proteins (Los et al., 2013). Therefore, the upregulation and downregulation of this gene under cold and heat stress, respectively, in this study corresponded with previous results. FBA is a key metabolic enzyme involved in glycolysis, gluconeogenesis, and the Calvin cycle (Flechner et al., 1999; Gross et al., 1999). Furthermore, it plays an important role in stress signaling in plants (Lu et al., 2012). Therefore, B. adriatica may be affected more by cold stress than that by heat stress.

The “FAD-binding domain-containing protein” belonging to Type 4 includes several genes playing key roles in many metabolic pathways, including the citric acid cycle, fatty acid beta-oxidation, and amino acid catabolism (Mansoorabadi et al., 2007; Kim and Winge, 2013). Moreover, “malate l-lactate dehydrogenase” belonging to Type 4 also functions in many metabolic pathways, including the tricarboxylic acid cycle, glyoxylate bypass, amino acid synthesis, and gluconeogenesis (Musrati et al., 1998; Minarik et al., 2002). They were upregulated only under heat stress, thus indicating that mechanisms for acclimation to temperature stress act differentially under cold and heat stress in B. adriatica metabolic pathways.

As “proline/betaine transporter” belonging to Type 6 is an amino acid transporter, its downregulation under heat stress suggests that the capacity to export amino acids in B. adriatica cell may reduce under heat stress (Li et al., 2021).

The unigenes encoding “chloroplast light harvesting complex (LHC) protein, partial” and “intraflagellar transport (IFT) protein 46-like” belonged to Type 7. LHC proteins function in light capture and photoprotection (Boldt et al., 2012). Three upregulated LHC genes were detected in Fugacium sp. under heat stress (Gierz et al., 2016). In contrast, under cold stress, LHC activity was reduced in C. reinhardtii and cereal grain Hordeum vulgare (Atienza et al., 2004; Valledor et al., 2013). Therefore, the results of this study are consistent with those of previous studies. The IFT system assembles and maintains eukaryotic cilia and flagella (Cole, 2003). The cilia and flagella are best known to play important roles as motile organelles, but they also function as sensory perception organelles that monitor the surrounding environment, such as chemicals, light, and temperature (Silflow and Lefebvre, 2001; Follit et al., 2006). Although B. adriatica speed was not measured during incubation, when observed under a microscope, B. adriatica cells swam more slowly under cold stress than under heat stress or optimal temperature.

The unigene encoding “DnaJ-like subfamily B member 6-A” was identified to belong to Type 8. Heat shock proteins (HSPs) are important chaperones that maintain cellular homeostasis under both favorable and unfavorable conditions (Wang et al., 2004). They are considered to be produced more under unfavorable conditions, such as temperature and light stress (Rajan and D’Silva, 2009). Unlike B. adriatica, which downregulates the unigene encoding DnaJ, DnaJ expression in other dinoflagellates was upregulated. The relative DnaJ transcription level in S. acuminata increased under both low and high temperatures than that under a moderate temperature; moreover, increased DnaJ expression in heat stress-exposed Cladocopium sp. and Fugacium kawagutii have been reported (Levin et al., 2016; Deng et al., 2019; Lin et al., 2019). In contrast, Gierz et al. (2017) reported that several transcripts encoding HSPs, including DnaJ, were identified in Fugacium sp. under heat stress, and the downregulated transcripts were more dominant than upregulated transcripts. Therefore, whether the expression of this gene is upregulated or downregulated is dependent on the dinoflagellate species.

The unigene encoding “troponin C, skeletal muscle” belonged to Type 9. Troponin C is a Ca²⁺ homeostasis protein of the calmodulin family (DeSalvo et al., 2010; Abassi et al., 2021). In the elkhorn coral Acropora palmata, this gene was downregulated under heat stress (DeSalvo et al., 2010). The downregulation of this gene under cold and heat stress in this study showed that Ca²⁺ homeostasis in B. adriatica may be affected by cold and heat stress. To elucidate this, more studies on the effect of temperature on the expression of this gene in dinoflagellates are required.

To better understand the response of dinoflagellate communities to diverse environmental changes, including cold waves, heat waves, ocean acidification, and local pollution, both eco-physiological and molecular characteristic changes should be explored together.