Temperature Differentially Affects Gene Expression in Antarctic Thraustochytrid Oblongichytrium sp. RT2316-13

Oblongichytrium RT2316-13 synthesizes lipids rich in eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA). The content of these fatty acids in the total lipids depended on growth temperature. Sequencing technology was used in this work to examine the thraustochytrid’s response to a decrease in growth temperature from 15 °C to 5 °C. Around 4% (2944) of the genes were differentially expressed (DE) and only a few of the DE genes (533 upregulated; 206 downregulated) had significant matches to those in the SwissProt database. Most of the annotated DE genes were related to cell membrane composition (fatty acids, sterols, phosphatidylinositol), the membrane enzymes linked to cell energetics, and membrane structure (cytoskeletal proteins and enzymes). In RT2316-13, the synthesis of long-chain polyunsaturated fatty acids occurred through ω3- and ω6-pathways. Enzymes of the alternative pathways (Δ8-desaturase and Δ9-elongase) were also expressed. The upregulation of the genes coding for a Δ5-desaturase and a Δ5-elongase involved in the synthesis of EPA and DHA, explained the enrichment of total lipid with these two long-chain fatty acids at the low temperature. This molecular response has the potential to be used for producing microbial lipids with a fatty acids profile similar to that of fish oils.


Sequence Analysis and Transcriptome Annotation
The total output of HiSeq 4000 for all the samples was approximately 174 million paired-end reads, with an average of 29 million reads for each sample. The quality filtering step retained over 99.5% of the reads (173.7 million), confirming an excellent sequencing performance and assuring a minimum of 3 Gbps of data (20 million reads) per sample (Table S1, Supplemental Material). The final transcriptome resulted in 118,595 transcript sequences that were contained in 69,220 genes. The quality metrics indicated that 50% of the assembled bases were contained in transcripts longer than 2518 bp (N50), with an average transcript length of 1243 bp, and a total size of 147 Mbp. The presence of several short transcripts in the assembly was evidenced by a median length metric of 554 bp (half of the assembled transcripts had a length of 554 bp or less) and this may have affected the average length of the transcriptome ( Table 1). The results from BUSCO (Benchmarking Universal Single-Copy Orthologs [17]) indicated a high degree of transcriptome completeness, resulting in 235 of the 255 genes used for evaluation being found ( Table 2).  [17]).
Alignment of high quality reads to the assembled transcriptome of RT2316-13 showed a high mapping rate: an average of 91.63% for the properly aligned paired reads, and an overall mapping rate > 99% for every sample. The overall alignment included reads that contained one mismatch in their sequences, and reads that were aligned separately (Table S2, Supplemental Material).

Differentially Expressed Genes
Using a false discovery rate (FDR) of 0.05 and a fold change (FC) of 4 as thresholds, 2944 genes were found to be differentially expressed (DE). Of these DE genes, 2062 were induced by the low temperature (LT) whereas 882 genes were repressed (Excel file Paredes_et_al.xlsx, Supplemental Material). The homology searches in the public database SwissProt indicated that relatively few of the DE genes (25.1%) had significant matches, probably because of a comparative rarity of microbes of the genus Oblongichytrium. Of the annotated genes (1999 transcripts), 533 were induced and 206 were repressed (Excel file Paredes_et_al.xlsx, Supplemental Material). Most of the DE genes belonged to the GO-slim (cut-down versions of the Gene Ontologies) cellular component category ( Figure 1a). The DE proteins were classified into 16 classes (Figure 1b). Among the induced proteins, the class "Nucleic acid binding" (PC00031, p = 0.0187) was depleted, while "Cytoskeletal protein" (PC00085, p = 0.0006) was overrepresented. In the repressed protein group the overrepresented classes were the "Translation protein" (PC00263, p < 0.0001) and its subclass "Ribosomal protein" (PC00202, p < 0.0001). The subclass "Oxidoreduction", whose parent was "Metabolite interconversion enzyme", was also overrepresented (PC00176, p = 0.0022). Genes relating to carbohydrate, amino acids, energy, and lipid metabolism represented 64.4% of the annotated genes ( Figure 1c). In the lipid metabolism group, the genes relating to fatty acid degradation, and metabolism of glycerophospholipid and glycerolipid, accounted for 44.5% of the annotated genes whereas 6.2% of the annotated genes related to fatty acid elongation (Figure 1d). Annotated genes related to carbohydrate metabolism (21%, Figure 1e) were the most abundant and within this group more than 52% corresponded to genes coding for enzymes related to glycolysis, inositol phosphate, citrate cycle, and pyruvate metabolisms.

