Identification of Targets of Transcription Factor WRINKLED1-Like Related to Lipid Biosynthesis From Marine Microalga Dunaliella parva

WRINKLED1 (WRI1) is an important transcription factor controlling lipid biosynthesis. To elucidate the function of Dunaliella parva WRI1-like (WRI1-like) (i. e., DpWRI1-like), the targets of DpWRI1-like were identified by chromatin immunoprecipitation sequencing. The results showed that DpWRI1-like regulated many target genes involved in carbohydrate metabolism, lipid metabolism, photosynthesis, and transcription factor. It was proposed that DpWRI1-like participated in a regulatory network controlling lipid biosynthesis. This work laid a good foundation for a deep understanding of the regulatory mechanism of DpWRI1-like in D. parva.


INTRODUCTION
The former studies indicated that nitrogen limitation could increase lipid content of microalgae (Yamaberi et al., 1998;Weldy and Huesemann, 2007;Wang et al., 2009;Moellering and Benning, 2010;Pancha et al., 2014). Therefore, nitrogen limitation was widely utilized to enhance the expression level and diversity of genes associated with the biosynthesis and catabolism of biofuel precursors. However, nitrogen limitation often resulted in the decrease of biomass, which confined its application. Our former work also showed that nitrogen limitation caused an increase of lipid content from 25% to 40% in D. parva. However, the biomass of D. parva had a decrease of 39% under nitrogen limitation condition compared with nitrogen sufficient condition (Shang et al., 2016). Up to now, studies on the response mechanism to nitrogen limitation was inadequate at the metabolic and molecular level in model microalgae such as Chlamydomonas reinhardtii, and a deeper exploration was required about the relevant metabolic pathways and gene regulatory networks in microalgae (Hu et al., 2008;Miller et al., 2010;Xiao et al., 2013).
Dunaliella parva is a unicellular halophilic green alga, which simultaneously accumulates lipids and large amounts of carotenoids (primarily beta-carotene), and shows a good resistance against unfavorable ambient conditions such as UV-A radiation, salinity up-shock, extreme pH and temperature, and nitrogen limitation (Mogedas et al., 2009;Borowitzka, 2013;Zhang et al., 2014;Srinivasan et al., 2015). In addition, D. parva has no rigid cell wall, which is favorable for genetic manipulation and extraction (Hosseini-Tafreshi and Shariati, 2009;Barzegari et al., 2010). These unique characteristics make D. parva more attractive as the feedstock for the production of biofuel and beta-carotene compared with model microalgae.
The application of genetic engineering method in lipid overproduction has recently emerged as a burgeoning trend in microalgae. Gene expression starts from the binding of various proteins to specific non-coding DNA sequences, such as motifs (Bylino et al., 2020). The regulation mechanism of gene expression is the focus of biological research. Identification of DNA regulatory elements, especially DNA motif, is of great significance for the study on regulation mechanism of gene expression.
Transcription factor (TF) can regulate the expression of multiple key enzyme genes related to the production of desired metabolites (Capell and Christou, 2004;Courchesne et al., 2009). AP2-type TFs are a family of TFs with AP2 domain. The former studies investigated different types of AP2-type TFs with various functions. The paralogous genes DORNRÖSCHEN and DORNRÖSCHEN-LIKE encoding AP2type TFs potentially promoted G1/S transition in Arabidopsis shoot meristem (Seeliger et al., 2016). An unusual AP2-type TF (SMALL ORGAN SIZE1) controlled organ size downstream of auxin signaling pathway (Aya et al., 2014). DREB2 class of AP2type TF in lily mediated heat-stress response (Wu et al., 2018). Four APB genes encoding AP2-type TFs determined stem cell identity in Physcomitrella patens (Aoyama et al., 2012).
WRINKLED1s are important AP2-type TFs, which are related to lipid accumulation (Ma et al., 2013). Researchers had conducted a series of in-depth studies about WRI1 in plants. WRI1 homologs and their target genes are largely expressed in the endosperm of oat and castor, but not in the starchy endosperm of wheat (Yang et al., 2019). GmWRI1a TF (one of three distinct orthologs of WRI1 in soybean) could positively regulate oil production in soybean seed through a complex gene expression network associated with fatty acid biosynthesis (Chen et al., 2018). WRI1 and WRI1-regulated genes related to fatty acid biosynthesis were upregulated, resulting in a corresponding increase in the rate of fatty acid biosynthesis (Adhikari et al., 2016). WRI1 positively regulated genes coding for fatty acid synthesis and BIOTIN ATTACHMENT DOMAIN-CONTAINING proteins, representing a coordinated mechanism to attain lipid homeostasis . A transcriptional regulation model supported a viewpoint that interaction of 14-3-3 protein with AtWRI1 enhanced the transcriptional activity through protein stabilization (Kong and Ma, 2018a). A negative feedback loop of WRI1 expression was an indirect interaction most possibly caused by downstream targets of WRI1 (Snell et al., 2019). The increasing levels of different WRI1 homologs (from Avena sativa L., Cyperus esculentus L., Arabidopsis thaliana and Solanum tuberosum L.) reduced the transcriptional activity of Arabidopsis WRI1 promotor region (Snell et al., 2019). WRI1 promoted the flow of carbon to lipid during seed maturation through directly activating genes associated with fatty acid biosynthesis and controlling genes related to assembly and storage of triacylglycerol in Arabidopsis (Maeo et al., 2009). The mRNA level of Prunus sibirica WRI1 was closely related to lipid accumulation in Siberian apricot kernel (Deng et al., 2019). There were a few reports about WRI1 in microalgae. AP2 type TF WRINKLED1 in Arabidopsis (AtWRI1) regulated a number of genes related to lipid biosynthesis in Nannochloropsis salina, resulting in the increase of neutral lipid and fatty acid methyl ester production (Kang et al., 2017). These results indicated that heterologous expression of AtWRI1 might be utilized for biodiesel production in industrial microalgae. In our former work, DpWRI1-like gene (GenBank KR185335) was identified, which was also associated with oil accumulation (Shang et al., 2016). In a series of recently published studies using plants, researchers indicated that WRI1, a well-known key TF related to the allocation of carbon into oil, had attracted much interest (Kong and Ma, 2018b;Kong et al., 2019). However, WRI1 gene was not well elaborated in microalgae (especially in non-model microalgae such as D. parva).
Previous works on the response of Dunaliella to different stresses have been carried out by transcriptome sequencing method, which provided some good ideas for the study on D. parva. The common differentially expressed genes (DEGs) were identified under hypersaline and normal conditions by Empirical Bayes method, and the novel hub genes were revealed in halophytic microalga Dunaliella salina to understand the underlying molecular mechanism (Panahi and Hejazi, 2021). The first full-scale transcriptome research of D. salina was performed under salt stress, and salt could regulate the expression of key genes (Hong et al., 2017). Transcriptome analysis of D. salina was performed under hypersaline condition (2.5 M NaCl) at different times (He et al., 2020). The heat stress response of Dunaliella bardawil was studied by transcriptome analysis (Liang et al., 2020). High salinity acclimatization could alleviate cadmium toxicity in D. salina by physiological and transcriptomic analysis .
Our former study identified the genes related to the biosynthesis and degradation of biofuel precursors (such as fatty acid, triacylglycerol, and starch) and DEGs, including the related TF gene DpWRI1-like under nitrogen limitation condition in D. parva (Shang et al., 2016). In addition, the promoter and fulllength cDNA of DpWRI1-like were cloned. Nitrogen limitation significantly induced the expression of DpWRI1-like and the increase in lipid content (from 25% to 40%) in D. parva (Shang et al., 2016). It was inferred that DpWRI1-like played an important role in the regulation of lipid biosynthesis. However, the target genes regulated by DpWRI1-like are still unclear in D. parva. It is of interest to study the target genes of DpWRI1-like at the molecular level, which will lay the foundation for a deep understanding of the molecular mechanisms of lipid biosynthesis in D. parva, and a great enhancement of lipid content through manipulating the expression of TF gene DpWRI1-like in microalgae.
Here, DpWRI1-like was overexpressed in E. coli Rosetta (DE3) strain, and the fused protein DpWRI1-like with His Tag was purified through Ni 2+ -NTA column. Then the purified DpWRI1-like was used to immunize rabbit to obtain rabbit anti-DpWRI1-like polyclonal antibody. After that, the purified polyclonal antibody was used for chromatin immunoprecipitation (ChIP). Then the enriched DNA fragments by TF DpWRI1-like were sequenced through ChIP sequencing (ChIP-Seq) technology. At last, the target genes regulated by DpWRI1-like were obtained and analyzed in detail.

