Genome-Wide Transcriptional Changes of Rhodosporidium kratochvilovae at Low Temperature

Rhodosporidium kratochvilovae strain YM25235 is a cold-adapted oleaginous yeast strain that can grow at 15°C. It is capable of producing polyunsaturated fatty acids. Here, we used the Nanopore Platform to first assemble the R. kratochvilovae strain YM25235 genome into a 23.71 Mb size containing 46 scaffolds and 8,472 predicted genes. To explore the molecular mechanism behind the low temperature response of R. kratochvilovae strain YM25235, we analyzed the RNA transcriptomic data from low temperature (15°C) and normal temperature (30°C) groups using the next-generation deep sequencing technology (RNA-seq). We identified 1,300 differentially expressed genes (DEGs) by comparing the cultures grown at low temperature (15°C) and normal temperature (30°C) transcriptome libraries, including 553 significantly upregulated and 747 significantly downregulated DEGs. Gene ontology and pathway enrichment analysis revealed that DEGs were primarily related to metabolic processes, cellular processes, cellular organelles, and catalytic activity, whereas the overrepresented pathways included the MAPK signaling pathway, metabolic pathways, and amino sugar and nucleotide sugar metabolism. We validated the RNA-seq results by detecting the expression of 15 DEGs using qPCR. This study provides valuable information on the low temperature response of R. kratochvilovae strain YM25235 for further research and broadens our understanding for the response of R. kratochvilovae strain YM25235 to low temperature.


INTRODUCTION
One of the crucial ecological factors associated with abiotic stress is low temperatures, which can disrupt microbial homeostasis and affect the biological functions of cells (García-Ríos et al., 2017;Wang et al., 2020). In addition, low temperatures can severely inhibit fungal growth and may also kill it (Wang et al., 2014;Long et al., 2015). Although 10-18 • C is considered as a low temperature range for fungus, these are permissive temperature (García-Ríos et al., 2017). Currently, the understanding of the molecular mechanisms of response in fungus to low temperature stress has been explored. For example, several fungal genera adapt to low temperature by changing the expression of different genes as part of the process known as the cold-shock response that has been investigated in Saccharomyces cerevisiae (Aguilera et al., 2007). During adaptation, several cold response-related genes, including those involved in energy preservation, detoxification, osmolyte production, protein folding support, and maintenance of membrane fluidity, are activated (Murata et al., 2006). Growth capacity diminishes as the temperature decreases below the optimum temperatures. Studies in Neurospora crassa have demonstrated a high accumulation of mRNAs of carotenoid genes to reduce oxidative stress under low temperature (Castrillo et al., 2018). These results suggest that differentially expressed genes (DEGs) play a key role under low temperature stress in fungi. Therefore, the identification of genes associated with low temperature response and intensive exploration for molecular mechanisms in response to low temperature at the level of gene expression in fungi are necessary.
Rhodosporidium species have been reported to be capable of synthesizing some value-added compounds with a wide industrial usage, such as lipid (Uprety et al., 2017), carotenoids (Buzzini et al., 2007), enzymes (MacDonald et al., 2016), polyunsaturated fatty acids (PUFAs) (Ageitos et al., 2011;Cui et al., 2016;Wang et al., 2017;He et al., 2019) and sugar alcohols including D-arabitol and galactitol (Jagtap and Rao, 2018;Jagtap et al., 2019). As a species of the genus Rhodosporidium, Rhodosporidium kratochvilovae strain YM25235 (isolated from Chenghai Lake, Yunnan, China) is a cold-adapted oleaginous fungal yeast strain that can grow at 15 • C (He et al., 2019), and it is capable of producing PUFAs (Cui et al., 2016;Wang et al., 2017;He et al., 2019). PUFAs are involved in the maintenance of optimal physical and biological properties of cell membranes, which is an essential adaptation strategy in microorganisms against cold stress (Clarke, 2001;Sampath and Ntambi, 2004). Certain fungi produce high levels of PUFAs in response to low temperature (Calvo et al., 2001;Wang et al., 2017). We had previously reported that the strain YM25235 at 15 • C produced high amounts of linoleic acid (LA, C18:2 9,12) and α-linolenic acid (ALA, C18:3 9,12,15), and the inhibition of PUFA biosynthesis negatively influenced the cold adaptation of YM25235 (Cui et al., 2016;Wang et al., 2017). In addition, cold-adapted fungal yeasts have attracted the wide attention of several scientists worldwide because of their significant potential for application in diverse industries (Margesin and Miteva, 2011;Carrasco et al., 2017). These fungal species have evolved physiological strategies to survive in cold climates, including modulation of enzyme kinetics and membrane fluidity (Firdaus-Raih et al., 2018).
For R. kratochvilovae, Miccoli et al. (2018) have completed whole genome sequencing of the strain LS11, and several other Rhodosporidium strains have also been sequenced (Hu and Ji, 2016;Zhang et al., 2016;Tran et al., 2019). In addition, multiomics analyses of the oleaginous fungal yeast Rhodosporidium toruloides were conducted to shield lights on Pi-limitationinduced lipid accumulation (Wang et al., 2018). By comparing microbial growth and gene expression patterns, several studies demonstrated that expression levels of genes that encode carotenoid biosynthesis were altered in response to light change (Pham et al., 2020). Dinh et al. (2019) reconstructed a genomescale model of R. toruloides IFO0880's metabolic network. However, changes in the gene expression profile of these fungal yeasts, especially R. kratochvilovae, under temperature stress have not been investigated.
We performed this study to gain a comprehensive understanding of the molecular adaptation mechanisms to low temperatures in R. kratochvilovae strain YM25235 from the perspective of gene expression. The whole genome of R. kratochvilovae strain YM25235 was sequenced. Subsequently, based on RNA-seq data and the above-sequenced genome of strain YM25235, DEGs were identified between R. kratochvilovae at low temperature (15 • C) and normal temperature (30 • C) groups. Furthermore, we performed annotation of R. kratochvilovae strain YM25235 genes and functional enrichment analysis [Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway] of DEGs. The DEGs were summarized, and several key pathways involved in low temperature stress tolerance were obtained from the list of KEGG terms enriched by the DEGs. This study will provide insights into the molecular adaptation mechanisms of R. kratochvilovae strain YM25235 under low temperature stress.