Cytoskeletal Proteins
Forty-eight genes coding for the cytoskeleton-related proteins were DE. The cytoskeleton of eukaryotes is a dynamic network of interlinking protein filaments that confer shape and mechanical resistance to cells, and allows them to move. Cytoskeleton is also involved in many cell signaling pathways and intracellular transport. The cytoskeleton is composed of microfilaments (linear polymers of G-actin proteins), intermediate filaments, and microtubules (MT, polymers of tubulin). In RT2316-13, the genes coding for tubulin α-chain (2 genes) (FC = 10.88; 16.46), tubulin tyrosine ligase (FC = 11.73), tubulin polymerization-promoting protein family member (FC = 17.02), and actin (2 genes) (FC = 8.47; 27.03) were upregulated.
Cilia are evolutionarily conserved, MT-based organelles that play important roles in cell movement, environment sensing, and signal transduction. Motile biflagellate zoospores are characteristic of most thraustochytrids [18], including RT2316-13 (microscopic observation). Eighteen genes coding for cilia and flagella-associated proteins were upregulated (FC in the range of 5.34 to 18.25). The assembly, maintenance, and functions of cilia depend on the intraflagellar transport (IFT),

Cytoskeletal Proteins
Forty-eight genes coding for the cytoskeleton-related proteins were DE. The cytoskeleton of eukaryotes is a dynamic network of interlinking protein filaments that confer shape and mechanical resistance to cells, and allows them to move. Cytoskeleton is also involved in many cell signaling pathways and intracellular transport. The cytoskeleton is composed of microfilaments (linear polymers of G-actin proteins), intermediate filaments, and microtubules (MT, polymers of tubulin). In RT2316-13, the genes coding for tubulin α-chain Cilia are evolutionarily conserved, MT-based organelles that play important roles in cell movement, environment sensing, and signal transduction. Motile biflagellate zoospores are characteristic of most thraustochytrids [18], including RT2316-13 (microscopic observation). Eighteen genes coding for cilia and flagella-associated proteins were upregulated (FC in the range of 5.34 to 18.25). The assembly, maintenance, and functions of cilia depend on the intraflagellar transport (IFT), the bidirectional Mar. Drugs 2020, 18, 563 5 of 16 trafficking of vesicles and proteins along the ciliary axoneme [19]. IFT is mediated by motor protein complexes, the cytoplasmic dyneins, and kinesins. In RT2313-13, the expression of genes coding for kinesin-related proteins (19 genes) (FC in the range of 4.41 to 15.84), and dynein-related proteins (24 genes) (FC in the range of 4.49 to 19.96) were upregulated.
The contribution of cytoskeleton to the cell's structural integrity, organization, and movement explained the significant number of the DE genes in RT2316-13 in response to a change in growth temperature. Most of the annotated genes were upregulated which suggested that the cells either needed to adapt or reshape their structure, or elevate the level of certain intracellular enzymes, to enable movement at the low temperature.