Culture and cDNA Synthesis
Dunaliella parva FACHB-815 was obtained from freshwater algae culture collection at the Institute of Hydrobiology (Wuhan, China). Cells of D. parva were cultured in Dunaliella medium with two nitrogen concentrations (nitrogen limited condition, 0.5 mM NaNO 3 , treated sample D_parva_T, and nitrogen sufficient condition, 5 mM NaNO 3 , control sample D_parva_C, respectively) in triplicate for 11 days. The composition of Dunaliella medium and culturing conditions were described by our former study (Shang et al., 2017). The dry cell weight of D. parva was calculated by measuring the weight difference of dried filtration membrane after drying algal cells. Cells of D. parva were collected by centrifugation (Centrifuge 5702R, Eppendorf, Germany) at 15,294 × g for 10 min at 4 • C. Total RNA was extracted from D. parva cells using RNAiso Plus reagent (Takara Bio, China). The cDNA of D. parva was synthesized using PrimeScript TM 1st Strand cDNA Synthesis Kit (Takara Bio, China).

Construction of Recombinant Expression Vector of DpWRI1-Like
Based on our published DpWRI1-like sequence (Shang et al., 2016), a pair of primers 30a-wri1(CF) and 30a-wri1(CR) ( Table 1) were synthesized and used for cloning of DpWRI1like cDNA into pET-30a vector (NdeI and XhoI site) to form pET-30a-DpWRI1-like. The pET-30a-DpWRI1-like was used to transform E. coli Rosetta (DE3) strain by chemical transformation method (Li et al., 2017). All primers used in this study are shown in Table 1.

