Transcriptome analysis reveals the genetic foundation for the dynamics of starch and lipid production in Ettlia oleoabundans

(1)

Contents lists available atScienceDirect

Algal Research

journal homepage:www.elsevier.com/locate/algal

Transcriptome analysis reveals the genetic foundation for the dynamics of starch and lipid production in Ettlia oleoabundans

Mark H.J. Sturme

^a,⁎,1

, Yanhai Gong

^b,1

, Josué Miguel Heinrich

^a

, Anne J. Klok

^a

, Gerrit Eggink

^a

, Dongmei Wang

^b

, Jian Xu

^b

, Rene H. Wij ﬀ els

^a,c

aDepartment of Bioprocess Engineering/AlgaePARC, Wageningen University, The Netherlands

bSingle-Cell Center, CAS Key Laboratory of Biofuels and Shandong Key Laboratory of Energy Genetics, Qingdao Institute of BioEnergy and Bioprocess Technology, Chinese Academy of Sciences, Qingdao, Shandong 266101, China

cFaculty of Biosciences and Aquaculture, Nord University, N-8049 Bodø, Norway

A R T I C L E I N F O

Keywords:

Microalgae Nitrogen starvation Starch

Lipids

Carbon-partitioning Transcriptome

A B S T R A C T

The oleaginous microalgaEttlia oleoabundansaccumulates both starch and lipids to high levels under stress conditions such as nitrogen starvation (N−). To steer biosynthesis towards starch or lipids only, it is important to understand the regulatory mechanisms involved. Here physiological and transcriptional changes under nitrogen starvation were analysed in controlledﬂat-panel photobioreactors at both short and long time-scales.

Starch accumulation was transient and occurred rapidly within 24 h upon starvation, while lipid accumulation was gradual and reached a maximum after 4 days. The major fraction of accumulated lipids was composed of de novo synthesized neutral lipids - triacylglycerides (TAG) - and was characterized by a decreased composition of the polyunsaturated fatty acids (PUFAs) C18:3 and C16:3 and an increased composition of the mono-unsaturated (MUFAs) and saturated (SFAs) fatty acids C18:1/C16:1 and C18:0/C16:0, respectively. RNA-sequencing revealed that starch biosynthesis and degradation genes show diﬀerent expression dynamics from lipid biosynthesis ones.

An immediate rapid increase in starch synthetic transcripts was followed by an increase in starch degrading transcripts and a decrease in the starch synthetic ones. In contrast, increased gene expression for fatty acid and TAG synthesis was initiated later and occurred more gradually. Expression of several fatty acid desaturase (FAD) genes was decreased upon starvation, which corresponds to the observed changes to higher levels of MUFAs and SFAs. Moreover, several homologs of transcription regulators that were implicated in controlling starch and lipid metabolism in other microalgae showed diﬀerential gene expression and might be key regulators of starch and lipid metabolism inE. oleoabundansas well. Our data provide insights into the genetic foundation of starch and lipid metabolism inE. oleoabundansunder nitrogen starvation and should facilitate metabolic engineering towards tailored strains with desired storage compound composition.

1. Introduction

With a growing world population and declining natural oil and gas reserves, the global demand for sustainable and renewable resources becomes ever more relevant. In this respect, microalgae are a promising source for sustainable production of food additives, chemical building blocks and cosmetic ingredients. Microalgae can accumulate high levels of storage compounds such as carbohydrates, lipids and pigments under stress conditions and can serve as a renewable source of proteins. They can be cultivated on marginal lands and in (semi-)arid regions and outdoor cultivation of marine and halotolerant microalgae will require a lower freshwater input. Thereby they compete less with agricultural

crop production and have a smaller eﬀect on drinking water supply.

Moreover, they are a particularly good substitute for some vegetable oils since the negative environmental impact of e.g. palm, one of the main oil sources currently used in food and cosmetics applications, can potentially be mitigated [1–3]. Nevertheless, the market for most microalgae-derived products currently is not economically competitive with those derived from existing plant-based and petrochemical resources. To overcome this gap, microalgal strain development in com- bination with improved bioreﬁneries is required [4,5]. Many studies for microalgal strain development have focused on metabolic engineering by the introduction or mutagenesis of single target genes involved in lipid, starch or pigment biosynthesis [6,7] and less on engineering of

https://doi.org/10.1016/j.algal.2018.05.004

Received 13 December 2017; Received in revised form 11 April 2018; Accepted 4 May 2018

⁎Corresponding author.

1Equal contribution.

E-mail addresses:[email protected](M.H.J. Sturme),rene.wijﬀ[email protected](R.H. Wijﬀels).

Available online 18 May 2018

T

(2)

multiple genes or complete pathways [8]. The selected gene targets for metabolic engineering have been identiﬁed by transcriptomics or were based on detailed biochemical and genetic information from other microalgal or plant models. Most of these studies focused on microalgae that accumulate only lipids or starch under stress conditions such as nitrogen or phosphorus starvation, while only a few have focused on transcriptomics and metabolic engineering of so-called“hybrid”producers: oleaginous microalgae that are capable of simultaneous accumulation of starch and neutral lipids such as triacylglycerides (TAG).

The oleaginous microalgaeEttlia oleoabundansandAcutodesmus obliquus and severalChlorellaspecies are known hybrid producers of industrial relevance that accumulate high levels of both starch and lipids under nitrogen starvation [9–11]. To achieve economically viable production of these hybrid producers multi-product bioreﬁneries and cultivation optimization are important factors, but for many industrial applications it is also desired that carbon-partitioning can be altered to produce mainly TAG or starch and to genetically control the switch between starch and lipid metabolism. The genetic foundation for this dual accumulation strategy however is unknown in most of these microalgae, due to a lack of well-annotated genome sequences and a limited number of transcriptome studies. It is also plausible that the regulatory mechanisms for storage compound accumulation in hybrid producers are more complex than in microalgae mainly accumulating a single storage compound.E. oleoabundansis a hybrid producer with several properties that make it a very suitable candidate for industrial production. It can accumulate up to 56% of its cell dry weight as lipids under nitrogen starvation, mainly in the form of TAG, and while isolated as a freshwater strain it can also grow and accumulate storage compounds under saline conditions and high pH as well [12–15]. Furthermore, continuous TAG production in continuously growing cells has been reported [16]. However, the genomic foundation for these favourable industrial traits ofE. oleoabundanscurrently is not yet established. The few studies that addressed the transcriptional changes upon nitrogen starvation in this species, have analysed either a single timepoint at an early-stage (< 24 h) [15] or a late-stage (11 days) during batch N- starvation [17] or investigated transcriptomes during nitrogen limita- tion under turbidostat cultivation [18]. For a full understanding of the cellular and transcriptional changes upon N-starvation it is necessary to track the temporal dynamics of metabolite concentrations, physiological parameters and transcript abundances over both a short and long timespan, to identify key regulatory genes controlling metabolic switch- points. Some recent studies in e.g.Chlamydomonas reinhardtii,Chlorella spp., Nannochloropsis spp. and Monoraphidium neglectum have performed this correlation analysis of storage compound accumulation and transcriptome patterns, which indicated possible transcription factors and metabolic nodes that are involved in or even control the“switch”

between starch and lipid metabolism [19–24]. In addition, integration of transcriptomics with genome-scale metabolic models could aid in pinpointing the critical metabolic nodes that should be targeted in metabolic engineering.