Inositol Phosphate Metabolism
Phosphatidylinositol (PI) is a family of lipids that are plentiful in cell membranes [20]. In one of its roles, PI modulates protein function, depending on the phosphorylation state of inositol. Deletion of INP51 (a gene coding a protein with 5-phosphatase activity against phosphatidylinositol 4,5-bisphosphate, PI(4,5)P2) in Saccharomyces cerevisiae resulted in accumulation of PI(4,5)P2 that correlated to cold-tolerance through an unknown mechanism [21]. In RT2316-13, the genes encoding phosphatidylinositol 4-phosphate 5-kinase (FC = 4.79; 7.92; 18.75) and phosphatidylinositol 4-kinase (FC = 8. 19), and the enzymes involved in PI(4,5)P2 synthesis, were upregulated, suggesting a role in cold-tolerance as previously observed in S. cerevisiae. PI(4,5)P2 is hydrolyzed by phospholipase C (PLC), an enzyme known to activate at low temperatures [22,23], to produce the secondary messengers diacylglycerol (DAG) and inositol-1,4,5-trisphosphate. In RT2316-13 a gene coding for a PLC (FC = 13.77) was upregulated. Once formed, DAG activated protein kinase C (PKC), an enzyme that has a central role in the cell wall integrity signaling pathway that responds to cell surface stress in S. cerevisiae [24]. In RT2316-13 a gene coding for a Ca 2+ -independent, phospholipid-and DAG-dependent serine/threonine-protein kinase (novel PKC) was upregulated (FC = 16.67). Upregulation of the genes relating to the synthesis of PI(4,5)P2, PLC and PKC suggested LT as a stress factor for RT2316-13.
Four genes coding for enzymes in the pentose phosphate pathway were DE in RT2316-13. Among these, a gene coding for NADPH-producing 6-phosphogluconate dehydrogenase, a source of NADPH for fatty acid biosynthesis, was downregulated (FC = −17.01). On the other hand, genes coding for a ribose-phosphate pyrophosphokinase 4 (FC = 4.30), involved in nucleotide biosynthesis, and ribokinase (FC = 11.30) which catalyzes the phosphorylation of ribose to ribose-5-phosphate, were upregulated. The upregulation of the ribokinase gene could be related to a possible degradation of ribosomes to release nutrients [26]. The LT treatment downregulated 33 genes coding for ribosomal proteins in RT2316-13. This was because a large number of ribosomes were needed only during rapid growth [27] and the slowed growth at low temperature promoted catabolic recycling of the older ribosomes to release amino acids and nucleotides. Nucleotides could be dephosphorylated, and hydrolyzed to ribose and nucleobases [28]. The action of ribokinase produces ribose-5-phosphate that can be converted to glycolytic intermediates (fructose-6-phosphate and glyceraldehyde-3-phosphate) in the nonoxidative pentose phosphate pathway.
The LT treatment downregulated four genes coding for enzymes in the citric acid cycle (fumarate hydratase, FC = −19.48; malate dehydrogenase, FC = −16.03; succinateCoA ligase, FC = −14.91; and aconitate hydratase, FC = −14.67), whereas a gene coding for isocitrate dehydrogenase [NADP+] (FC = 6.38) was upregulated. These results showed that the LT treatment depressed the activity of the TCA cycle, reducing energy supply and the production of citric acid. The decreased activity of the TCA cycle and NADPH-producing 6-phosphogluconate dehydrogenase explained the reduced lipid content (from 33.5% to 15.2% w/w) in RT2316-13 biomass grown at 5 • C [14].
Other similar observations have been reported. An earlier study with Aurantiochytrium found a reduced supply of biochemical energy from glycolysis and TCA cycle, as a consequence of cold stress (15 • C compared with 25 • C) [29]. The authors further observed that in the logarithmic growth phase, the expression of genes coding for a ribose 5-phosphate isomerase and ribose-phosphate pyrophosphokinase, the enzymes required for nucleotide synthesis, was upregulated.