Culture Conditions and Cell Sample Preparation
For whole genome sequencing, R. kratochvilovae strain YM25235 was grown in yeast extract peptone dextrose (YPD: 1% yeast extract, 2% peptone, and 2% glucose) broth at 30 • C to logarithmic growth phase (OD 600 = 1.60). For RNA-seq, R. kratochvilovae strain YM25235 was pre-incubated in YPD broth at 30 • C for 24 h and then was exposed to 15 and 30 • C for another 8 h, with the final OD 600 values of ∼2.1 and ∼2.3, respectively, with three technical replicates collected for each sample. The culture broth was centrifuged at 5,000 × g for 10 min at 4 • C. For cell disruption, the culture pellets were ground in liquid nitrogen.

DNA Extraction, Genome Sequencing, and Assembly
Total DNA of R. kratochvilovae strain YM25235 was extracted using Wizard Genomic DNA Purification Kit (Promega, Madison, WI, United States). The next experimental procedures were performed according to the standard protocols provided by Oxford Nanopore Technologies (ONT). To optimize the sequencing experiments and improve the throughput, a library was constructed using the SQK-LSK109 Ligation Sequencing Kit. Then, the library was sequenced on the Nanopore PromethION platform. After sequencing, the downstream sequencing data were analyzed by basecalling programs using the Albacore software from the MinKNOW package to convert the raw sequencing data from FAST5 format to FASTQ format (Wick et al., 2019). Further filtering for the adaptor, low quality, and short reads (<2,000 bp in length) resulted in total dataset clean reads. Canu v1.5 software was used to correct the filtered subreads (Koren et al., 2017). Next, we assembled the subreads after error correction using the wtdbg software to obtain the final genome with high accuracy (Ruan and Li, 2019). The BUSCO v2.0 software was used to assess the completeness of the R. kratochvilovae strain YM25235 genome assembly (Simão et al., 2015).