In this study, we therefore set out to perform an in-depth analysis of the temporal dynamics of storage compound accumulation and transcriptome changes during nitrogen starvation inE.oleoabundans, with the aim to identify the metabolic genes and transcriptional factors that are involved in the switch between starch and lipid metabolism in this microalga. Our results revealed the diﬀerences in temporal dynamics of gene expression for the starch and lipid pathways upon nitrogen starvation, and identiﬁed transcriptional regulators that might be involved.

These insights are valuable for rational metabolic engineering in this important industrial microalga.

2. Materials & methods

2.1. Culture conditions and nitrogen starvation experiment

Ettlia oleoabundans UTEX 1185 (culture collection of Algae,

University of Texas, Austin) was grown in Bold's Basal Medium (BBM) at pH 7.5 and 25 °C, with 25.2 mM KNO3 (N-replete conditions) or without nitrate (N-starvation). Culture purity during experiments was monitored by microscopical observation and analysis on a Multisizer™3 Coulter Counter®(Beckman Coulter). Experiments were performed in flat panel air-lift photobioreactors (Labfors 5 Lux, Infors HT, Switzerland) with a working volume of 1.7 L, at pH 7.5 and 25 °C, and sparged with 2% CO2at an airflow of 1.2 L min⁻¹ and photonflux density (PFD) of 800μmol s⁻¹m⁻²at a 12:12 light:dark cycle.

The pH was maintained at pH 7.5 by the on-demand supply of 0.1 M HCl. A nitrogen-replete batch pre-culture ofE. oleoabundans(250 mL) in exponential growth phase was used to inoculate all three photobioreactors, and microalgae were further batch cultivated until a biomass concentration of 2 g cell dry weight (CDW)/L was reached.

Subsequently, continuous nitrogen-replete cultivation (N+) was con- tinued for several days until a stable biomass concentration was maintained. For the nitrogen starvation experiment biomass from the nitrogen-replete continuous culture was collected, spun down, rinsed and resuspended in nitrogen-free BBM medium. The reactors were rinsed with distilled water,ﬁlled with nitrogen-free BBM medium and inoculated with the same volume of the collected microalgae. During the nitrogen-depleted phase (N−), the reactors were run in non-steady state continuous mode. Medium was fed continuously at 800 mL/d to compensate for the volume-loss due to sampling. Samples for analysis of biomass and cell parameters from the nitrogen-replete continuous phase (N+) were collected at four timepoints in the light-phase preceding nitrogen starvation (−25, −23, −19 and −15 h; Fig. 1).

Samples for biomass, cell parameter and molecular analysis from the onset of nitrogen starvation (N−) were collected at time-points at the start, middle and end of each light-phase and at the end of each dark- phase over a 4-day period (0, 2, 6, 10, 23, 29, 33, 46, 53, 71, 77 and 99 h;Fig. 1). Samples were centrifuged, and cell pellets snap frozen in liquid nitrogen and stored at−80 °C. The experiments were performed as three biological replicates.

2.2. Biomass analysis

2.2.1. Cell dry weight, cell count and cell parameters

Cell dry weight (CDW) concentrations were determined as described by Kliphuis et al. [25], by filtering culture broth (around 10 mg of biomass) through pre-dried (100 °C overnight) and pre-weight Whatman glassfibre filter paper (GF/F; Whatman International Ltd., Maidstone, UK). The filter was washed twice with filtered deminer- alized water to remove adhering inorganic salts and trace elements and subsequently dried overnight at 100 °C before weighing. Cell numbers were determined using a Multisizer™ 3 Coulter Counter® (Beckman Coulter). Quantum yield (QY) was determined using an Aquapen-C (Photon Systems Instruments, Czech Republic).

2.2.2. Total fatty acids

Extraction and quantiﬁcation of total fatty acids (TFA) were adapted from Breuer et al. [26]. Around 20 mg of pellet was transferred to bead beating tubes (Lysing Matrix E; MP Biomedicals, Santa Ana, CA, USA) and lyophilized overnight. Freeze-dried cells were disrupted by a bead beating step in a Precellys®24 bead beater (Bertin Technologies) for 2× (3 × 60 s with 120 s pause in between at 2500 rpm), followed by 1× (2 × 60 s with 120 s pause in between at 2500 rpm) in the presence of a chloroform:methanol mixture (1:1.25) to extract the lipids from the biomass. The internal standards C15:0 (tripentadecanoin - T4257;

Sigma-Aldrich, St Louis, MO, USA) and C19:0 (trinonadecanoin - T4632; Sigma-Aldrich, St Louis, MO, USA) were added to the extraction mixture to enable fatty acid quantiﬁcation. Methylation of the fatty acids to fatty acid methyl esters (FAMEs) and the quantiﬁcation of the FAMEs by GC/MS analysis were performed as described by Breuer et al.

[26]. TFA concentration was calculated as the sum of triacylglycerides (TAG) and polar lipids (PL).

(3)

2.2.3. Triacylglycerides and polar lipids

Quantiﬁcation of TAG was similar to the TFA analysis method with some modiﬁcations. After TFA extraction, the chloroform-methanol mixture was evaporated under N2gas and the TFA fraction dissolved in 1 mL hexane and separated based on polarity using a Sep-Pak Vac silica cartridge (6 cm³, 1000 mg; Waters, Milford, MA, USA) equilibrated with hexane. The neutral TAG fraction was eluted with 5 column volumes (10 mL) of hexane-diethyl ether (7:1 v/v). Subsequently the polar lipid fraction containing glycolipids and phospholipids was eluted from the silica cartridge with 5 column volumes (10 mL) of methanol:acetone:- hexane (2:2:1 v/v/v). All solvents were then evaporated from the fractions under N2gas. TAG and PL fractions were then methylated and analysed as described in the TFA analysis section.

2.2.4. Total carbohydrates

Total carbohydrates were extracted and analysed by a colorimetric assay according to DuBois et al. [27]. A biomass dry-weight of∼20 mg was hydrolysed in 5 mL 2.5 M HCl in a heat block at 100 °C for 3 h.

Samples were neutralized using 5 mL of 2.5 M NaOH. Then 0.5 mL of 5% w/w phenol and 2.5 mL of concentrated sulfuric acid was added to 0.5 mL of hydrolysed sample. The samples were incubated at 35 °C for 30 min before reading of the absorbance at 485 nm against a blank of 0.5 mL 5% w/w phenol, 2.5 mL concentrated sulfuric acid and 0.5 mL of deionized water. Glucose was used as a standard.

2.2.5. Starch

Starch content was analysed by enzymatic degradation of starch to glucose using the thermostable α-amylase and amyloglucosidase enzymes from the Total Starch Assay Kit (AA/AMG) (Megazyme International, Ireland) using the protocol described by de Winter et al.

[28]. Around 10 mg of freeze-dried cells were disrupted by bead beating in bead beating tubes (Lysing Matrix E; MP Biomedicals) in the presence of 80% ethanol. Starch was then converted to glucose using the α-amylase and amyloglucosidase enzymes from the kit. Subse- quently glucose was oxidised to gluconate, with the release of 1 mol of hydrogen peroxide (H2O2) that can be measured quantitatively by a

colourimetric assay, by using peroxidase and measuring the production of the dye quinoneimine. Absorbance for all samples was measured against aD-glucose calibration control series at a wavelength of 510 nm.

2.2.6. Protein

Protein concentration was determined on 5 and 10 mg of freeze- dried cells using the Lowry DC protein assay (BioRad).