Respiratory Chain and Oxidative Phosphorylation
Several genes of the respiratory chain were downregulated by the LT treatment. These were the genes coding for subunits of NADH dehydrogenase (ubiquinone) (
A homology search for genes in the PKS system involved in the EPA and/or DHA biosynthesis (pfaA, pfaB, pfaC, and pfaD) in the transcriptome of RT2316-13 was performed. The one significant hit after being extracted and queried by homology on GenBank and UniProt databases matched to a polyketide synthase (PKS) gene related to antibiotic biosynthesis. Thus, under the growth conditions, LCPUFA synthesis in RT2316-13 was not due to the PKS system, unlike in the related thraustochytrids Schizochytrium sp. [33] and Thraustochytrium sp. SZU445 [34].
Twenty-five genes coding for enzymes related to glycerolipid metabolism and biosynthesis of fatty acids and steroids were DE in RT2316-13. Among the DE genes coding for enzymes involved in the synthesis of PUFA were those coding for: a very-long-chain 3-oxoacyl-CoA reductase (FC = 8.32), a long-chain fatty acid CoA ligase (2 genes) (FC = 4.35; 6.08), ELOVL6 (FC = 6.41), ELOV2 (FC = 5.29), and a ∆8-fatty-acid desaturase (FC = 5.10). The latter gene was found to be highly similar to a ∆5-desaturase, a key enzyme in the synthesis of EPA and DHA. Upregulation of this gene together with upregulation of the gene coding for ELOV2 (a ∆5-elongase) explained the enrichment of the total lipid in EPA and DHA at LT (Figure 2). Two genes coding for diacylglycerol O-acyltransferase 2, the enzyme required for the synthesis of triacylglycerols (TAG), were DE (FC = 14.62; −6.88), whereas two genes coding for monoacylglycerol lipase (MAGL) were upregulated (FC = 5.63; 12.24). MAGL catalyzes the last step in the hydrolysis of storage TAG and has been postulated to have a role in the detoxification of deleterious membrane lipids [38] which could form through oxidation of membrane PUFA because of the increased solubility of oxygen in water at LT [39,40].
In RT2316-13, a gene coding for cardiolipin synthase A was upregulated (FC = 4.60). Cardiolipins are phospholipids of the inner mitochondrial membrane required for the optimal activity of oxidative phosphorylation complexes and the ATP/ADP carrier [41]. Interactions of sterols with other membrane components contribute not only to regulate membrane fluidity and permeability, but also to the activity of membrane-bound enzymes and transporters [42]. Disruption of the SYRl/ERG3 gene that encoded sterol C-5 desaturase, an enzyme involved in the ergosterol biosynthesis, inhibited the growth of S. cereviciae at LT possibly because sterol was required for proper functioning of the tryptophan transporter at LT [43]. A similar function could explain the upregulation (FC = 20.09) of a ∆7-sterol-C5(6)-desaturase, involved in the biosynthesis of sitosterol and campesterol in RT2316-13.

Validation of Gene Expression Profiles by qRT-PCR
The amplification of four selected genes that coded for acetyl-CoA carboxylase (AAC-5, AAC-15), malic enzyme (Me-5, Me-15), subunit A of PKS (Pfaa-5, Pfaa-15), and subunit C of LCPUFA PKS (Pfac-5, Pfac-15), in cDNA of the biomass grown at 5 • C and 15 • C was evaluated. The genes for AAC and Me were selected because the first of these is involved in the synthesis of fatty acid precursors, and the second one is involved in the synthesis of NADPH required for the synthesis of PUFA. Genes coding for acetyl-CoA carboxylase and malic enzyme were expressed at both the low and control temperatures ( Figure S1). The lack of expression of the genes coding for the two subunits of LCPUFA PKS ( Figure S1, Supplemental Material) was consistent with the results of transcriptome analysis. Five other genes coding for acetyl-CoA carboxylase and four other genes coding for malic enzyme were annotated in the transcriptome of RT2316-13, but these were not differentially expressed.

Microorganism and Culture Conditions for RNA Extraction
The thraustochytrid Oblongichytrium sp. RT2316-13 used in this work had been isolated from seawater samples collected near the Antarctic Base Professor Julio Escudero (62 • 12 57 S, 58 • 57 35 E). Details of isolation and identification were previously published [15]. Briefly, the samples were supplemented with sterile pine pollen (~15 mg per 200 mL) and kept at 5 • C. After 15 days, pollen grains were collected by filtration (0.45 µm pore size) and dispersed on a solid medium supplemented with streptomycin sulfate and penicillin G (0.3 g/L of each) (Sigma, St. Louis, MO, USA) to prevent proliferation of bacteria. Pure isolates were obtained through repeated subculture of individual colonies on solid medium. Cell DNA was isolated and the sequence of the 18S rRNA gene was amplified using the method described previously [15]. Results of the Sanger sequencing were assembled using the Geneious 4.8.4 program (Biomatters Ltd., Auckland, New Zealand) and compared with sequences in EMBL/DDBJ/PDB/GenBank databases by BLASTN 2.2.21. Pure stock cultures were kept frozen at -80 • C in 50% glycerol (v/v).
Cells were grown aseptically. The culture medium comprised of glucose (20 gL −1 ; Merck, Darmstadt, Germany), yeast extract (10 gL −1 ; BBL™, Becton, Dickinson and Co., Sparks, MD, USA), and monosodium glutamate (0.6 gL −1 ; Merck) in a 1:1 by volume mixture of artificial sea water (ASW) [44] and distilled water. Trace elements and vitamins were added to the medium as previously specified [45]. Prior work had shown this medium to support growth with a highest attainable total lipid content in the biomass (35% of dry weight) at 15 • C [14].
Inoculum was prepared by transferring 1 mL of the thawed stock to 100 mL of the earlier specified liquid medium in a 250 mL Erlenmeyer flask. The flask was incubated (15 • C) on an orbital shaker at 150 rpm for 7 days. This culture (5 mL) was used to inoculate 100 mL of the same fresh medium and incubated (15 • C, 150 rpm) for 2 days. Afterwards, 50 mL aliquots of the grown culture were transferred to two flasks, each containing 50 mL of the same fresh medium. One of the flasks was incubated at 15 • C (control culture) and the other was incubated at 5 • C (low temperature (LT) culture) in an orbital shaker at 150 rpm. Both flasks were harvested after 3 days.
Biomass was harvested by centrifugation (7000× g, 4 • C, 10 min). The recovered cell pellet was suspended in 12 mL of a 19:1 (by vol) mixture of aqueous phenol (0.8 M) and absolute ethanol, and stored at −20 • C until the extraction of the nucleic acids. Three replicate experiments were used to produce the control and LT biomass samples.

RNA Extraction
Cells from the previous step were lysed in a high-pressure homogenizer (Stansted Fluid Power Ltd., London, England) (one pass, 20 bar) under sterile conditions. Cell debris was removed from the lysate by centrifugation (10,000× g, 4 • C, 15 min). Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific Inc., Waltham, MA, USA) according to the manufacturer's instructions. The collected RNA was purified using RNA purification columns (Thermo Fisher Scientific Inc., Waltham, MA, USA) following the manufacturer's instructions.
Sample degradation and integrity were monitored by agarose gel electrophoresis (1% w/v, 100 mV, 30 min). The purity and recovery yield of RNA were calculated using the absorbance values measured at 260 and 280 nm using a UV spectrophotometer (BioTek Synergy HT Microplate Reader; BioTek Instruments, Inc., Winooski, VT, USA).

