Triacylglycerols (TAGs) are neutral lipids that are naturally produced via fatty acid biosynthetic pathways. TAGs are used as foods, animal feeds and feedstocks for chemical syntheses1. Recently, derivatives of TAGs, such as long-chain fatty acid methyl esters2 and alkanes3, have been explored as drop-in biofuels. However, traditional technologies that are based on oil crops and plants are limited in terms of their production capacity and rate. The harnessing of microbial fatty acid biosynthetic machineries may circumvent those limitations4 5. Some oleaginous species are exceptionally robust in converting carbohydrates into lipids, such as TAGs6. Yeast species of the anamorphic genus Rhodotorula have often been used for transformation of renewable resources into microbial oils7 8 9. Rhodosporidium toruloides (formerly known as Rhodotorula glutinis or Rhodotorula gracilis), which is a non-pathogenic, pink-coloured basidiomycetous fungus, can accumulate lipids to more than 70% of its dry cell weight6 10. This 'red' yeast is also a good producer of carotenoids11 and biotechnologically important enzymes, such as cephalosporin esterase12 and epoxide hydrolase13. More importantly, R. toruloides has excellent tolerance of inhibitory compounds that are found in biomass hydrolysates14. Thus, R. toruloides is a unique yeast strain of great biotechnological potential. Although a copious literature presenting the results of lipid production by oleaginous species exists6 15, the understanding of the molecular basis of microbial oleaginicity and the sparse capacity for rational strain engineering is limited, largely because genetic backgrounds of naturally oleaginous species remain poorly documented. Next-generation DNA sequencing technologies enable rapid genome decoding, high-throughput transcriptome sequencing (RNA-seq)16 and gene expression profiling (3′-Tag sequencing)17. Furthermore, advances in mass spectrometry-based proteomic analysis allow for the identification of thousands of proteins in a single experiment. The mystery of microbial oleaginicity may be unearthed using a natural lipid producer and an in-depth study that integrates the above-mentioned state-of-the-art technologies18 19 20. Here we present the results of genomic sequencing of R. toruloides, and differential transcriptomic and proteomic analyses of samples that were prepared under various culture conditions. We show that lipid accumulation under nitrogen-limited conditions is tightly connected with cellular processes that are related to lipogenesis, macromolecule metabolism and autophagy. We also find novel genes and pathways (such as fatty acid synthesis and lipid storage) that are committed to oleaginicity. Much information is consistent at different omic levels, and such multi-omic data significantly enriches our understanding of the molecular basis of lipid accumulation. This should facilitate the engineering of a natural lipid producers and other microorganisms for improved efficiency, robustness and economics of the microbial production of fatty acid derivatives. Results Genomic and transcriptomic sequencing and annotation The genome of the haploidic strain R. toruloides NP11 was sequenced using the Illumina sequencer platform (Supplementary Table S1), and the de novo sequence assembly resulted in a 20.2-Mb genome with an average coverage of 96× (Supplementary Fig. S1) and 94 scaffolds (>1 kb, L50=575 kb, Supplementary Table S2). After PCR-based gap closure, the genome was narrowed down to 34 scaffolds, which included one ribosomal DNA unit sequence (Supplementary Fig. S2 and Supplementary Table S2). The genome was organized into 16 chromosomes that ranged from 650 to 1900 kb (Supplementary Fig. S3 and Supplementary Methods). To validate the completeness of the genome assembly and to define the transcriptional landscape of R toruloides, the transcriptome was sequenced using RNA-seq, which generated 6.7 M paired-end reads (Supplementary Table S1). It was found that 79% of the sequenced complementary DNA (cDNA) fragments and 98.6% of the assembled expressed sequence tags could be mapped to the genome. Moreover, except for KOG1185 and KOG1394, which are absent in basidiomycetous fungi, 246 of the 248 core eukaryotic genes21 were found. These results indicated that the assembled genome encompassed a very high percentage of the gene space. The repetitive sequences constituted 1.66% of the current assembly (Supplementary Table S3), and RNA-seq alignment combined with de novo gene prediction, homology search and manual curation resulted in the prediction of 8,171 protein-coding genes (Supplementary Fig. S4 and Supplementary Table S4). The gene-coding regions were considerably dense with an average of 2.6 kb per gene. In addition, similar to other basidiomycetous fungi, split genes, which have short and multiple introns, were predominant, and 7,938 genes were found to have at least one intron. Moreover, most genes with 5′- and 3′-untranslated regions were annotated with the assistance of RNA-seq (Supplementary Table S4). It was found that 90% of the gene models had more than 30% of their length covered by RNA-seq reads, which were mapped to the annotated open-reading frames using the Bowtie software22, and 38% of the annotated proteins had at least one unique peptide that was present at the proteomic level (Supplementary Table S5). RNA-seq data indicated that 1,371 genes encoded two or more transcript isoforms, suggesting that the transcriptome was more complex than that of the static genomic data set (Supplementary Fig. S5). For the functional annotation, 79% of the R. toruloides proteins had homologous proteins or conserved domains of known proteins (Supplementary Table S5). Approximately 60% of proteins were annotated with Gene Ontology (GO) terms, and 32% of proteins were found with KEGG (Kyoto Encyclopedia of Genes and Genomes) orthology numbers. On the basis of these annotations, we reconstructed the central and lipid metabolism of this organism (Fig. 1). In comparison with Saccharomyces cerevisiae, R. toruloides contains ATP:citrate lyase (ACL), and the mitochondrial β-oxidation (MBO) and carotenoid biosynthetic pathways. Although cytoplasmic acetyl-CoA synthase can produce acetyl-CoA from acetate, in oleaginous fungi and higher eukaryotes, acetyl-CoA, which results from the cleavage of citrate by ACL, is a major source for lipid synthesis23. Unlike β-oxidation in peroxisomes, the MBO pathway recovers more energy through the degradation of fatty acids by using the flavin adenine dinucleotide (FAD) cofactor, which couples with the respiratory chain to produce ATP. Moreover, MBO-associated enzymes are required to metabolize branched-chain amino acids to produce acetyl-CoA24 25. Thus, MBO may provide alternative acetyl-CoA and energy sources for fatty acid synthesis by degrading amino acids and membrane lipids (vide infra). Two genes that code putative carotenoid synthesis-related enzymes, phytoene synthase (PSY1) and phytoene dehydrogenase (CRTI), were identified, revealing the genetic basis for the formation of the pink-coloured pigments by this organism (Supplementary Fig. S6 a). A novel fatty acid synthase (FAS) system was found in R. toruloides (Fig. 2). FASs are typically classified into two variants, the dissociated type II system and the integrated type I multi-enzyme. In fungi, FASs are usually composed of eight distinct domains and organized into two subunits or one polypeptide26. Although split FASs are common in fungi, FAS in R. toruloides is split into a novel form that consists of two subunits: the β-subunit (Fas1), which contains acyltransferase and enoyl reductase domains, and the α-subunit (Fas2), which includes all other domains (Fig. 2a). Two acyl carrier protein (ACP) domains of 76% sequence identity are present in Fas2, and these two ACPs share high similarity with ACPs from other species (Supplementary Figure S7). Unique peptides from different ACPs were also found by LC-MS/MS (liquid chromatography–mass spectrometry) during proteomic analysis (Fig. 2b; Supplementary Fig. S8), which unambiguously indicated the expression of the FAS with two ACPs. As ACP has a key role in iterative substrate shuttling in the chain-elongation process, the presence of tandem ACPs may improve the efficacy of fatty acid biosynthesis by providing a higher intermediate concentration in the reaction chamber of FAS27. Transcriptomic analysis Early biochemical studies demonstrated that nitrogen limitation was a major regulator for initiating lipid overproduction (Supplementary Fig. S6 d,f,g)6 10. In oleaginous fungi, mitochondrial NAD+-dependent isocitrate dehydrogenase (IDH) is dependent on adenosine monophosphate (AMP), and under nitrogen limitation conditions the cellular AMP level is lowered because of AMP deaminase activity, which leads to an impaired activity of IDH. Consequently, mitochondrial citrate accumulates and permeates out from the mitochondria to the cytosol where it serves as a substrate for ACL to form acetyl-CoA for fatty acid biosynthesis. In addition, cytosolic NADP+-dependent malic enzyme (ME) is considered as the key enzyme that supplements nicotinamide adenine dinucleotide phosphate (NADPH) for de novo lipogenesis. In addition to lipid accumulation, nitrogen limitation is expected to mediate a more global response, such as the release of nitrogen catabolite repression (NCR)28 and the induction of autophagy by the target of rapamycin (TOR) signalling cascades29 30. We grew chemostat cultures of R. toruloides to examine the influence of nitrogen starvation on this organism (Table 1). In a minimal medium with abundant ammonium (MM), cells accumulated lipids to 22.8%, and the residual ammonium concentration was 47±1 mM. In minimal medium with trivial amounts of ammonium (MM-N), the cellular lipid content was 33.3%. Furthermore, the oil yield, which was based on consumed sugar (Y SP), in the MM-N medium was 3.6-fold higher than that in the MM medium. These results clearly indicated that nitrogen limitation promoted lipid accumulation and that the metabolic flux was repartitioned for lipid biosynthesis. The 3′-tag digital gene expression (DGE) profiling experiment produced more than three million tags from each RNA sample that was prepared from cells cultured in the MM or MM-N media (Supplementary Table S1). Of those tags, 78% from the MM sample and 63% from the MM-N sample mapped to annotated transcripts, and 89% from the MM sample and 74% from the MM-N sample mapped to the genome (Supplementary Table S6). In addition, 5.7% of those tags mapped to antisense transcripts of predicted genes, which led to the observation that 74% of the total genes exhibited antisense transcription (Supplementary Table S6). These data indicated that antisense transcription is widespread in R. toruloides. Approximately 12% of the reads that uniquely mapped to sense transcripts were not aligned to the 3′-most NlaIII site (Supplementary Figure S9), which suggested that transcription at the 3′-end was heterogeneous. For transcript quantification, the tags that mapped to sense sequences were counted and compared. Of the 7,045 genes that were expressed in these two samples (Supplementary Table S6), 1,177 and 853 genes were significantly upregulated and downregulated, respectively (Fig. 3a). It was obvious that nitrogen limitation influenced a variety of cellular processes such as metabolism, localization and transportation (Fig. 3b). NCR-related genes were induced (Fig. 3c), and MEP2, which is the orthologue of the S. cerevisiae high-affinity ammonia transporter, was upregulated tenfold. In addition, three genes that coded amino acid permeases (RHTO_00398, RHTO_00931 and RHTO_01882) were upregulated. The central nitrogen metabolism machinery, which is composed of glutamate dehydrogenase (GDH), glutamine synthase (GLN) and glutamate synthase (GOGAT/GLT), is the hub of cellular nitrogen utilization, and it was discovered that the transcription of GDH1, GDH2 and GLN1 were upregulated under nitrogen-limited conditions. It was also found that the expression of genes involved in urate and urea degradation was upregulated, whereas the expression of genes related to the utilization of other poor nitrogen sources was not significantly altered (Supplementary Table S7). Furthermore, genes coding proteases (PRB1 and PEP4) and amino acid permeases (RHTO_05155 and RHTO_02128) that are linked to vacuolar protein degradation were induced. As Pep4 activates other vacuolar hydrolases, this elevated expression level would enhance the degradation of macromolecules and organelles. Indeed, autophagy-related genes (ATG1, ATG2, ATG20 and ATG8) and vacuolar-type ATPases (VMA5, VMA9 and VMA11), vacuolar lipase (ATG15) and α-mannosidase (AMS1) were also upregulated (Supplementary Table S7), which suggested that the cells initiated the turnover and recycling of intracellular components. The downregulated genes were related to primary metabolism, macromolecule metabolism and cellular metabolism, especially nucleobase, nucleoside, nucleotide and nucleic acid metabolic processes (GO: 0006139) and ribonucleoprotein complex biogenesis (GO: 0022613). Ribosome biogenesis is responsive to nutrient starvation through an evolutionarily conserved TOR signalling mechanism31, and for the 65 genes that were sampled by DGE and coded cytoplasmic ribosomal proteins, 39 genes were downregulated (Fig. 3d, blue), which suggested that inhibition of the TOR complex suppressed the transcription of ribosomal components. Together, nitrogen limitation triggered a set of rescue pathways that included the induction of nitrogenous compound uptake, nitrogen assimilation, alternative nitrogen-source utilization, vacuolar protein degradation and autophagy, and inhibition of nitrogenous macromolecule biosynthesis. The regulation of central metabolism and lipid metabolism at the transcriptional level was also investigated (Fig. 1; Supplementary Table S8). FAS in the MM-N sample was upregulated by 4.7- and 1.8-fold for FAS1 and FAS2, respectively, but there were no significant differences in the expression of genes that were involved in the late stage of TAG biosynthesis. Isoprenoids are a variety of compounds (for example, steroids, carotenoids and terpenes) that are derived from the common precursor dimethylallyl diphosphate and its isomer 3-isopentenyl diphosphate. Four genes that are involved in isoprenoid backbone synthesis were upregulated in the MM-N sample (Fig. 1). Indeed, the pink colour development of the culture was an indication of carotenoid biosynthesis. As lipids are highly reduced compounds, the supplement of a reducing power, mainly NADPH, is equally important to that of the precursor acetyl-CoA6. Activities of glucose-6-phosphate dehydrogenase (ZWF1), 6-phosphogluconate dehydrogenase (GND1), NADP+-dependent isocitrate dehydrogenase (IDP1) and ME (ME1) can generate NADPH. Although ZWF1 and IDP1 were upregulated by 1.9- and 4.8-fold, respectively, the expression of GND1 and ME1 were downregulated by 6.8- and 4.5-fold, respectively. The reduced transcriptional level of ME1 was consistent with early observations that indicated the expression of ME was suppressed during nitrogen exhaustion32 33. Interestingly, the expression of some genes that were involved in lipid degradation (that is, lipolysis, the β-oxidation cycle, and glyoxylate shunt) were also upregulated (Fig. 1), which might be caused by free fatty acids being released by an elevated autophagy process34. We also compared the highly expressed genes of cells that were cultured in yeast-extract peptone dextrose (YPD), MM or MM-N media (Fig. 4; Supplementary Table S10). It was found that the numbers of highly expressed genes were similar, and approximately half of those genes were shared by the three samples (Fig. 4b, group G). The shared genes included those coding ribosomal proteins, actin, ACP (RHTO_01048), TAL1 of the phosphopentose pathway, MDH1 and IDH2 of the tricarboxylic acid (TCA) cycle and GAPDH, ENO1 and FBA1 of glycolysis. Forty-two genes were highly expressed exclusively in the MM-N sample (group C) and included those encoding FAS1& 2, PYC1, GDH1 and acyl-CoA-binding protein (RHTO_03853). Genes coding ACL1 and pyruvate dehydrogenase components (LAT1 and LPD1) were highly expressed in the MM and MM-N samples (group E), where lipid contents were greater than 20%. Comparative proteomic analysis In an earlier study, we performed comparative proteomic analysis of R. toruloides based on the proteome database of S. cerevisiae. As a result, 184 proteins were identified, and 46 of those proteins were found to change significantly during the lipid-production stage35. Here we re-analysed raw LC-MS/MS data that was based on the newly annotated protein data set that had been generated from the genomic data of R. toruloides. Although the previous analysis only cumulatively observed 184 proteins, in this study we observed 3,108 proteins, which amounted to roughly a 17-fold increment over the previous data set (Supplementary Table S11). Thus, 38% of annotated genes were confirmed at the proteomic level. We were able to identify 2,057 proteins, which included over 32-fold more proteins than our previous data set (Fig. 5a; Supplementary Table S11). The intersection between the '24 h' and the '96 h' samples included approximately 100 more proteins than those of the other two intersections, which indicated that samples at the lipid-production stages more physiologically resembled each other. Further semi-quantitative analysis revealed 538 proteins that significantly varied in abundance between any two samples (Fig. 5b and Supplementary Data 1). A good correlation was observed between protein expression levels of the '24 h' and the '96 h' samples, which had a Pearson's correlation coefficient of 0.79 (Fig. 5c) and suggested that many proteins were constantly expressed during lipid accumulation under nitrogen-limited conditions. Significantly higher levels were found for proteins linked to nitrogen metabolism under nitrogen-limited conditions. Specifically, proteins involved in nitrogenous compound transportation (Mep2, RHTO_07825, RHTO_07915, RHTO_00398, Dur3 and Nrt2), alternative nitrogen source utilization (proline, urate, urea and nitrate) and central nitrogen metabolism (Gdh1, Gln1, Gln2 and Glt1) were found at higher levels in the lipid-production samples than those in the 'Seed' sample (Fig. 5d), which suggested that nitrogen limitation alleviated the NCR process. On the other hand, the levels of ribosomal proteins, aminoacyl transfer RNA synthases, translation initiation and elongation factors, and protein folding-related proteins were reduced (Supplementary Data 1), which indicated that the protein biosynthesis machinery was suppressed. The increased levels of vacuolar degradation- and autophagy-related proteins, such as vacuolar proteases (Cpb1, Rpb1 and Pep4), vacuolar-type ATPases (Vma2, Vma4, Vma13 and Stv1) and α-mannosidase (Ams1) suggested that an activated autophagy process was involved in responding to nitrogen starvation. The above-mentioned proteomic observations were also in good agreement with the differential transcriptomic data. Proteins involved in lipogenesis, especially for de novo fatty acid biosynthesis, were found at elevated levels (Fig. 5e). These proteins included Acl1, Me1, Acc1, Fas1 and Fas2 for fatty acid synthesis and glycerol-3-phosphate dehydrogenase isoform 2 (Gpd2) and glycerol-3-phosphate acyltransferase (Gat1) for glycerolipid synthesis. Cytosolic NADP+-dependent ME is the key enzyme that supplements NADPH for de novo lipogenesis6 32 33. Interestingly, although the transcription of ME1 was downregulated, the protein level of Me1 was significantly higher in the lipid production samples, which suggested that the regulation of ME activity was complicated. The level of acetyl-CoA synthase (Acs1) was increased by 2.2- and 2.1-fold in the '24 h' and '96 h' samples, respectively (Fig. 6a); however, levels of acetyl-CoA C-acetyltransferase (Erg10) and acetyl-CoA hydrolase (Ach1), which are the branch point enzymes for acetyl-CoA, were decreased. These data demonstrated that the acetyl-CoA flux was switched to the synthesis of fatty acids but not of isoprenoids. In fungi, a perilipin-like protein was characterized as a lipid-droplet protein that protected lipids from degradation36. RHTO_05627 was a perilipin family protein37 that showed significantly higher levels at lipid-production stages, and the levels of proteins that were related to MBO were also increased, which was consistent with the upregulation of β-oxidation at the transcription level. The levels of proteins associated with the TCA cycle (that is, Idh1, Sdh1, Sdh2 and Mdh1) decreased in the lipid-production samples (Fig. 5f), which was in agreement with early biochemical observations6, and the levels of proteins involved in pyruvate metabolism were significantly altered (Figs 5f and 6b). In the lipid-production samples, higher levels of pyruvate kinase (Pyk1) and pyruvate dehydrogenase components (Pdb1 and Lpd1) were found, whereas a lower level of pyruvate decarboxylase (Pdc1) was observed, which suggested that pyruvate was channelled to acetyl-CoA. Pyruvate carboxylase (PYC) functions as a key enzyme for the anaplerosis of the TCA cycle and gluconeogenesis. PYC was also proposed to be an essential part of the transhydrogenation machinery, which is composed of PYC, ME and malate dehydrogenase for the conversion of NADH to NADPH6. As the level of phosphoenolpyruvate carboxykinase (Pck1) that is responsible for the transformation of pyruvate to phosphoenolpyruvate was reduced, an increased level of Pyc1 should be mainly related to the NADPH supply for lipogenesis (Fig. 6b). Interestingly, enzymes catalysing these irreversible reactions in glycolysis (Hxk2 and Pyk1) were elevated. Moreover, the key enzyme in glycolysis, phosphofructokinase (Pfk1), was also augmented by 1.4-fold (P-value 2, a P-value 1.35 for the G-test (corresponding to 3.9 for raw spectral counts) of each of the two samples were assigned as differentially expressed proteins, and this analysis led to 538 proteins that showed significant level changes. Author contributions Z.K.Z. conceived the project. Z.K.Z. and Z.Z. designed the study. Z.Z. performed most of the experiments and analyses unless specified otherwise. S.Z. isolated the haploid strains and participated in sample preparation, genome annotation and data analysis; H.L., M.Y. and H.Z. participated in the proteomic data analysis; H.S. prepared the chemostat cultures; X.L. and F.Y. participated in the genome annotation; Y.Z. participated in data analysis; G.J. prepared the scanned electrical microscope images; and Z.Z. and Z.K.Z. wrote the manuscript. Additional information Accession codes: This Whole Genome Shotgun Sequencing project has been deposited into DDBJ/EMBL/GenBank under the accession number ALAU00000000. The version that is described in this paper is the first version, ALAU01000000. Raw reads of the WGS sequencing and RNA-seq analysis have been deposited into the NCBI Sequence Read Archive (SRX156633–SRX156637). DGE profiling data have been deposited into NCBI's Gene Expression Omnibus and was given the accession number GSE39023. The proteomic data that is associated with this manuscript may be downloaded from ProteomeCommons.org Tranche using the following hash: HzAX+3G+zoH9LngqJ+29NyszSltcjGUyxKWjNhCIRKOab/G0LG9Kw15quSpAteJOio7+V2tR9O3vt9pyHyu1O9/I54gAAAAAAAA43A==. All data are also available at http://www.bioconversion.dicp.ac.cn/EWEB/ZZW-data.html. How to cite this article: Zhu, Z et al. A multi-omic map of the lipid-producing yeast Rhodosporidium toruloides. Nat. Commun. 3:1112 doi: 10.1038/ncomms2112 (2012). Supplementary Material Supplementary Information Supplementary Figures S1-S14, Supplementary Tables S1-S11, Supplementary Methods and Supplementary References Supplementary Data 1 Differentially expressed proteins between "seed", "24 h" and "96 h"