2.2.7. Biomass productivity/production rates

Biomass productivity,rx, was computed according to Eq.(1). In Eq.

(1),Δx is the biomass concentration changes in the PBR during the periodΔt; x is the averaged biomass concentration (considering the initial andﬁnal biomass concentration of each period);QAis the media ﬂow rate,VRis the reactor volume; and,rxis the actual biomass productivity.

= + ≈ = +

dx dt

Q

V x r x

t Q

V x r Δ

Δ

A R

x A

R

x (1)

The productivities of each metabolite were calculated from Eq.(1), considering their intracellular concentrations.

2.3. RNA isolation, RNA-Seq library preparation and sequencing Microalgal cells for RNA extraction were collected from photobioreactors and harvested by centrifugation. The supernatant was dis- carded and cell pellets immediately snap-frozen in liquid nitrogen and stored at−80 °C until further processing. Total RNA was extracted from frozen microalgal cells (~100–200 mg wet weight) by cell disruption by bead-beating in bead-beater tubes (Lysing Matrix E; MP Biomedicals, Santa Ana, CA, USA) for 1 × 10 s at 5500 rpm in a Precellys®24 bead beater (Bertin Technologies), in the presence of RLT buffer with β- mercaptoethanol (Qiagen). Subsequently, tubes were placed on ice and spun down for 2 min at 13.000 ×g. RNA was further cleaned using the RNeasy Plant Mini Kit (Qiagen). In short, supernatant was transferred to a Qiashredder column, and columns centrifuged for 2 min at 13.000 ×g.Supernatant was transferred to a 2 mL Eppendorf vial, 0.5 V Fig. 1.Experimental design for analysis of the temporal dynamics of cellular content, storage compound accumulation and transcriptional response in nitrogen- replete (N+) and nitrogen-starvation (N−) conditions inE. oleoabundans. Nitrogen starvation was initiated att= 0 h.E. oleoabundanswas cultivated in triplicate in flat-panel air-lift photobioreactors at a photonflux density (PFD) of 800μmol m⁻²s⁻¹and a 12:12 light:dark (LD) cycle. Alternating 12 h light phases (white background) and 12 h dark phases (grey background) are indicated. Individual sampling time points for biomass analysis (open circles) and RNA-seq analysis (open diamonds) are depicted and were taken in the light phase (marker at top of light phase) or in the dark phase (marker at bottom of dark phase).

(4)

of 96% ethanol was added and mixed. The mixture was then processed on Qiagen RNeasy columns, including on-column DNaseI treatment, according to the manufacturers' instructions. Further preparation of RNA and RNA sequencing was performed by BaseClear (Leiden, The Netherlands).

RNA quality control and quantiﬁcation were performed using a BioAnalyzer and an Illumina TruSeq strand-speciﬁc mRNA library was prepared (average mRNA library size of 345 bp) for the construction of the reference transcriptome using 125 bp paired-end sequencing (PE 125) from a pool of all RNA samples (N-replete and N-deplete).

Individual sample libraries were prepared by 50 bp single-end read sequencing (SR 50). Paired-end and single-end sequencing reads were generated using the Illumina HiSeq2500 instrument. FASTQ sequencing files were generated using the Illumina CASAVA software and internally qualityfiltered based on Illumina's chastityfiltering procedure, then PF reads (raw reads) were used for further analysis. Time series data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [29] and are accessible through GEO Series accession number GSE104807 (https://www.ncbi.nlm.nih.gov/geo/

query/acc.cgi?acc=GSE104807).

2.4. De-novo transcriptome assembly and functional annotation

RNA-seq raw reads were quality controlled using FASTQC v.0.11.3 and Trimmomatic [30] (automatically invoked by Trinity software [31]). The reference transcriptome was assembled with Trinity v.2.0.6 using the following parameters:–seqType fq–SS_lib_type FR–trimmomatic –quality_trimming_params SLIDINGWINDOW:20:25 MINLEN:80 and have been deposited at DDBJ/EMBL/GenBank under the accession GFXW00000000.