Conversion of RNA to cDNA
The extracted RNA was treated with DNase to eliminate any contamination. One µL of random primers, mainly random hexanucleotides (500 µg mL −1 ; Promega Corporation, Madison, WI, USA) were added. A conversion mix (M-MLV RT 5×, M-MLV RT transcriptase, RNasin and dNTPs) was added. The first PCR program (MultiGene optiMAX; Labnet International Inc., Edison, NJ, USA) was: 5 min at 20 • C; 45 min at 37 • C; 5 min at 95 • C; and 4 • C until storage. A PCR using the primers for 18S rRNA subunit (Table S3, Supplemental Material) was performed to verify the correct conversion. Conditions of the second PCR program were: 3 min at 95 • C; 30 cycles of 1 min each at 94 • C followed by 1 min at 61 • C, and 1 min at 72 • C; 10 min at 72 • C; and 4 • C until storage. The cDNA was stored for validation assays.

Sequencing and Quality Control of RT2316-13 Transcriptome
A total of six samples were sequenced (Genoma Mayor, Santiago, Chile) in an Illumina HiSeq 4000 (Illumina, Inc., San Diego, CA, USA) platform for 150 cycles in paired-end mode. The sequencing outputs, in FastQ file format, were processed using TrimGalore software to remove adaptors and reads that were likely to be sequencing errors. The reads, which on average had a Phred score of less than 30 (1 error in 1000 bp), were removed from downstream analysis.
The remaining high quality reads, selected during the quality control step, were used to construct a de novo assembly of RT2316-13 transcriptome using the Trinity v2.9.1 software [46]. The software used de Bruijn graphs to reconstruct transcripts and their variants using an extensive K-mer search strategy. The assembly integrity was assessed using a set of conserved orthologous gene sequences contained in the Eukaryota OrthoDB v10 database (www.orthodb.org), which were queried by similarity to the assembled transcripts, as in the BUSCO approach (Benchmarking Universal Single-Copy Orthologs) described by Simão et al. [17].