Functional Annotation of the Genome
Because of the relatively low conservation of repeat sequences among different species, predictions of repeat sequences for a particular species require the construction of specific repeat sequence databases. Therefore, we used four software, including LTR_FINDER v1.05, MITE-Hunter, RepeatScout v1.05, and PILER-DF v2.4, to construct a repeat sequence database for the genome of R. kratochvilovae strain YM25235 based on the methods of structural and de novo predictions (Edgar and Myers, 2005;Price et al., 2005;Xu and Wang, 2007;Han and Wessler, 2010). Classification of the repeat sequences was performed using the PASTE Classifier software (Wicker et al., 2007). The type of the classified repeat sequences was subsequently annotated using the Repbase database (Jurka et al., 2005). Repeat sequences in the genome of R. kratochvilovae strain YM25235 were searched in the RepeatMasker v4.0.6 database (Chen, 2004).
To obtain functional annotation information of R. kratochvilovae strain YM25235 genes, the predicted gene sequences were annotated by searching them in the following gene collections using the BLAST tool: Cluster of orthologous groups for eukaryotic complete genomes (KOG) (Tatusov et al., 2000), KEGG (Kanehisa et al., 2004), Swiss-Prot (Boeckmann et al., 2003), and non-redundant protein sequence database (Nr) (Deng et al., 2006). Based on the searched results from the Nr database, the GO annotation information of genes was extracted using the in-house Perl script (Ashburner et al., 2000;Conesa et al., 2005). In addition, domains of the protein sequence of each gene were searched in the Pfam database using the Hmmer software (Eddy, 1998;Finn et al., 2016).

Extraction of RNA, Library Construction, and Sequencing
Total RNA of R. kratochvilovae strain YM25235 was purified using RNeasy Mini Kit (Qiagen, Valencia, CA, United States). After extraction of total RNA from the samples, the quality of total RNA [i.e., RNA integrity number (RIN) values, 28S/18S and concentration] was determined using the Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit). The purity of the samples was assessed using a UV spectrophotometer NanoDrop TM , ND-1000 (Thermo Scientific, United States). Next, the cDNA library was constructed. In brief, R. kratochvilovae strain YM25235 mRNA was enriched using magnetic beads with oligo (dT). The resulting mRNA was fragmented by adding an appropriate amount of the shearing reagent under high temperature conditions so that the interrupted mRNA was used as a template for the first-strand cDNA synthesis. Subsequently, the second strand cDNA was synthesized, purified using commercial kits, and sticky ends were repaired. The base "a" was added to the 3 end of the cDNA and adaptors were ligated, followed by fragment size selection and PCR amplification. The quality of the constructed libraries was checked using the Agilent 2100 Bioanalyzer (Agilent Technologies, United States) and ABI StepOnePlus Real-Time PCR System. Raw reads generated by RNA sequencing (RNA-seq) were filtered to obtain clean data by removing low-quality reads, adapter contamination, and reads with excessive unknown bases (Cock et al., 2010). Clean reads were mapped to the  (Kim et al., 2015).

Quantification of Gene Expression and Identification of Differentially Expressed Genes
Clean reads were mapped to the reference gene sequences using Bowtie 2 (Langmead and Salzberg, 2012). Next, the expression of genes was assessed based on Fragments Per Kilobase per Million (FPKM) values calculated by RSEM (RNA-Seq by Expectation Maximization) methods (Li and Dewey, 2011). DEGs were identified using the DESeq2 software, an R program suitable for the identification of DEGs from high-throughput sequencing data between the control and treatment groups (Love et al., 2014). The DEGs with fold change (FC) ≥ 2 (| log 2 ratio| ≥ 1) and p-values (Wald test in Deseq2) corrected by a false discovery rate (corrected p-values) < 0.05 were identified as DEGs.

Functional Enrichment Analysis of Differentially Expressed Genes
To further elucidate the biased biological function of DEGs and signaling pathways, the functional analyses of DEGs were performed using the Blast2GO pipeline (Conesa et al., 2005). The annotated functional terms of DEGs were clustered through the Functional Annotation Clustering tool in the KOBAS software (Xie et al., 2011). The KEGG pathways from each clustered set with an FDR-value < 0.01 were considered to be significantly enriched. Furthermore, the redundant GO terms were systematically discarded using GO trimming (Jantzen et al., 2011).