Expression and Purification of Recombinant DpWRI1-Like
The positive transformant of E. coli Rosetta (DE3) strain containing pET-30a-DpWRI1-like was inoculated into 40 mL LB medium containing 50 µg/mL kanamycin and incubated at 37 • C, 220 rpm until OD 600 value reached 0.6-0.8. Then 40 mL culture was inoculated into 4-L LB medium containing 50 µg/mL kanamycin. The 4 L culture grew at 37 • C, 220 rpm until OD 600 value reached approximately 0.6, and then was induced with isopropyl β-D-thiogalactopyranoside (IPTG) (final concentration of 0.5 mM) for 4 h.
Cell pellet was harvested after induction by centrifugation (Centrifuge 5702R, Eppendorf, Germany) at 15,294 × g for 10 min at 4 • C. Cell pellet was resuspended in lysis buffer consisting of 8 M urea, 20 mM phosphate buffered saline (PBS), 300 mM NaCl, 0.1% TritonX-100, 1 mM dithiothreitol (DTT), pH 8.0. Cell resuspension solution was disrupted by sonication for 20 min on ice (2 s ON and 6 s OFF) with 400 W power. After sonication, the supernatant was obtained by centrifugation at 15,294 × g for 20 min at 4 • C. Then the supernatant was further purified using Ni 2+ -NTA column (QIAGEN, Germany). The fusion protein DpWRI1-like from recombinant pET-30a vector contains a His-Tag, which is often used for affinity purification

Primer
Sequence The underlined sequences CATATG, CTCGAG, and AAGCTT indicated the recognition sites of the respective restriction enzymes NdeI, XhoI, and HindIII. and identification of the purified fusion protein via Western blot (Zhao et al., 2013). The purification of recombinant DpWRI1-like by Ni 2+ -NTA column was performed using the method described by the former study (Amin and Latif, 2017). The binding buffer consists of 8 M urea, 20 mM PBS, 300 mM NaCl, 0.1% TritonX-100, 1 mM DTT, pH 8.0. The wash buffer and elution buffer were prepared by adding 50 mM and 500 mM imidazole to binding buffer, respectively. The eluted components with better purity were dialyzed against solution consisting of 1xPBS, 0.1% N-lauroylsarcosine sodium salt, pH 8.0. Then the dialyzed DpWRI1-like was filtered by filter membrane (0.45 µm) to obtain the purified DpWRI1-like.

Sodium Dodecyl Sulfate-Polyacrylamide Gel Electrophoresis, Western Blot, and the Measurement of Concentration for the Purified DpWRI1-Like
The purified DpWRI1-like was monitored by SDS-PAGE. The concentrations of acrylamide were 12% and 5% in separation gel and stock gel, respectively. The operation procedure were as follows. The protein went through stock gel electrophoresis under 80 V for 20 min, then separation gel electrophoresis under 120 V for 60 min. The resulting gel was stained by Coomassie brilliant blue solution for 20 min and decolorized by 10% acetic acid.
Western blot assay of purified DpWRI1-like was performed to check the authenticity of recombinant DpWRI1-like by standard procedure (Halabian et al., 2009). The concentrations of separation gel and stock gel were 12% and 5%, respectively. The first antibody and the second antibody were anti-6 × His rabbit polyclonal antibody and HRP-conjugated goat anti-rabbit IgG (Sangon Biotech, China), respectively. The chromogenic reagent of Western blot assay was 3,3,5' ,5'-tetramethyl benzidine (TMB, Sangon Biotech, China). The protein concentration of purified DpWRI1-like was determined using non-interference protein assay kit (Sangon Biotech, China) and standard protein bovine serum albumin (BSA).