Differential Gene Expression Analysis
The relative abundance was calculated using RSEM v1.2.26 [47]. The resulting abundance for each sample was merged in a matrix, and analyzed with the Bioconductor package edgeR in the R statistical environment [48]. The significance of changes in gene expression was judged using a False Discovery Rate value (FDR) of less than 0.05 and a minimum fold change (FC) of 4 as thresholds.

Transcriptome Annotation and Enrichment Analysis
All resulting genes were aligned into the UniProt/SwissProtKB database using BLAST+ with an e-value of 1e −10 as the threshold. Functional annotation and ontology assignments were performed using the PANTHER classification system [49] with gene lists obtained from blast results (top hit) aligned with the reference proteome database (version 2018_4) from the European Molecular Biology Laboratory (EMBL). To search for enrichment or depletion of an ontology term among differentially expressed genes, they were analyzed using GOSeq [50].

Validation of the Transcriptome Results
Conventional PCR was performed to confirm or discard the presence of genes associated with LCPUFA synthesis. Five pairs of primers were used: Subunits A and C of PKS, malic enzyme, acetyl-CoA carboxylase, and, as a housekeeping gene, the subunit 18S rRNA. The sequences were obtained from Ma et al. [29] (Table S1, Supplemental Material). Green Master Mix (Promega, Madison, WI, USA) was used for PCR in accordance with the manufacturer's instructions. The reaction was carried out in a final volume of 20 µL using 4 µL of cDNA synthesized as described above. Conditions of PCR program 2 (described above) were used. Results were verified by agarose gel (1% w/v) electrophoresis.

Data Availability
The raw reads of the Sequence Read Archive (codes from SRR12600987 to SRR12600992) have been deposited at NCBI BioProject (accession number PRJNA642733).

Conclusions
This study contributes to knowledge of the cold water marine thraustochytrid Oblongichytrium RT2316-13. This microorganism grows abundantly at 5 • C, a feature that is at least partly attributed to its ability to synthesize long chain omega-3 PUFA. Transcriptome analysis of Oblongichytrium RT2316-13 revealed a significant number of differentially expressed (DE) genes (2944) in response to a decrease in temperature. Most of the DE genes related to the composition of the cell membrane (fatty acids, sterols, phosphatidylinositol), the membrane-associated enzymes, and the membrane-linked cytoskeletal proteins and enzymes. This suggests that in this thraustochytrid, the cell membrane has a role in "sensing" temperature and is an active part of the temperature signal transduction system.
Thraustochytrids are mainly studied because of their capability to synthesize docosahexaenoic acid (DHA). However, strains such as RT2316-13 produce lipids rich in both EPA and DHA, fatty acids with known health effects. An understanding of temperature modulation of the genes related to the biosynthesis of these fatty acids is important for designing efficient production processes. The present work shows that in RT2316-13 the synthesis of EPA and DHA occurs through the action of elongases and desaturases and thus differs from synthesis in Aurantiochytrium species. The upregulation of genes coding for a ∆5-desaturase and ∆5-elongase involved in the synthesis of EPA and DHA explained the enrichment of total lipid with these two long chain polyunsaturated fatty acids at the low temperature.
Supplementary Materials: The following are available online at http://www.mdpi.com/1660-3397/18/11/563/s1, Figure S1: Amplification of cDNA segments. AAC, acetyl-CoA carboxylase; Me, malic enzyme; Pfaa, subunit A of PKS gene; Pfac, subunit C of PKS gene; 18s, 18S subunit of RNA. xx-5 and xx-15 denote cDNA of cells grown at 5 • C and 15 • C, respectively, Table S1: Sequencing results. The initial amount of reads (total reads) and the number of reads that passed quality control filters (HQ reads); Table S2: Total number of high quality (HQ) reads that were mapped to a sequence in RT2316-13 transcriptome, concordantly one or more times (properly aligned), and the total number of reads aligned (overall aligned), with their respective mapping rate (%); Table S3: Primers used for identifying genes related to production of fatty acids, Excel file: Paredes_et_al.xlsx.
Author Contributions: C.S. conceived the research project leading to this publication. G.L. performed bioinformatics. P.P., A.L., L.F., and C.G.I. performed the experiments and prepared the initial draft of the manuscript. C.S., J.A.A., and Y.C. discussed the results and wrote the final manuscript. All authors have read and agreed to the published version of the manuscript.