Trinotate v.2.0.2 (https://trinotate.github.io) was used for transcriptome annotation. Transcripts and their predicted ORFs were queried with BLASTx and BLASTp against both the SwissProt [32] and UniRef90 [33] databases and only the top-scoring hit was retained.

HMMER [34] was used to search for conserved protein domains in predicted ORFs against the pfam-A database [35]. BLAST homologs and Pfam domain entries were loaded into the pre-formatted Trinotate SQLite database which contained UniProt-associated annotation information. The subcellular localizations of each protein were predicted using a standalone version of TargetP v1.1 [36] in plant mode and without cut-oﬀs (winner takes it all).

2.5. Transcript abundance and diﬀerential expression analysis

Scripts bundled with Trinity software were mainly adopted to quantify the transcripts and ﬁnd the diﬀerentially expressed subset.

Transcript-level rather than unigene-level expression was investigated to retain isoform-speciﬁc information. Reads from individual libraries were aligned to the reference transcriptome with bowtie v.1.1.2 [37]

and quantified using RSEM v.1.2.19 [38]. A table of TMM-normalized FPKM expression matrix and a separate table of raw fragment counts were generated for further analysis and visualization. Differentially expressed (DE) transcripts were identified from raw counts with the Bioconductor package EdgeR v.3.1 [39]. Three biological replicates for each condition were provided. The most significant differentially expressed transcripts (FDR < 0.001 and FC > 4) were extracted for further analysis. A hierarchically clustered heatmap was generated from the Pearson correlation matrix of pairwise sample comparisons based on the most significant DE subset. The most significant DE subset was also partitioned into expression clusters by cutting hierarchical clustering trees of transcripts using MeV v.4.9.0 (https://sourceforge.

net/projects/mev-tm4/).

2.6. Gene ontology analysis

GO terms were extracted from transcriptome annotation using

scripts from Trinotate and visualized using BGI's WEGO [40]. GO functional enrichment tests on diﬀerentially expressed transcripts were performed with scripts from Trinity v2.4.0 that wrapped the Bio- conductor package GOseq [41]. Blast2GO v.4.0.7 [42] was used to map GO terms to plant GO Slim [43] terms in order to obtain a broad overview of the transcriptome. For each time point, enriched GO terms were visualized using GOplot v1.0.2 [44]. Most signiﬁcant over- and under-represented GO terms (top 10 for each time point, then combined) were extracted using custom scripts and visualized using the GOplot package (http://wencke.github.io/).

2.7. KEGG pathway analysis

A standalone version of KOBAS 2.1.0 [45] was used to functional annotate the transcriptome with KO terms, the corresponding pathways and enzymes were extracted from information retrieved using KEGG API. The TMM-normalized FPKM expression matrix was ﬁltered by removing low expressed transcripts (average expression value below 1) and transcripts without proper KO terms, then used as input for Path- view ([46]; online version) to visualize signiﬁcant pathways (internally use GAGE [47] for pathway enrichment analysis) and some manually selected pathways. Pathways with q-values smaller than 0.05 were combined and shown as heatmaps.

3. Results

E. oleoabundanswas grown in continuous mode in controlledﬂat- panel photobioreactors and subjected to full nitrogen-starvation during a 4-day period. Before and during this starvation period, culture parameters were monitored and physiological and transcriptional changes in culture samples measured to determine the temporal dynamics of the nitrogen-stress response. For the pre-cultures grown in nitrogen-replete conditions (+N) samples were taken at four timepoints in the light- phase preceding nitrogen starvation, while from the onset of nitrogen starvation (−N), time-series samples were taken at the start, middle and end of the light-phase and at the end of the dark-phase over a 4-day period (Fig. 1).

3.1. Biomass dynamics under nitrogen starvation

Upon nitrogen starvation,E. oleoabundanscould sustain an increase in biomass concentration (cell dry weight) during the light phase in the ﬁrst two days (Fig. 2A). Overall, the biomass concentration gradually decreased, which coincided with a decreased biomass productivity over time (Fig. 2D). The protein content and quantum yield (QY) also showed a steady decline throughout the starvation period (Fig. 2D and Supplementary Fig. A.1).

Nitrogen starvation induced the fast accumulation of carbohydrates, increasing from 15% at the start to a maximum content of 33% of the cell dry weight after 10 h starvation. Within the same time frame and with the same accumulation pattern, starch content increased from 8%

at the start to a maximum content of 28% of the cell dry weight (Fig. 2B). This indicated that the increase in total carbohydrates was mainly due to starch accumulation. The starch content subsequently decreased again to N-replete levels of 18% in the light phase at day 4. In contrast, accumulation of total fatty acids (TFA) was more gradual from a content of 16% at the start to a maximum of 41% after 4 days. The increase in lipid content was mainly due to TAG accumulation, which increased from a content of 0.4% to 30% after 4 days, while the polar lipid (PL) content was stable around 10–13% (Fig. 2C). This indicates that TAG biosynthesis occurred mainly de novo,although some interconversion from starch to TAG or from membrane lipids to TAG cannot be excluded. During the dark phases starch was consumed and decreased to 15% content, while a decrease of fatty acids (indicating fatty acid consumption), was not observed. In the hybrid producer E.

oleoabundansstarch therefore is the primary storage compound that is

(5)

accumulated in the light phase and that serves as a transitory energy source that is respired in the dark phase.

Concomitant with the accumulation of storage compounds the fatty acid composition of the diﬀerent lipid fractions also changed during nitrogen starvation. A clear change from polyunsaturated fatty acids (PUFAs) towards mono-unsaturated and saturated fatty acids was observed in all fractions (Fig. 3). Under N-replete conditions a large part of all lipid fractions was composed of C18:3 and C16:3 fatty acids (approximately 30% and 15%, respectively). These species decreased during N-starvation (to 5% and 2.5%, respectively), while there was an increase in the (mono-)unsaturated fatty acids C18:1, C16:2 and C16:1 and the saturated fatty acid C18:0 respectively. The biggest change was observed for C18:1 that increased from 5% to almost 40% in the lipid fraction. PUFAs with higher carbon chain lengths (e.g. C20:3 and C20:5) were not detected in any of the lipid fractions in N-replete and N-starvation conditions.

Production rates for total carbohydrates (TCH) and starch showed an identical rapid increase during the ﬁrst day of N-starvation. This indicates that increased TCH levels in the ﬁrst day can mainly be

ascribed to starch production. Subsequently, both TCH and starch production rates decreased to levels below those for N-replete growth, although more severe for TCH. This means that production rates of carbohydrates other than starch were reduced much stronger. In the lipid fraction TAG productivity increased slowly in the ﬁrst 2 days following N-starvation, and then decreased again to N-replete rates. The total lipid productivity remained stable on theﬁrst day of N-starvation and then started to decrease (Fig. 2D).

3.2. Diﬀerential gene expression dynamics under nitrogen starvation A total of 25.6 million high-quality RNA-seq clean reads were de novo assembled into 61,010 transcripts belonging to 42,338“genes”

with a N50 length of 1838 bp and average GC content of 58.9%

(Table 1). Functional annotation of the transcriptome was carried out with the Trinotate pipeline, and 44,789 transcripts (73.4%) were annotated using UniRef90 database (Table 1). The high-quality reference transcriptome served as a solid foundation for subsequent gene expression analysis.

Fig. 2.Dynamics of cellular content inE. oleoabundansgrown in N-replete (+N) and N-starvation (−N) conditions in air-liftﬂat-panel photobioreactor. (A) Biomass concentration [CDW g/L] and protein content [g/gCDW]. (B) Total carbohydrates (TCH) and starch content (% of CDW). (C) Total fatty acids (TFA), triacylglycerides (TAG) and polar lipids (PL) content (% of CDW). (D) Productivities of biomass (rCDW), total carbohydrates (rTCH), starch (rStarch), total fatty acids (rTFA) and triacylglycerides (rTAG). Results are shown as mean values ± standard deviations (n= 3). The onset of N-starvation is att= 0 h. The 12:12 light:dark cycles used for cultivation are shown as alternating 12 h light phases (white background) and 12 h dark phases (grey background).

(6)

Fig. 3.Fatty acid composition of (A) total fatty acids, (B) TAG and (C) polar lipids fractions ofE. oleoabundansgrown in N-replete (N+) and N-starvation (N−) conditions in air-liftﬂat-panel photobioreactors. Mean values (n= 3) are shown. Onset of N-starvation was att= 0 h.

(7)

3.2.1. Global gene expression analysis

Time series experiments were performed (7 time-points, 3 replicates; GEO accession GSE104807), reads from each RNA-seq library were individually aligned back to the reference transcriptome and the abundance of each transcript was determined. The sequencing depth for individual samples ranged from 8.5 to 19 million reads, and over 97%

of the reads could be mapped to the reference transcriptome for all the samples. The transcriptome coverage, defined here as the percentage of genes with TPM (Transcripts Per Kilobase Million) abundance bigger than one, was 68.2–80.5% (Supplementary Table 5). The hierarchical clustering of transcripts from samples obtained from triplicate photobioreactors indicated a good experimental reproducibility (Supplementary Fig. A.2). Differentially expressed transcripts were then determined from raw count values for each time point from the onset of N-starvation. A total of 12,975 transcripts were either up-regulated or down-regulated compared to control samples (0 h) and the number varied between different time points (2014-6589) (Supplementary Table 1).

Hierarchical clustering of expression profiles of each transcript re- sulted in 23 distinct expression clusters (Fig. 4and Supplementary Table 1). The largest clusters were cluster 21, 15, 11 and 4, which contained 2938, 1834, 1145 and 1081 differentially expressed gene isoforms respectively. Many transcripts related to fatty acid, TAG and starch metabolism belong to these clusters. The diverse expression clusters to which the genes from a metabolic route belong to, indicates that different genes from these metabolic pathways are most probably under the control of different transcriptional regulators and therefore do not follow simultaneous activation. Some transcriptional regulators group into these clusters as well and might play a regulatory role in the nitrogen-starvation activated response of fatty acid and starch biosynthesis genes. These include transcriptional regulators such as the 5′- AMP-activated and SNF1-related protein kinases (TR3188|c0_g1 and_g2, and TR24203|c0_g1), a dual specificity tyrosine-phosphorylation-regulated kinase (DYRK3: TR935|c0_g2_i1), PHR1-LIKE 1 protein (PHL1: TR19823|c0_g1) and a nitrogen regulatory protein P-II homolog (GLNB: TR228|c0_g1 and TR228|c0_g3). Clusters 1–3, 12, 19 and 20 showed very distinct expression peaks and interestingly several transcripts related to reactive oxygen stress grouped to the clusters 12, 19 and 20. Detailed description of the transcriptional response to nitrogen starvation is discussed below.

The GO terms for upregulated transcripts were enriched with terms related to fatty acid and lipid metabolism as well as carbohydrate and starch metabolism, as is expected in view of the observed biomass changes (Supplementary Fig. A.4 and Supplementary Table 2). In addition, GO terms for ammonium uptake and incorporation were enriched (direct response to absence of nitrogen) as well as GO terms related to reactive oxygen species (ROS) stress, indicating that N-starvation might cause the production of ROS. The GO terms for down- regulated transcripts were enriched with terms related to the

chloroplast and photosynthesis, which makes sense in view of the breakdown of pigments (cell bleaching) and reduced photosynthetic capacity (e.g. quantum yield) due to the lack of nitrogen (Supplementary Fig. A.4 and Supplementary Table 3). Taken together, GO analysis indicates that most of the transcripts can be grouped into a small number of broad yet distinct functional categories that correlate well to the observed biomass and physiological changes. GOBubble plots based on plant GO Slim terms (Supplementary Fig. A.4) showed that the highest number of increased and decreased GO Slim terms was found after 6 to 10 h of nitrogen starvation. This indicates that most transcriptional changes upon nitrogen starvation occur relatively rapid.

The broader GO Slim terms indicated that “carbohydrate metabolic process”(GO:0005975) was the only GO Slim term that for all time points signiﬁcantly increased compared toT= 0 (negative z-score H0 versus other time points). Other GO Slim terms like“catabolic process”

(GO:0009056) and “response to stress” (GO:0006950) also increased for an extended period. In contrast, the GO Slim terms for “photosynthesis” (GO:0015979), “plastid” (GO:0009536) and “thylakoid”

(GO:0009579) were clearly decreased from 6 to 29 h (positive z-score H0 versus other time points), indicating a reduction in the photosynthetic capacity.

KEGG pathway enrichment was performed using GAGE (Supplementary Table 4) and data were mapped on KEGG pathways to generate graphs showing detailed temporal changes (Supplementary Fig. A.5). The most striking observations from the KEGG pathway enrichment analysis was an up-regulation of the TCA cycle (ko00020), Proteasome (ko03050), Protein processing in endoplasmic reticulum (ko04141) and Phagosome (ko04145) pathways, and a down-regulation of the Photosynthesis (ko00195) pathway (Supplementary Table 4). This indicates increased protein and organelle degradation and correlates with the increased GO Slim terms such as“catabolic process”(GO:0009056),“generation of precursor metabolites and energy”(GO:0006091) and“protein metabolic process”(GO:0019538), as well as the decreased photosynthesis, plastid and thylakoid GO Slim terms mentioned above (Supplementary Fig. A.4) and decrease in quantum yield(Supplementary Fig. A.1).

3.2.2. Transcriptional changes in starch metabolism

The transient and rapid increase in starch productivity and accumulation in thefirst 10 h of N-starvation coincided with increased expression of several genes involved in starch biosynthesis (Fig. 5and Supplementary Table 1). In particular, starch biosynthesis genes encoding phosphoglucomutase (conversion of glucose-6-phosphate to glucose-1-phosphate) showed a strong up-regulation in expression from the onset of N-starvation, while the ADP-glucose phosphorylase (conversion of glucose-1-phosphate to ADP-glucose) was only up-regulated from 10 h on. Multiple starch synthase genes showed different responses. While the soluble starch synthase SSY3 and the granule-bound starch synthases SSG1 and SSG2 showed a fast and transient up-regulation from 2 to 6 h from the start, the soluble starch synthase SSY2 showed a strong transient increase in expression from 6 to 10 h. Like- wise, the 1,4-alpha-glucan-branching enzyme (GLGB1) that modifies amylose to starch was only up-regulated in thefirst 6 h. Among the starch degradation genes, gene expression of the alpha-glucan water dikinase (GWD1) was increased from 6 h on. This enzyme increases the phosphorylation-level of starch and thereby stimulates the degradation of starch in higher plants andChlamydomonas[48]. Starch breakdown involves several amylases and we found that multiple genes for α- amylase (AMY1 and AMY3),β-amylase (AMYB, AM1, BAM1) and isoamylase (ISOA1 and ISOA3) showedfluctuating expression responses.

Overall though, they exhibited down-regulation at the start of starch accumulation and up-regulation in the later stages, from 6 to 10 h on.

Hexokinase performs the ﬁrst step in glycolysis, the conversion of glucose to glucose-6-phosphate, and hexokinase genes were upregulated from 10 h starvation on. This indicates breakdown of glucose, which for a large part could be starch-derived glucose, as starch Table 1

Summary of de novo transcriptome assembly and annotation for E. oleoabundans.

Features Results

Reads Raw reads 27.19 M

Clean reads 25.58 M

Assembly Number of transcripts 61,010

Transcriptome size (Mbp) 70.8

Mean (bp) 1161

Median (bp) 800

N50 length (bp) 1838

GC (%) 58.9

Annotation Swiss-Prot 27,765 (45.5%)

UniRef90 44,789 (73.4%)

GO 21,591 (35.4%)

KO 15,336 (25.1%)

(8)

degradation increases from this time on as well. Overall, these results suggest that the increased production and accumulation of starch can mainly be attributed to the strong up-regulation of the starch biosynthesis genes and less to a decrease in starch degradation. The subsequent decrease in starch productivity and accumulation from 29 h on concurred with decreased gene expression for the soluble starch synthase SSY2 as well, which supports this hypothesis.

3.2.3. Transcriptional changes in fatty acid and lipid metabolism Transcriptional changes for genes related to lipid biosynthesis

(Fig. 6and Supplementary Table 1) also correlated with the observed increase in productivity and accumulation of TFA and TAG (Fig. 2).

Several genes required for fatty acid biosynthesis showed increased expression within 2 to 10 h from N-starvation. Thefirst step in fatty acid biosynthesis requires the precursor acetyl-CoA, which can be derived from different biochemical conversions. Afirst route is via the pyruvate dehydrogenase complex (PDH) that converts pyruvate to acetyl-CoA. A second route is via acetyl-CoA synthetase (ACS), that converts acetate to acetyl-CoA. Finally, a third route is via the TCA cycle where ATP citrate lyase (ACLY) can convert citrate to acetyl-CoA and oxaloacetate.

Fig. 4.Hierarchical clusters of gene expression proﬁles during nitrogen starvation. Expression levels are shown as median-centered log2(fpkm + 1) values for the time series H0-H99 and the number of transcript isoforms in each cluster are indicated.

(9)

Gene expression of the PDH subunits ODPA1 and ODPB1 (TR1574|c0_g2 and TR13897|c0_g2), several acetyl-CoA synthetases (TR16088|c0_g1/g2, TR18695|c0_g1, TR2994|c0_g1/g2/g3 and TR26776|c0_g1) and ATP citrate lyases (TR25574|c0_g2 and TR26430|c0_g2) clearly increased within 2 to 10 h from N-starvation and remained so until 99 h of starvation. Increased expression of PDH, ACS and ACLY genes can result in a higher supply of acetyl-CoA, the precursor for fatty acid biosynthesis, and thereby increase the carbon- ﬂux towards lipid synthesis (Fig. 6). The subsequent steps in fatty acid biosynthesis involve acetyl-CoA carboxylase (ACCase) and malonyl- CoA-ACP transacylase (FABD), and gene expression for these enzymes also increased in the ﬁrst 10 h from N-starvation. The fatty acid synthesis (FAS) complex genes were up-regulated as well, including the ones encoding chloroplastic 3-oxoacyl-ACP synthase I (KASC1), 3-oxoacyl-ACP reductase (FABG), NADPH-dependent enoyl-ACP reductase (FABL) and 3-hydroxyacyl-ACP dehydratase (FABZ), as well as the acyl- ACP thioesterase (FATA). Overall however, this response was more delayed than that for starch biosynthesis.

For subsequent biosynthesis of triacylglycerides (TAG), transcripts encoding glycerol kinase (GK), glycerol-3-phosphate dehydrogenase (GPDA), two glycerol-3-phosphate O-acyltransferases (GPAT) and the phosphatidic acid phosphatase (PAP) were strongly up-regulated from 2 h after the onset of N-starvation, while gene expression of lysophosphatidic acid acyltransferase (LPAT) - which is responsible for the intermediate conversion of lysophosphatidic acid to phosphatidic acid - only increased from 29 h on. These genes perform biosynthesis of the lipid backbone glycerol-3-phosphate (G3P) and the subsequent incorporation of free fatty acids up to diacylglycerol (DAG). The diacylglycerol O-acyltransferases (DGAT) genes required for the ﬁnal conversion of diacylglycerol (DAG) to triacylglycerol (TAG) also showed an initial decrease in expression, but were clearly up-regulated only from 29 h of N-starvation on, indicating that the DAG pool is built up before the ﬁnal DAG to TAG conversion. Both a DGAT-1 and a

DGAT-2 gene were up-regulated and appear to be responsible for the final conversion of DAG to TAG. The transcriptome annotation predicted a single DGAT-1 gene andfive DGAT-2 gene sequences, and our data indicate that only one of thefive DGAT-2 genes responds to N- starvation. The subcellular locations of these DGAT enzymes and the corresponding TAG biosynthesis pathways inE. oleoabundansare still to be determined. Interestingly, a gene encoding a putative lipid body- associated caleosin (peroxygenase) protein (TR11980|c0_g1_i1) was also up-regulated from 29 h on, indicating that increased lipid body formation is concomitant with DAG to TAG biosynthesis. Incorporation of free long-chain fatty acids into lipids requires the activation to their acyl-CoA form, which is usually performed by long-chain acyl-CoA synthetases (LACS). A gene isoform encoding a chloroplastic long- chain-fatty-acid CoA ligase (TR19322|c0_g1_i1) showed up-regulation from 53 h on and might serve this role.

TAG lipase activity negatively affects TAG accumulation, while phospholipases are involved in liberation of free fatty acids from (plastid) membrane lipids and interconversion into TAG. Transcript levels of two TAG lipase genes (SDP1: TR9756|c0_g3 and TR22634|c0_g2_i1) were decreased in the first 10 h of N-starvation, while at the end of N-starvation a marked increase was observed (Supplementary Table 1). The late upregulation of TAG lipases follows the earlier up-regulation of TAG biosynthesis genes, which could ex- plain the initial steady increase in TAG accumulation andfinal level- ling-off in TAG productivity. Besides de novo TAG synthesis, interconversion from membrane lipids to TAG is a common mechanism under stress conditions in microalgae [49–52], and the enzyme phos- pholipid:diacylglycerol acyltransferase (PDAT) plays an important role in this [51,53]. A PDAT homolog was not identified in our reference transcriptome, but several chloroplastic phospholipase A1 genes (DAD1 homologs) did show increased expression at different timepoints during N-starvation (Supplementary Table 1). Phospholipase activity results in the conversion of phospholipids to lysophospholipids and free fatty Fig. 5.Differential expression of starch biosynthesis genes during nitrogen starvation. Heatmaps of expression levels per transcript are shown as median-centered log2(fpkm + 1) values sequentially for time series T0-T99. The inset shows the log2 colour coding scale. PGM, phosphoglucomutase; AGPase, ADP-glucose phosphorylase; SS, starch synthases; HXK, hexokinase; PHS1, alpha-glucan phosphorylase; GLGB, 1,4-alpha-glucan-branching enzyme; ISA, isoamylase;α-AMY, alpha- amylase;β-AMY, beta-amylase; GWD, alpha-glucan water dikinase;α-1,4-Glc, alpha-1,4-glucan. Starch^℗indicates starch targeted for degradation through increased phosphorylation by GWD.

(10)

acids. Lysophospholipase subsequently converts lysophospholipids further to free fatty acids and glycerol-phosphate derivates [54,55]. A lysophospholipase gene (TR21369|c0_g1_i5) showed increased expression from 29 h on. Reacylation of lysophospholipids to phospholipids is performed by lysophospholipid transferases, and two lysophospholipid transferase genes (LPT1: TR19784 and TR9894|c0_g2_i1) were down- regulated within theﬁrst 10 h. The gene expression proﬁles of the (lyso) phospholipases and lysophospholipid transferases could indicate that free fatty acids are liberated from membrane phospholipids upon nitrogen starvation, which subsequently can be used in TAG biosynthesis.

The decrease in chloroplast-related GO Slim terms in the period before the clear increase in DGAT expression at 29 h also supports that TAG biosynthesis initially was not only de novo, but might also be the result of some interconversion from membrane lipids to TAG. Altogether, there are several indications that in the earlier stages of N-starvation there is active recycling of plastidial membrane phospholipids taking place, therefore incorporation of membrane-derived fatty acids into TAG cannot be excluded and awaits further investigation, e.g. by using carbon radiolabelling studies. Expression of a chloroplastic monogalactosyldiacylglycerol synthase gene (MGDGS: TR23849|c0_g1_i) was down-regulated from 6 to 53 h (Supplementary Table 1). MGDG synthase couples UDP-galactose and DAG to form the galactoglycerolipid monogalactosyldiacylglycerol (MGDG) which is abundantly present in the photosynthetic membranes of chloroplasts in higher plants and microalgae [56–58]. It thereby competes with DGAT enzymes for the chloroplastic DAG pool and decreased expression of MGDG synthase indicates reduced production of chloroplast-MGDG, leaving more DAG available for TAG synthesis. Finally, gene expression of chloroplast lipid transporters (TRIGALACTOSYLDIACYLGLYCEROL TGD1 and 2:

TR2378|c0_g1_i1 and TR29717) was up-regulated. TGD proteins are

present in the chloroplast membrane inChlamydomonas, and facilitate the transport of lipids derived from the endoplasmic reticulum (ER) into the chloroplast [58]. Up-regulation of these genes under nitrogen starvation inE. oleoabundanscould indicate that in this way starved cells try to compensate for the decrease in chloroplast lipids as we hypothesize above.

Expression of several chloroplast-specific fatty acid desaturase (FAD) genes was repressed from 6 h on (Fig. 7), which explains the increased incorporation of mono-unsaturated (C18:1 and C16:1) and saturated (C18:0) fatty acids in the different lipid fractions and decrease of C18:3 and C16:3 fatty acids. The combined reduction in expression of FAD6C (C16:1➔C16:2 and C18:1➔C18:2) and FAD3C (C16:2➔C16:3 and C18:2➔C18:3) supports the increase in the C16:1, C16:2 and C18:1 fatty acid fractions. Concomitant with this, expression of the elongase-complex genes 3-ketoacyl-CoA synthase (KSC1), very-long- chain (3R)-3-hydroxyacyl-CoA dehydratase (HACD) and very-long- chain enoyl-CoA reductase (ECR) showed increased expression. This could result in higher C16:0 to C18:0 conversion, and higher C18:0 levels could thereby increase C18:1 production as well. Both increased C18:0 and C18:1 levels were observed, though C16:0 levels were not clearly reduced. This indicates that most probably the observed increase in C18:1 and C18:0 is mainly caused by the reduction in fatty acid desaturases gene expression. In support of this we did not observe significant differential gene expression of two stearoyl-ACP desaturase genes (SAD3 and SAD5) in our reference transcriptome, which are responsible for desaturation of C16:0 and C18:0. From the start of N- starvation a chloroplastic FAD4 desaturase gene expression was also reduced (TR22593|c0_g2_i1), which converts [sn-2-C16:0]-containing glycerophospholipids into [sn2-C16:1]-containing glycerophospholipids, while at the same time expression of a chloroplastic Fig. 6.Differential expression of fatty acid and triacylglyceride biosynthesis genes during nitrogen starvation. Heatmaps of expression levels per transcript are shown as median-centered log2(fpkm + 1) values for the time series T0-T99. The inset shows the log2 colour coding scale. PDH, pyruvate dehydrogenase complex; ACS, acetyl-CoA synthase; ACCase, acetyl-CoA carboxylase; FABD, malonyl-CoA-ACP transacylase; KAS, 3-oxoacyl-ACP synthase; FABG, 3-oxoacyl-ACP reductase; FabZ, 3- hydroxyacyl-ACP dehydratase; FABL, NADPH-dependent enoyl-ACP reductase; FATA, acyl-ACP thioesterase; LACS, long-chain acyl-CoA synthetase, G3PDH, glycerol-3-phosphate dehydrogenase; GK, glycerol kinase; GPAT, glycerol-3-phosphate O-acyltransferases; LPAT, lyso-phosphatidic acid acyltransferase; PAP, phosphatidic acid phosphatase; DGAT, diacylglycerol O-acyltransferase; TAG lipase: triacylglyceride lipase.

(11)

palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase (ADS3, TR16760|c0_g1_i1) was upregulated in the early phases of N-starvation.

This indicates decreased C16:0 to C16:1 conversion in the glycer- ophospholipid fraction, while in the glycolipid fraction (speciﬁcally MGDG) there was an initial increased conversion of C16:0 to C16:1.

Overall, the combined eﬀect of changes in expression of several desaturase genes involved in C16:0 to C16:1 conversion, led to an increased C16:1 content in the lipid fractions. Finally, expression of a palmitoyl- protein thioesterase gene (TR646|c0_g3_i1) was upregulated in the late stages of N-starvation. This could supply additional free palmitate (C18:0) derived from lipid-containing proteins for incorporation into e.g. TAG.

3.2.4. Transcriptional changes in photosynthesis, carbon ﬁxation and oxidative stress

Expression of genes for the photosystems proteins initially increased (chlorophylla/bbinding proteins), which coincided with a decreased expression of chlorophyllase (Supplementary Table 1). This could indicate that as N-starved cells cannot synthesize novel (N-containing) photosystems due to the lack of nitrogen, cells try to maintain the level of active photosystems by reducing chlorophyll degradation and maintaining its incorporation chlorophyll into existing photosystem proteins. However, at 10 to 29 h post-starvation this can no longer be sustained, indicating a decrease in photosystem biosynthesis and (most probably) photosynthetic eﬃciency (as can also be observed from the decreased quantum yield). The observed decrease of chloroplast-related GO Slim terms in this period is in concordance with this.

Many genes involved in CO2fixation (RuBisCo, PGK, PEP carboxylase, PEPC kinase) and CO2concentrating mechanisms (carbonic an- hydrases) show a rapid up-regulation in expression in thefirst 2 h, a decrease afterwards and increase again after 99 h (Supplementary Table 1). This correlates to the initial sustained expression of photosystem genes and indicates that N-starved cells initially try to maintain photosynthesis and carbon fixation. Ultimately, when because of N- starvation cells cannot maintain active photosystems they can experi- ence oxidative stress by the production of reactive oxygen species (ROS) such as hydrogen peroxide (H2O2) and superoxide (O2%⁻) that can cause cell damage, as has been observed in e.g. Chlorella and Acutodesmus species [59–62]. In response to N-starvation, E. oleoabundans showed increased transcript abundance for superoxide

dismutase (SOD), catalase, peroxiredoxin and several peroxidases which indicates that these ROS were produced in this species as well (Supplementary Table 1).

3.2.5. Transcriptional changes in nitrogen uptake and incorporation As expected, N-starvation induced changes in expression of genes involved in nitrogen metabolism. The depletion of nitrate led to reduced gene expression for nitrate transporters, nitrate reductase and ferredoxin-nitrite reductase. At the same time, increased expression started for ammonium transporter genes and the glutamine synthetase (GS) and NADH-dependent glutamate synthase (GOGAT) pathways.

This indicates that cells switch to an alternative mode of nitrogen incorporation via ammonium uptake (Supplementary Table 1).

3.3. Transcriptional regulators of starch and lipid accumulation

The patterns in starch and lipid accumulation suggest the existence of a regulatory mechanism for redirecting carbon-partitioning (switching) between the associated metabolic pathways. Identiﬁcation of the associated“switch”regulators would be of great interest to improve existing metabolic engineering strategies that currently are often based on increased expression of a speciﬁc metabolic gene or pathway, or introduction of heterologous transcription regulators from higher plants [6].

In this study, we identiﬁed several transcriptional regulators with signiﬁcantly altered gene expression upon N-starvation, some of which have been implicated in the nitrogen-starvation triggered lipid accumulation in microalgae (Supplementary Table 1). Genes encoding 5’- AMP-activated and SNF1-related protein kinases (TR3188|c0_g1 and g2, and TR24203|c0_g1, respectively) clearly showed increased expression from 10 h N-starvation on, and very strongly from 53 h on. This is in line with a previous transcriptome study inE. oleoabundans[17]

where after 11d of N-starvation the 5’-AMP-activated protein kinase SnRK1 showed up-regulation. SnRK1 is a global regulator of carbon metabolism in plants [63,64].

InNannochloropsis,bZIP-type regulators have been implied to control several steps in the synthesis of TAG [22]. We identified two such bZIP regulators with opposite expression profiles in our transcriptome set: TR2010|c0_g1 and TR14421|c0_g1. TR14421|c0_g1 expression was activated 2 h after N-starvation, but showed down-regulation from 29 h Fig. 7.Differential expression of fatty acid desaturase genes during nitrogen starvation. Heatmaps of expression levels per transcript are shown as median- centered log2(fpkm + 1) values for the time series T0-T99. The inset shows the log2 colour coding scale. KSC1, 3-ketoacyl-CoA synthase; KCR1, (very long-chain) 3-ketoacyl-CoA reductase; HACD, very- long-chain (3R)-3-hydroxyacyl-CoA dehydratase;

ECR, very-long-chain enoyl-CoA reductase; SAD, stearoyl-ACP desaturase; FAD, fatty acid desaturase;

ADS3, palmitoyl-monogalactosyldiacylglycerol

delta-7 desaturase; MGDG, mono-

galactosyldiacylglycerol.^#No signiﬁcant diﬀerential gene expression was observed.

(12)

on, while TR2010|c0_g1 was down-regulated until 53 h, after which the gene was activated. We also identified a transcript (TR935|c0_g2_i1) with homology to the dual specificity tyrosine-phosphorylation-regulated kinase (DYRKP) from Chlamydomonas, a negative regulator of carbon storage that inhibits both starch and lipid accumulation and photosynthetic efficiency under N-starvation in this algae [65]. Ex- pression of this transcript was markedly decreased from the start but increased from 29 h starvation on. Differentially expressed transcripts encoding two homologs of PHR1-like regulators from higher plants were also found (TR19823|c0_g1 and TR25002|c0_g1). PHR1 and its homolog inChlamydomonasPSR1 are Phosphate Starvation Response transcription factors, and PSR1 was found to be essential for lipid and starch accumulation during phosphate starvation in Chlamydomonas [66]. A homolog of the nitrogen regulatory protein P-II (GLNB:

TR228|c0_g1 and g3) was upregulated from 6 h starvation on. The P-II regulator from Chlamydomonaswas found to negatively control TAG accumulation in lipid bodies during acclimation to nitrogen starvation in this green alga [67].

If and how these regulators aﬀect carbon partitioning for the accumulation of starch and lipids remains to be investigated through functional analysis.

4. Conclusions and discussion

Metabolic engineering of microalgae towards increased lipid or starch productivity requires deeper insight into the genetic details of the metabolic pathways involved and the regulatory mechanisms controlling these pathways. For the oleaginous green microalgaeE. oleoa- bundansthis genetic information is still very limited, therefore we set out to study the biomass and transcriptional changes that occur upon nitrogen starvation, tofind genetic triggers of starch and lipid accumulation in this industrially relevant microalga. Starch and lipid accumulation patterns showed different dynamics over the applied 4-day starvation period. Whereas the starch content increased strongly and reached its peak within thefirst day upon nitrogen starvation, lipids (mainly TAG) accumulated more gradually (Supplementary Fig. A.2). A pattern where starch accumulates rapidly ahead of lipid accumulation is observed in other hybrid oleaginous microalgae [9–11]. Starch seems to be the preferred primary storage compound in these species, which is partially due to the fact that its synthesis is biochemically and en- ergetically more favourable than fatty acid biosynthesis [10,68]. On the other hand, lipids have a higher energy density and therefore for the longer term might be the preferred storage compound. Several studies indicate that starch and lipid accumulation are both dependent on de novo biosynthesis and therefore competition for common carbon precursors affects carbon-partitioning between starch and lipids. In- creasing the carbon precursor supply can thereby be a strategy to in- fluence carbon-partitioning and improve lipid accumulation in some cases [69,70]. Nitrogen starvation also affected the fatty acid composition of the lipid fractions, with a distinctive decrease of the polyunsaturated fatty acids C18:3 and C16:3, concomitant with a large increase of the mono-unsaturated fatty acid C18:1, and small increases of the mono-unsaturated fatty acid C16:1 and the saturated fatty acids C18:0 and C16:0. These changes in fatty acid profile correspond to observations from earlier studies on N-starvation in E. oleoabundans [71,72], and was observed in other hybrid oleaginous microalgae as well [9,73]. Besides these metabolic factors, several biosynthetic genes and transcriptional factors have been identified that play a role in accumulation and carbon-partitioning of starch and lipids [6,20,22,53,65]. We therefore determined the transcriptional changes upon N-starvation through RNA-sequencing and observed that gene transcripts for starch- and lipid-related metabolic pathways correlated well with the observed biomass changes. The accumulation of starch in thefirst day could mainly be ascribed to a general increase in transcript abundance for several starch biosynthesis genes. The peak in starch accumulation correlated particularly well with the rapid increase and

subsequent decrease in expression of starch synthase genes (SSY2, SSY3, SGG1, SSG2), while the ultimate return to N-replete starch levels also appears to be caused by an increase in transcript levels of diﬀerent amylase genes (α-,β- and isoamylases) in the later stages.

Expression of fatty acid biosynthesis genes and TAG biosynthesis genes was in general more delayed compared to those for starch biosynthesis and remained at increased levels until 4 days of starvation.

Transcripts for the enzymatic reactions performed by PDH, ACS, ACLY and ACCase, which supply the fatty acid biosynthesis precursors acetyl- CoA and malonyl-CoA respectively, were clearly increased, but espe- cially PDH was constitutively upregulated from 6 to 99 h. Likewise, for TAG biosynthesis the genes for G3PDH and glycerol kinase (glycerol-3- phosphate synthesis) and GPAT (incorporation offirst acyl chain to glycerol-3-phosphate backbone) were in general up-regulated constitutively from 2 h starvation on, while other TAG biosynthesis genes had a more limited time-frame where they showed upregulation. In both pathways (fatty acid and TAG biosynthesis) there appears to be a trend where continuous supply of pathway precursors is important for the gradual accumulation of lipids/TAG. Metabolic engineering of microalgae for increased lipid production has focused on early and late genes in both pathways, either based on the hypothesis that the supply of fatty acid precursors is key for an increasedflux through the FA biosynthesis or on the hypothesis that the ultimate conversion(s) steps to TAG are rate-limiting. Effective and ineffective metabolic engineering for increased production of lipids/TAG has been reported for both thefirst (e.g. [6,74]) and the second approach (e.g., [6,75]), yet it still remains unclear whether the key metabolic nodes for increasing lipid/TAG production are conserved among the plethora of oleaginous microalgae. In this respect, a recent study has demonstrated that overexpression of a native DGAT-2 gene inE. oleoabundanscan more than double the lipid and TAG content and productivity, as well as alter the saturated fatty acid composition of cells [76]. This indicates that metabolic engineering focused on thefinal conversion of DAG to TAG is a successful strategy for increased lipid production inE. oleoabundans.

The DGAT-2 gene that was used for overexpression is not the one showing diﬀerential gene expression in our study (TR20062|c0_g1), but another one of the four DGAT-2 genes found in the reference transcriptome (TR12237|c0_g1, g2 and g3). On the other hand, the location of de novo TAG biosynthesis remains unclear at this moment. While most of the fatty acid biosynthesis transcripts were predicted to be lo- cated in the chloroplast, the subcellular location of the TAG biosynthesis transcripts was variable or unspeciﬁed (Supplementary Table 1).

Whether TAG biosynthesis occurs both in the chloroplast and on the endoplasmic reticulum inE. oleoabundans, as is the case in many microalgae [77], therefore needs to be veriﬁed by further biochemical experimentation.

Several transcriptional regulators have been implicated so far to play a central role in lipid accumulation in microalgae. These include the bZIP-type regulators that are suggested to control several steps of TAG synthesis in Nannochloropsis [22], and the nitrogen-responsive regulator (NRR1) and phosphorus stress response regulator (PSR1) in Chlamydomonas[20,53]. While these putative switch regulators were identified through transcriptomics and in silico studies, the direct role of thePSR1gene in lipid accumulation was confirmed by mutagenesis and overexpression. The kinase DYRKP has been implicated to repress both starch and lipid accumulation inChlamydomonasand to result in sustained photosynthetic efficiency under N-starvation [65]. From these described transcriptional regulators, we identified several homologs inE. oleoabundansthat showed significant differential gene expression upon nitrogen starvation. The actual role of these regulators in starch and lipid accumulation and carbon partitioning now needs to be further analysed through e.g. gene inactivation and gene overexpression approaches.

Finally, it remains unclear whether accumulation of these storage compounds inE. oleoabundansis controlled by genetic switches only (i.e. via speciﬁc activating or repressing regulators) or that