Preparation of Polyclonal Antibody of Purified DpWRI1-Like
The purified DpWRI1-like was used to immunize rabbit. Then the serum was collected and rabbit anti-DpWRI1-like polyclonal antibody was purified. The rabbit anti-DpWRI1-like polyclonal antibody was tested via Western blot for D. parva cells expressing natural DpWRI1-like protein. The prepared rabbit anti-DpWRI1-like polyclonal antibody might serve as a tool in ChIP assay to detect the binding fragments of TF DpWRI1-like. The details were as follows.
Polyclonal antibody in this study was obtained by immunizing rabbit with the purified DpWRI1-like as the immunogen. The experiment was performed using New Zealand white rabbit (female, 4 months old, 2.1 kg, Beijing HFK Bioscience, China). The DpWRI1-like solution was emulsified with Freund's complete adjuvant (FCA) and injected into the rabbit. FCA was used as an immune response enhancer in primary immunization of DpWRI1-like (Barinova et al., 2017). The DpWRI1-like solution was emulsified with Freund's incomplete adjuvant (FIA) in secondary immunization, the third immunization, and the fourth immunization. FCA and FIA were purchased from Sigma-Aldrich (St. Louis, MO, United States). After the fourth immunization, a small amount of blood was collected to provide a continuous IgG. The serum was extracted from the blood by centrifugation. The serum contained different types of immunoglobulins and several other proteins. Therefore, IgG from the serum was purified by affinity column chromatography to obtain pure IgG polyclonal antibody for ChIP.
ELISA for antibody titer determination was conducted. The antigen (purified recombinant DpWRI1-like) was dissolved in 0.05 M carbonate (pH 9.6) and filled into microtiter plate with 100 µL/well (0.2 µg/well) before overnight incubation at 4 • C. The microtiter plate was washed three times (3 min/time) with 250 µL/well PBST solution (PBS with 0.05% Tween-20). Then the blocking agent (5% non-fat dry milk in 0.01 M PBS, pH 7.4) was filled into the plate with 100 µL/well before incubation for 1 h at 37 • C. After that, the wells were washed three times (3 min/time) with 250 µL PBST/well. The purified antibody from the fourth immunization was filled into the microtiter plate (100 µL/well) as the first antibody and incubated for 1 h at 37 • C. After that, the wells were washed three times (3 min/time) with 250 µL PBST/well. The diluted HRP-labeled Goat Anti-Rabbit IgG (H + L) solution (1:8000 in 0.01 M PBS, pH 7.4) was added to the plate with 100 µL/well as the second antibody. Incubation was done for 45 min at 37 • C, followed by washing the wells five times (3 min/time) with 250 µL PBST/well. After that, TMB was added to each well with 100 µL/well and incubated for 15 min at room temperature. At last, 2 M sulfuric acid solution was added to the wells with 100 µL/well to terminate the reaction. Then the absorbance was determined at 450 nm using microplate spectrophotometer (ST-360, Shanghai Kehua Bio-Engineering Co., Ltd., China).
After obtaining polyclonal antibody of purified DpWRI1like, the specificity of the antibody was checked by Western blot. Proteins were extracted from treated sample D_parva_T and control sample D_parva_C using One Step Plant Active Protein Extraction Kit (Sangon Biotech, China). Western blot assay of D. parva proteins was performed by standard procedure (Halabian et al., 2009).

Chromatin Immunoprecipitation of DpWRI1-Like and Chromatin Immunoprecipitation Sequencing
Chromatin immunoprecipitation sequencing is a wellestablished method to identify the binding sites of DNA binding proteins such as TFs (Brdlik et al., 2014;Yao et al., 2017). The simple process of ChIP was as follows. DNA binding proteins of D. parva cells were crosslinked to DNA with formaldehyde (final concentration of 1%) for 10 min. Then D. parva cells were lysed with lysis buffer containing protease inhibitor, and DNA was sheared by micrococcal nuclease at 37 • C for 15 min. DNA fragments bound to a protein of interest, TF DpWRI1-like, were immunoprecipitated using the purified rabbit anti-DpWRI1-like polyclonal antibodies (5 µg) and protein A/G plus agarose. Crosslinking was then reversed through extensive heat incubation (65 • C for 1.5 h) and digestion of protein component with proteinase K, and the precipitated DNA fragments were isolated using cation-exchange resin Chelex 100. Normal rabbit IgG and anti-RNA polymerase II antibody were used as negative control and positive control to validate the reliability of ChIP assay, respectively. The purified DNA fragments were sequenced by ChIP-Seq after library construction. Sequenced reads were aligned with the database (C. reinhardtii genome) to obtain the corresponding genes. The negative control was used to discard and assess false-positive sequences in ChIP-Seq analysis. The details of library preparation were as follows.
The DNA library was prepared using paired-end sample preparation kit (Illumina, United States) according to Illumina's protocol. The ChIPed DNA fragments were blunt ended and phosphorylated, and a single ' A' nucleotide was added to the 3 ends of the fragments in preparation for ligation to an adapter that had a single-base 'T' overhang. Adapter ligation at both ends of genomic DNA fragment conferred different sequences at the 5 and 3 ends of each strand in genomic fragment. The products of this ligation reaction were purified and size-selected by agarose gel electrophoresis. Size-selected DNA was amplified to enrich for fragments that had adapters on both ends. The resulting sample library was again purified and size-selected by agarose gel electrophoresis. The final purified product was then quantitated prior to seeding clusters on a flow cell. At last, the library was sequenced by BGI Genomics using Illumina HiSeq 2000 platform (Illumina, United States). CHIP-Seq analysis with three replicates was performed.

Determination of Protein and Transcript Levels of DpWRI1-Like
DpWRI1-like transcript levels of D_parva_T and D_parva_C were measured in our previous study (Shang et al., 2016). DpWRI1-like protein levels of D_parva_T and D_parva_C were measured as follows. Firstly, proteins from D_parva_T and D_parva_C were extracted using One Step Plant Active Protein Extraction Kit (Sangon Biotech, Shanghai, China). Secondly, protein content was quantified by BCA Protein Assay Kit (Sangon Biotech, Shanghai, China) using BSA as standard protein. Three replicates were used for the measurement of protein levels.