Validation of RNA-Seq Data and Gene Expression Analysis by qPCR
To validate the results of RNA-seq, 15 DEGs were randomly selected for downstream quantitative real-time PCR (qPCR) analysis, with three biological replicates (Table 1). Parallelly, we obtained the corresponding control sample in the same manner.
Purification and evaluation of total RNA were performed using the above-mentioned methods ("Extraction of RNA, Library Construction, and Sequencing"). This was performed thrice independently to generate three biological replicates. Specific primers used in the qPCR analysis were designed automatically using the Beacon Designer 7 (Supplementary Table 1). For the preparation of cDNA libraries, the first-strand cDNA was synthesized from total RNA using the PrimeScript RT reagent kit with gDNA eraser (TaKaRa, Tokyo, Japan). Next, we used RNase-free water to dilute the preliminary cDNA solution (final concentration, 100 ng/µL). Next, qPCR was performed on the Light Cycler 480 Real-Time PCR System (Roche, Basel, Switzerland) using the Bestar SYBR Green qPCR Master mix (DBI Bioscience, Shanghai, China). The qPCR reaction of each gene in each sample was independently repeated thrice, and experiments were performed in three biological replicates. The qPCR amplification steps were a holding stage at 95 • C for 5 min, followed by the cycling stage with 40 cycles of 94 • C for 10 s and 60 • C for 40 s. The subsequent melt curve stage consisted of 95 • C for 15 s, 60 • C for 1 min, and 95 • C for 15 s. Small subunit rRNA (SSU rRNA) was used as an internal control. The relative expression of target genes was calculated by comparison with SSU rRNA using the 2 − Ct method. The final results are shown as means ± standard deviations (SD) (Livak and Schmittgen, 2001).

Summary and Assembly Information of Genome Sequencing Data
A total of 4,795,537,118 bp raw data of R. kratochvilovae strain YM25235 were obtained from the Nanopore sequencer. Furthermore, 4,578,684,148 bp clean data were obtained. Sequencing depth was 193.04× and the total length of the assembled complete genome was 23.71 Mb, which is higher than that of R. kratochvilovae strain LS11 (22.11 Mb), R. toruloides   (Figure 1), SwissProt, and TrEMBL databases, respectively. Furthermore, 2,288 genes could be classified into three GO categories: cellular component, biological process, and molecular function (Figure 2A). To gain further insights into the gene function of R. kratochvilovae strain YM25235, we analyzed 3,143 genes and plotted the KEGG annotation classification statistics ( Figure 2B). Summary for number of genes annotated in each database is shown in Supplementary Table 2. In addition, a repeat sequence of 557,107 bp was obtained, accounting for 2.35% of all sequences ( Table 4).

Summary of Transcriptome Sequencing Data
After filtering the raw data, a total of 32.34 and 31.12 million clean reads were obtained for the treatment and control groups, respectively. The percentage of reads mapping to the R. kratochvilovae strain YM25235 genome was 81.1% in the treatment group and 75.9% in the control group. Uniform alignment rates across samples indicate the comparability of data between the samples ( Table 5). The number of predicted novel genes was 40, and a total of 8,048 expressed genes were detected, including 8,008 known genes and 40 predicted novel genes. A total of 5,227 novel transcripts were detected, of which 4,840 belonged to the novel alternatively spliced isoforms of known protein-coding genes, 40 transcripts of novel  protein-coding genes, and the remaining 347 to long noncoding RNAs. For the functional annotation, the reference gene set was aligned with sequences from major databases, including Swiss-Prot, Nr, GO, and KEGG. We performed principal component analysis (PCA) among all the samples (Supplementary Figure 1).

Identification of Differential Gene Expression
Overall, 1,300 genes were significantly determined to be DEGs, which comprised 553 upregulated and 747 downregulated DEGs (Figure 3). They occupied 16.19% of all genes that were identified in this study. The top 10 DEGs with the most dramatic expression changes were listed in Table 6. The top three most upregulated DEGs were histone-binding protein RBBP4, molecular chaperone GrpE, protein of pyridoxal phosphate phosphatase-related family (PLPP). The three most downregulated DEGs were NADH oxidase (Nox), glycosyltransferase family 31 protein (GT31), and esterase/lipase (LIP). In addition, a minority of DEGs were annotated to hypothetical proteins, indicating their unknown functions. The qPCR analysis for all 15 DEGs represents a close correlation (Pearson's correlation coefficients = 0.815, p < 0.01) in fold changes for DEGs between RNA-seq and qPCR (Figure 4), suggesting the accuracy of RNA-seq and DEG analyses.

GO and KEGG Enrichment Analyses of DEGs
All DEGs were divided into three GO categories: biological process (19 terms), molecular function (18 terms), and cellular component (6 terms) (Figure 5). Among the BP terms, the three GO terms containing the highest number of DEGs were metabolic process, cellular process, and cellular component organization or biogenesis. Moreover, some BP terms involving response to low temperature were detected, such as biological process, negative regulation of biological process and growth. In the MF category, catalytic activity, binding, and transporter activity were the top three terms, and for the CC category, cell, cell part, and membrane were the top three terms. In addition, the membrane part ranked fourth.

DISCUSSION
Compared with several Rhodosporidium strains that have been sequenced, R. kratochvilovae strain YM25235 was isolated from Chenghai Lake in Yunnan, China. This plateau lake locates at the border of Tibet Plateau and Yunnan plateau, where the water temperature ranges from 2 to 31 • C, with an average temperature of about 15 • C (Wan et al., 2005). The high-quality genomic and transcriptomic data of R. kratochvilovae strain YM25235 were obtained and the molecular mechanism regarding low temperature response of R. kratochvilovae strain YM25235 was explored in this study. The comparative transcriptomic analysis of R. kratochvilovae strain YM25235 treated by low temperature and control temperature identified several DEGs primarily involved in BP functions (e.g., metabolism, regulation processes, and stress response). It was established that low temperature induced changes in the biological process of R. kratochvilovae strain YM25235. In addition, the terms associated with metabolism were overrepresented (77, 71.96% of all terms) among the enriched KEGG pathways and adjustment of metabolism function due to gene expression changes primarily contributed to the low temperature response of R. kratochvilovae strain YM25235.
In particular, steroids can serve as important components of cell membranes, which can alter the fluidity of cell membranes (Fernández-Cabezón et al., 2018). Moreover, triterpenoid serves as a precursor for steroid biosynthesis (Davis and Croteau, 2000). In this study, several KEGG pathways involved in sesquiterpenoid and triterpenoid biosynthesis were significantly enriched. This evidence suggested that steroids could act as a crucial factor for changes in the protective ability FIGURE 3 | A volcano plot of differentially expressed genes (DEGs). The red spots indicate significant upregulation, whereas blue spots indicate significant downregulation. The black spots signify no significantly changed genes. of R. kratochvilovae strain YM25235 against low temperature by altering the fluidity of cell membranes. Therefore, genes involved in steroid products are useful targets in genetic manipulation to improve the cold resistance of R. kratochvilovae strain YM25235. Another analysis revealed that the DEGs of the biosynthesis of steroid pathway were overwhelmingly upregulated (Supplementary Figure 2), further supporting that R. kratochvilovae strain YM25235 adapts to the low temperature environment by adjusting the biosynthesis of steroids. Ubiquinone and other terpenoid quinone are involved in MVA (Dallner and Sindelar, 2000). The MVA pathway, in turn, is upstream of not only steroid biosynthesis but also carotenoid biosynthesis (Katsuki and Bloch, 1967;Groeneveld, 1999). However, carotenoids could modulate membrane fluidity because carotenoids can stabilize the membrane to respond to low temperature (Chattopadhyay et al., 1997;Jagannadham et al., 2000). Low temperature stress also leads to increased intracellular reactive oxygen species (ROS) (Gasch, 2003), more carotenoids biosynthesis can mitigate cellular damage caused by ROS (Breierová et al., 2018). In this study, ubiquinone and other terpenoid quinone biosynthesis were significantly enriched. Certain amino acids, including alanine, aspartate, glutamate, and phenylalanine metabolic pathways were enriched, suggesting an altered utilization of amino acids by R. kratochvilovae strain YM25235 at low temperature, similar to the findings reported by Beltran et al. (2007). In addition, glutathione metabolism was significantly enriched, glutathione (GSH) protects the cells by neutralizing (i.e., reducing) ROS (Wu et al., 2004), Thus, GSH might protect R. kratochvilovae strain YM25235 from increased ROS at low temperature. In addition, the DNA replication, mismatch repair, and nucleotide excision repair pathways were enriched in the current study. Previous studies have reported genes whose expression is affected by low temperature associated with DNA replication (Driessen et al., 2014). Thus, it can be inferred that the changes occurring in the levels of DNA replication might be involved in the important biological processes in R. kratochvilovae strain YM25235 acting at low temperature. Moreover, most of the downregulated DEGs were FIGURE 6 | Top 20 differential gene pathway enrichment results. The X-axis represents the enrichment factor value, and the Y-axis represents the pathway name. Color represents the Q value, with smaller values representing more significant enrichment results. The size of the dots represents the number of DEGs. Rich factor refers to the enrichment factor value, and a larger value indicates a more pronounced enrichment result.
related to propanoate metabolism, valine, leucine, and isoleucine degradation, pentose, and glucuronate interconversion. We speculate that low temperature limits the progression of these biological responses (Supplementary Figure 2). MAPKs are serine/threonine kinases that, when phosphorylated, enter the nucleus and phosphorylate several transcription factors, enzymes, and other proteins that modulate the cellular activity (Cowan and Storey, 2003). The MAPK gene family is well known for its function in transmitting stress signals from the environment to the cell nucleus (Cowan and Storey, 2003). Significant enrichment of the MAPK pathway in this study indicated the implication of MAPKs in transmitting stress signals from the environment to the cell nucleus in R. kratochvilovae strain YM25235 at low temperature. Traditionally, the HOG pathway, a conserved MAPK cascade, has been considered a specific signaling pathway, responding to changes in external solute concentrations (Hohmann, 2002;Westfall et al., 2004). However, the HOG pathway was demonstrated to function as a transducer of cold stimuli (Panadero et al., 2006). Hog1, as a key gene of the HOG pathway was implicated in response to low temperature stress in our lab study (unpublished). However, Hog1 was not identified as a DEG in this study. We speculate that this process is achieved at the protein level or by certain other mechanisms and is not manifested by differences in the transcript levels.
The top 10 identified DEGs might play a key role for R. kratochvilovae strain YM25235 when grown at low temperature. The RBBP4 gene has been reported to regulate signaling through the MAPK pathway . The activity of GrpE, a nucleotide exchange factor for the heat shock protein DnaK, is downregulated in response to increasing temperature (Bracher and Verghese, 2015). In this study, the expression of GrpE significantly increased at low temperature, indicating that GrpE is a key molecule participating in the response of R. kratochvilovae strain YM25235 to low temperature. The nucleus export protein Brr6 is involved in mRNA and protein export from the nucleus and contribute to mitosis (Lo Presti et al., 2007). The translational activity of mRNAbased protein biosynthesis greatly decreased at low temperatures (Jones and Inouye, 1996). Vitamin B6 as oxygen reactive species scavengers and factors was able to increase resistance to biotic and abiotic stress (Ehrenshaft et al., 1999;Bilski et al., 2000). This vitamin comprises six interconvertible pyridine compounds (vitamers) such as pyridoxal 5'-phosphate (PLP) and five others (ShuoHao et al., 2019). Gene that codes a protein of pyridoxal phosphate phosphatase-related family (PLPP) were significantly upregulated in this study, indicating that PLPP plays a key role in responding to low temperature. In this study, the transcript levels of certain genes involved in regulating the transcription were significantly changed at low temperature. HTH and GATA were detected with significant changes at the transcript level. Most of these proteins were transcription activators and are known to negatively regulate their expression (Baumeister et al., 1992). Thus, the results revealed potential key roles of these genes for R. kratochvilovae strain YM25235 to grow at low temperature. Notably, the functions of several top DEGs remain largely unexplored, and they cannot be linked well to the response of R. kratochvilovae strain YM25235 to low temperature. However, this study is the first to reveal their acute response of R. kratochvilovae strain YM25235 to low temperature.
Membrane rigidity increases when cells are incubated at low temperatures in vitro (Hayashi and Maeda, 2006). Microorganisms overcome the membrane stiffness at low temperatures by adapting to the environment to keep the membrane fluidity constant to a degree (Morgan-Kiss et al., 2006). To accomplish this, cold stress triggers a coordinated response that induces fatty acid desaturases and dehydrogenases, which in turn increase the ratio of polyunsaturated to saturated fatty acids and/or actually decrease the length of fatty acid chains in the membrane (Dickens and Thompson, 1981;Avery et al., 1995;Carty et al., 1999). The PUFAs association with low temperature response in R. kratochvilovae strain YM25235 has been confirmed in our previous study (Cui et al., 2016). However, no pathway and no DEGs related to PUFA biosynthesis were significantly enriched in this study. This is because the response to low temperature stress is a complex process. We speculate that this process may be achieved at the protein level or by certain other mechanisms and is not manifested by differences in the transcript levels.

CONCLUSION
In this study, the genome-wide transcriptional changes of R. kratochvilovae strain YM25235 at low temperature was obtained by whole genome sequencing of R. kratochvilovae strain YM25235 and transcriptome sequencing from low temperature (15 • C) and normal temperature (30 • C) groups. Likewise, this study also showed that there exist genes whose functions are still not known and certain mechanisms of response to low temperature cannot be elucidated. Further experimental verifications are needed. Notably, despite the DEGs obtained in this study are related to low temperature response, further verification for DEGs directly involved in the response is required in the future.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA739038 (genome sequencing and assembly) and PRJNA739250 (Transcriptome).