Quantitative Real-Time PCR Analysis
The GAPDH promoter was tested using a pair of primers GAPDH-F and GAPDH-R (

Validation of Target Gene MME5
The interaction of DpWRI1-like and the promoter of target gene MME5 from D. parva was verified by Matchmaker Gold Yeast One-Hybrid Library Screening System kit (Takara, China) to prove TF activity of DpWRI1-like and ChIP analysis output. The upstream promoter sequence (1.5 kb) of MME5 gene encoding D. parva malic enzyme was amplified by primers ME-PF and ME-PR (Table 1) using D. parva genomic DNA as template. Then PCR product was purified, digested with HindIII and XhoI, and then cloned into HindIII and XhoI sites of pAbAi vector to form bait plasmid. The BbsI-linearized bait plasmid was integrated into the genome of Y1HGold strain on solid agar synthetic defined (SD) medium-Ura by PEG/LiAc method to form bait yeast.

Construction, Expression, and Purification of Recombinant DpWRI1-Like
The positive E. coli Rosetta (DE3) strain containing expression construct pET-30a-DpWRI1-like was sequenced. Blastn search revealed that 1,023 bp cDNA was same as our published DpWRI1-like sequence. The construction of recombinant vector is shown in Figure 1.
Then the induction of DpWRI1-like gene was achieved. Through the induction of expression, the fusion protein with His-Tag (i.e., recombinant DpWRI1-like) was accumulated. The detailed process of recombinant DpWRI1-like production is shown in Supplementary Figure 1.
The fusion protein DpWRI1-like contained a His-Tag which was often used for affinity purification and identification of purified fusion protein via Western blot (Debeljak et al., 2006;Hwang et al., 2014;Wingfield, 2015). After affinity purification by Ni 2+ -NTA column, SDS-PAGE was used to examine purified recombinant DpWRI1-like (Figure 2). The result showed that the molecular weight of purified DpWRI1-like was about 39 kDa, which was consistent with the predicted size.
To further validate the band as DpWRI1-like, Western blot was conducted based on His-Tag. The color reaction of TMB (Figure 3) verified the existence of a His-Tag protein.
Then the purified DpWRI1-like (1.72 mg/mL) was used to immunize rabbit.

Preparation of Polyclonal Antibody of Purified DpWRI1-Like
Freund's incomplete adjuvant can avoid inflammation and lesions resulting from the formation of granulomas (Talib et al., 2018).  After the fourth immunization, the serum was extracted and purified to obtain pure IgG polyclonal antibody.
The quantity of antibody identifying a specific antigen was evaluated by antibody titer (Delahaut, 2017). As shown in Table 2, the antibody titers of purified anti-DpWRI1-like polyclonal antibody were as follows: A ≥ 256 K, B ≥ 512 K, C ≥ 512 K, D ≥ 512 K. Therefore, antibody titers can meet the demand for further experiment (ChIP).
Then the specificity of antibody against D. parva proteins was tested. Anti-DpWRI1-like polyclonal antibody could specifically bind to native DpWRI1-like protein in treated sample D_parva_T, control sample D_parva_C, and recombinant DpWRI1-like protein. Therefore, the specificity of antibody can meet the demand for ChIP. The transcript and protein levels for DpWRI1-like under N-deficient condition were 19.85-fold and 12.5-fold of that under N-sufficient condition at 9 days based on our previous study (Shang et al., 2016) and the measurement of this study (protein levels under N-deficient and N-sufficient conditions were 52.5 µg/mL and 4.2 µg/mL), respectively. Chromatin immunoprecipitation was performed using anti-DpWRI1-like rabbit polyclonal antibody on sheared chromatin from D. parva cells using Pierce Agarose ChIP Kit (Thermo Fisher Scientific, United States). Fragment size of sheared DNA mainly was 100 bp-500 bp, which is shown in Supplementary  Figure 3. Normal rabbit IgG was used as a negative control. Anti-RNA polymerase II antibody was used as a positive control. Normal rabbit IgG could not bind to RNA polymerase II, therefore it did not pull down the promoter of GAPDH gene from D. parva. In contrast, anti-RNA polymerase II antibody could bind to RNA polymerase II, therefore it could pull down the promoter of GAPDH gene from D. parva. The results of the negative control and positive control (Figure 4) showed that ChIP experiment was effective.

Analysis and Annotation of Chromatin Immunoprecipitation Sequencing Data
The ChIP-Seq data of treated sample D_parva_T were assembled according to the resulting sequencing reads. Sequencing in sample D_parva_T generated a total of 19,733,651 clean reads after removing reads with adaptors and unknown nucleotides, overall low quality reads. The clean rate and clean reads Q20 rate were 99.93% and 95.87% for D_parva_T, which indicated that ChIP-Seq data quality was qualified. The clean reads were aligned with the genome of C. reinhardtii using SOAPaligner/SOAP2 (version 2.21t). The assembly was better using model organism C. reinhardtii as the reference genome from a data point of view (Beijing Genomics Institute, China). After alignment, there were 159,448 mapped reads and 44,048 uniquely mapped reads in sample D_parva_T.
Then peak calling was conducted by MACS (Model-based Analysis of ChIP-Seq) software to obtain significant regions (i.e., A to D indicated four tubes of antibodies. 1K indicated that antibody was diluted at a ratio of 1:1,000. Negative control used non-immunized rabbit serum. Blank group used PBS buffer. When antibody titer was more than 2.5-fold of negative control, the titer was considered as effective. peaks) which were likely to be the binding sites of TF DpWRI1like. There were 2,204 peaks with average length of 349 bp. The distribution of the peak size in sample D_parva_T is shown in Figure 5. The X-axis and Y-axis indicated length distribution and the number of ChIPed DNA fragments in sample D_parva_T, respectively. The sequencing and mapping depths in sample D-Parva_T are shown in Supplementary Figure 4. According to annotated information of UCSC database, genes were divided into different functional elements such as intergenic, introns, downstream, upstream, and exons. A total of 2,204 genes in total had DpWRI1-like binding peaks in sample D_parva_T. The distribution of peaks in different functional elements is shown in Figure 6. The result indicated that the distribution of peaks was primarily in introns (50.3%) and intergenic (24.4%). The position of the peak was in the exon of one gene, and this position probably was in the intron of another gene. Therefore, one peak probably corresponded to multiple positions, which resulted in the high proportion of intron functional elements in CHIP-Seq data.
WRINKLED1 of A. thaliana was bound to AW-box sequence [CnTnG](n) 7 [CG], which was conserved among proximal upstream region of the genes involved in fatty acid biosynthesis such as Pl-PKβ1, KAS1, BCCP2, and SUS2 (Maeo et al., 2009). The results indicated that the conserved [CnTnG] and [CG] core motifs were important for WRI1 binding (Maeo et al., 2009). In our present study, two motifs (TGTGTGTGTGTG and ACACACACACAC) were predicted by the Motif Alignment and Search Tool (E-value less than 10, position p-value less than 0.0001) (Bailey et al., 2009). Two identified motifs in our present study were different from the AW-box sequence (Maeo et al., 2009), which would be favorable to study the regulation mechanism of DpWRI1-like in D. parva. The binding of DpWRI1-like to identified motifs will be examined by electrophoretic mobility shift assay in the future.
In addition, many genes regulated by DpWRI1-like were identified. The complete gene information is shown in Supplementary Table 1. Some important genes regulated by DpWRI1-like were discussed in detail (Table 3). In addition, the expression levels of some target genes are also shown in Table 3 based on our previous study (Shang et al., 2016). WRINKLED1 acted as a transcriptional enhancer of PKp-β1 and BCCP2, two genes encoding proteins of glycolytic and fatty acid biosynthetic pathways, respectively in Arabidopsis (Baud et al., 2009). The key regulatory elements were identified in the promoters of some target genes of WRI1, including BCCP2 encoding biotin carboxyl carrier protein, PKp-β1 encoding pyruvate kinase beta subunit 1, genes encoding PDH-E1β and PDH-E1α subunits of pyruvate dehydrogenase (PDH), and PhoGM encoding phosphoglycerate mutase (Baud et al., 2007(Baud et al., , 2009. Then there was the occurrence of a 15-bp motif (consensus sequence: CAAAAGTAGGGGTTT) identified within the above five promoter sequences of target genes of WRI1 (Baud et al., 2007(Baud et al., , 2009. A set of genes encoding plastidic enzymes related to late glycolysis or de novo fatty acid biosynthesis, such as pyruvate kinase, PDH E1α subunit, acetyl CoA carboxylase BCCP2 subunit, an enoyl-ACP reductase, and a phosphoglycerate mutase, were the target genes of WRI1 in Arabidopsis (Baud et al., 2007). LAS encoding lipoic acid synthase and BIO2, which catalyzed the last steps of lipoic acid and biotin biosynthesis, was also assumed as target gene of WRI1 in Arabidopsis (Baud et al., 2007). In our present study, the above target genes encoding pyruvate kinase and  T indicated sample D_Parva_T, C indicated sample D_Parva_C. The log 2 (T/C) indicated the difference of gene expression in mRNA level between D_Parva_T and D_Parva_C. Up or down indicated compared with D_Parva_C, mRNA level was upregulated or downregulated in D_Parva_T. The ID indicated the corresponding gene ID in C. reinhardtii. The expression data originated from our previous work (Shang et al., 2016).
BCCP subunit were also found. In addition, four genes PFK2, PFK3, FBA3, and PGK2, encoding phosphofructokinase family protein, phosphofructokinase, fructose-bisphosphate aldolase 1, and phosphoglycerate kinase in glycolysis pathway were identified as target genes of DpWRI1-like (Table 3). In these four enzymes, phosphofructokinase family protein and phosphofructokinase encoded by PFK2 and PFK3 were key enzymes of the glycolysis pathway. Transient expression of LUC reporter genes using the proximal sequences upstream from the start codon ATG of genes for a subunit of pyruvate kinase (Pl-PKβ1), a subunit of acetyl-CoA carboxylase (BCCP2), and ketoacyl-acyl carrier protein synthase (KAS1) in protoplasts was upregulated by co-expression of WRI1, which was confirmed by in vitro binding of WRI1 to the above upstream sequences (Maeo et al., 2009). In our present study, KAS1 gene associated with fatty acid biosynthetic pathway from D. parva was also identified as target gene of DpWRI1like, which was consistent with the result of Maeo et al. (2009). In addition, a set of genes related to lipid metabolism such as SQD2, DES6, CGLD24, PLAP5, PLAP8 and SSD3, were identified as target genes of DpWRI1-like (Table 3). WRI1 promoted the carbon flow to lipid during seed maturation through directly activating genes involved in fatty acid synthesis and controlling genes related to assembly and storage of triacylglycerol (TAG) in A. thaliana (Maeo et al., 2009). Our present results were consistent with the inference of Maeo et al. (2009).
DES6 catalyzes the desaturation of monoenoic to dienoic acids in chloroplasts. The mRNA level of DES6 increased under nitrate limitation condition, especially at 0.11 mM NaNO 3 in Messastrum gracile SE-MC4 (Anne-Marie et al., 2020). However, our previous study (Shang et al., 2016) found that the expression of DES6 decreased under nitrogen limitation condition (0.5 mM NaNO 3 ) in D. parva. The difference might be related to the species of microalgae. Some target genes (SQD2, CGLD24, PLAP5, PLAP8, SSD3, and KAS1) were not found as DEGs under nitrogen limitation condition in our previous study (Shang et al., 2016), which might be related to the requirement of transcriptome analysis for DEGs (≥ 2 fold change in mRNA level) or the imperfection of transcriptome sequencing. KAS1 encodes a fatty acid elongase (3-ketoacyl-CoA synthase) affecting wax biosynthesis in A. thaliana (Todd et al., 1999). Diacylglycerol acyltransferase encoded by CGLD24 catalyzes the last step in TAG biosynthesis, specifically, the esterification of fatty acid with diacylglycerol. KAS1 and CGLD24 are important genes, and further research will be needed in the future.
In our present study, a set of genes (GWD2, AMYB1, STA6, SSS3, SBE1, PUL1, ISA3) encoding R1 protein/alphaglucan water dikinase, beta-amylase, glucose-1-phosphate adenylyltransferase, soluble starch synthase, starch synthase III, starch branching enzyme, pullulanase-type starch debranching enzyme, isoamylase/starch debranching enzyme from D. parva (Table 3) were identified as target genes of DpWRI1-like, which were related to the synthesis and degradation of starch. The former studies suggested that regulating genes coding enzymes of starch metabolism, such as the small subunit of ADP-glucose pyrophosphorylase, could efficiently divert carbon from starch to TAG biosynthesis (Sanjaya et al., 2011;Marchive et al., 2014). Considering the close relationship between starch and lipid metabolism, these target genes for starch metabolism might play an important role in the regulation of oil biosynthesis. Meantime, many genes related to photosynthesis (such as psaA, psaB, psbB, psbD, PRK1, rbcL, LCI30, and MBB1) were identified as target genes ( Table 3). These genes, which are related to carbohydrate production, could also play an important role in the regulation of oil biosynthesis.
Nicotinamide adenine dinucleotide phosphate (NADP) malic enzyme encoded by MME5 catalyzes the reduction of NADP + to NADPH (Chang and Tong, 2003). The increased NADPH yield, as a result of malic enzyme overexpression, was utilized by enzymes associated with lipid biosynthesis, which led to increased lipid production (Wynn et al., 1999). However, our previous study (Shang et al., 2016) showed that the expression of MME5 decreased under nitrogen limitation condition in D. parva. Of course, there were many factors leading to the increase of oil production. In storage roots of GWD-RNAi lines of cassava, maltose content dramatically decreased and starch with lower phosphorylation level drastically reduced β-amylolytic rate (Zhou et al., 2017). The mRNA level of GWD2 increased under nitrogen limitation condition in D. parva (Table 3), which might contribute to oil biosynthesis. The β-amylase encoded by AMYB1 produces maltose from starch. The decrease of AMYB1 expression level (Table 3) negatively affected oil biosynthesis in D. parva. Glucose-1-phosphate adenylyltransferase encoded by STA6 catalyzes the rate-limiting step in the biosynthesis of glycogen. MaSSIII-1 (SSIII from Musa acuminata) was a key gene in banana fruit, and it catalyzed amylopectin biosynthesis (Miao et al., 2017). The decrease of STA6 and SSS3 from D. parva might contribute to oil biosynthesis because of competitive relationship between oil and glycogen/amylopectin. Phosphoglycerate kinase encoded by PGK2 catalyzes the formation of ATP from ADP and 1,3diphosphoglycerate. Phosphofructokinase encoded by PFK is the rate-limiting enzyme of glycolysis pathway. As shown in Table 3, the decrease of expression levels of PFK2/PFK3/PGK2 might cause the decrease of acetyl CoA (the substrate of fatty acid biosynthesis), which might lead to the decrease of oil biosynthesis.
As shown in Table 3, the expression of four genes (psaA, psbB, psaB, and psbD) related to photosystem was greatly upregulated under nitrogen limitation condition in D. parva, which contributed to photosynthesis, leading to the increase of oil biosynthesis because photosynthesis was the source of oil biosynthesis (Azarin et al., 2020). To confirm the expression levels of the above-mentioned genes, real-time PCR was conducted (Supplementary Figure 5). The results showed that the expression levels of four genes were significantly upregulated under nitrogen limitation (D-Parva_T, all p-values < 0.01). The increase of rbcL expression level (Table 3) might help to increase dark reaction of photosynthesis, leading to the increase of oil biosynthesis. Low-CO 2 -inducible protein (encoded by LCI) played an important role in carbon-concentrating mechanism, and could increase inorganic carbon uptake and intracellular inorganic carbon accumulation under low CO 2 conditions (Atkinson et al., 2016;Kono and Spalding, 2020). As shown in Table 3, the decrease of LCI3 expression might weaken photosynthesis and oil biosynthesis.
In addition, four genes, RWP9, TF2F, RWP11, and TF2H2, encoding RWP-RK TF, TF IIF, RWP-RK TF, and the general TF IIH subunit, were identified as target genes of DpWRI1-like. WRI1 target genes, altering other physiological processes such as hormone homeostasis, were important for the application of WRI1 to increase oil content (Kong et al., 2017;Kong and Ma, 2018b). These TFs and their target genes remained unknown, which excited our interest.
The enrichment of several important genes (PK7, PK8, PK11, KAS3, and BCCP) identified by ChIP was tested to check the reliability of ChIP-Seq. All DNA fragments corresponding to these five genes could be detected. In addition, the enrichment was more significant in D_Parva_T compared with D_Parva_C (Figure 8  In summary, a hypothesis on the regulation of lipid biosynthesis by DpWRI1-like was established in D. parva. DpWRI1-like regulated target genes related to carbohydrate metabolism (including photosynthesis, glycolysis pathway, tricarboxylic acid cycle, starch metabolism) to enhance the conversion of sugar into lipid. DpWRI1-like regulated target genes related to lipid metabolism (especially fatty acid biosynthesis such as KAS gene) to promote the production of lipid.

Validation of Target Gene MME5
The interaction of DpWRI1-like and the promoter of MME5 from D. parva was checked by Yeast One-Hybrid technology. Target sequence (bait, MME5 promoter) was cloned into pAbAi vector to form pBait-AbAi vector. Then the linearized pBait-AbAi was integrated into the ura3-52 locus of Y1H Gold to create Y1HGold[Bait] reporter strain. The prey protein encoded by DpWRI1-like cDNA was screened by cotransformation (with pGADT7-Rec) and in vivo homologous recombination in Y1HGold[Bait] yeast according to the growth of transformed  cells on SD/-Leu/+AbA plate. The linearized pAbAi vector was used as the negative control to eliminate false positive condition (transformed cells could not grow on SD/-Leu/+AbA plate). Figure 9 showed that PCR amplification by 30a-wri1(CF) and 30a-wri1(CR) using positive yeast colony as template had the expected size (about 1 kb). After sequencing, the fragment was DpWRI1-like cDNA. The result confirmed the interaction of DpWRI1-like and MME5, TF activity of DpWRI1-like, and ChIP analysis output.

CONCLUSION
WRINKLED1 is a key regulator of oil biosynthesis. The metabolic pathways involved in the biosynthesis and degradation of fatty acid, TAG, and starch were constructed, and Dpwri1-like encoding TF DpWRI1-like had been identified in D. parva in our former study. However, target genes of DpWRI1like remain unknown. DpWRI1-like targets were identified by ChIP-Seq to elucidate the function of DpWRI1-like in this study. The results showed that many target genes related to carbohydrate metabolism, lipid metabolism, photosynthesis, and transcription factor were identified. A hypothetical diagram about the mechanism of oil accumulation induced by DpWRI1like under nitrogen limitation condition was constructed based on our present study (Figure 10). This study laid a good foundation for further understanding of regulatory mechanism of oil biosynthesis related to DpWRI1-like.

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: NCBI [accession: PRJNA698763].

ETHICS STATEMENT
About ethics, the experimental material was microalga Dunaliella parva FACHB-815 which was purchased from Freshwater Algae Culture Collection at the Institute of Hydrobiology. In addition, this article had no academic misconduct. Therefore, this research had no problem about ethics.

AUTHOR CONTRIBUTIONS
CS contributed to the conception, design, data acquisition, and drafting of the article. BP and HY contributed to data acquisition and drafting of the article. SG and YL contributed to drafting of the article. All authors read and approved the final manuscript.