Team:Evry Paris-Saclay/Model
Computational analyses and Modeling were very important parts of our project as they allowed us to address two fundamental questions:
(1) How can we produce jacaric acid in the absence of any information on the enzyme catalysing its synthesis?
(2) Which Yarrowia lipolytica strain is the most adapted / efficient / robust for the Conjugated Linolenic Acids (CLnAs) production like punicic and jacaric acids?
To find answers to the above questions, we (i) undertook a substantial experimental and computational effort to identify the sequence of the FadX gene directly from Jacaranda seeds, and (ii) carried out Flux Balance Analysis to identify the Yarrowia lipolytica strains best suited for the production of rare fatty acids of our interest. We describe our work below in detail.
Uncovering the Jacaranda mimosifolia FadX gene
In the scope of the iGEM project, we needed to identify and characterize the enzyme responsible for catalyzing the synthesis of jacaric acid, a Conjugated Linolenic Acid (CLnA) with anti-tumoral, anti-inflammatory and anti-obesity properties. This compound is naturally found at a high concentration in the seeds of the plant Jacaranda mimosifolia (hereafter called Jacaranda), however, the genome of this species is not sequenced. To identify which enzyme would catalyze such reaction in Jacaranda, we decided to perform Exome-sequencing on both fresh and germinated seeds.
RNA sequencing of Jacaranda mimosifolia and data pre-processing
We performed RNA extraction for the fresh and germinated seeds of Jacaranda using the NucleoSpin® RNA Plant and Fungi kit (Macherey Nagel) and, after a quick analysis on an ethidium bromide stained agarose gel followed by a Qubit RNA assay (Thermo Fisher Scientific) samples were prepared for sequencing using the TruSeq Stranded mRNA kit (Illumina). The libraries were analyzed on a HiSeq 4000 system (Illumina) and the obtained data were preprocessed and cleaned by the sequencing platform (Genoscope).
Illumina Filter
Raw reads were filtered to remove clusters that had too many bases with ambiguous signal. The purity of the signal was analyzed on the first 25 cycles of each cluster by assigning a score to each of the cycles: Chastity = (Maximum Intensity) / (Sum of the two highest intensities). The filter used in basic calling allows at most a cycle with a chastity value of less than 0.6.
Genoscope Filter
Adapters and primers sequences were eliminated from the readings as well as low-quality bases (Q <20) at 2 extremities. Reads of less than 30 nucleotides were discarded. Finally, PhiX reads from Illumina internal spiking were discarded.
Filtering of ribosomal sequences
The cleaned reads corresponding to ribosomal RNA were separated from other reads.
Methods
For the analyses of the exome sequencing, we decided to follow two independent approaches: alignment to a reference genome and de novo transcriptome assembly.
Alignment using a reference genome
In this first approach, we used the genome of Handroanthus impetiginosus [1], a Bignoniaceaea plant as a reference genome for the alignment of the reads. We choose it following the discussions we had with Dr. Florian Jabbour, senior lecturer and collection manager in the field of morpho-anatomy and plant development at the Institute of Systematics, Evolution, Biodiversity of the National Museum of Natural History of Paris (for details, see our Human Practices page on this wiki). We used the Galaxy platform [2] to analyze the fastq files provided by the Illumina platform. We uploaded four fastq files containing the paired-end sequences for the germinated and fresh seeds (two for each seed type) as well as the reference genome [1]. First, we used FastQ Groomer to ensure that our fastq files used standard quality format. Next, we used FastQC to perform a quality control of our files. Before mapping the obtained sequences against a reference genome, we visualized the statistics corresponding to the raw data for each paired-end sequence, we did so for the fresh seeds (referred to as COS) and the germinated ones (BOS).
We then proceed to the mapping of our different fastq files (alignment of the reads to a reference genome) creating bam files using BWA tool and obtained bam files for both seed exomes. We sorted those files according to their coordinates on the reference genome. At this step, we used bcftools mpileup tool to create a VCF file with the variant calls. We then used bcftools consensus in order to create a consensus sequence, substituting the reference genome bases for the variants of the VCF file. Finally, the consensus sequences obtained for the fresh seeds and germinated seeds were translated into peptide sequences according to three forward frames using EMBOSS Transeq [3]. The corresponding peptide sequences were then aligned along with nine sequences of FadX enzymes, five Fad2 enzymes and the two peptide sequences of Handroanthus impetiginosus with the best similarity.
De novo transcriptome assembly
In this second approach, we performed a de novo transcriptome assembly from Illumina reads of 150 nucleotides length to obtain longer contigs which would cover the entire sequence of our enzyme. We then aligned those transcripts to the FadX enzymes (table 1) to identify the closest transcripts. We performed the de novo transcriptome assembly using the galaxy platform [1]. Briefly, we used the Trimmomatic tool [4] with a sliding window of four nucleotides to ensure the quality of the reads. Then we used the Trinity tool [5] on the paired reads to create longer contigs and obtain two libraries of roughly 135000 and 200000 sequences (assembly quality reports: germinated seeds; fresh seeds.
As the genome of Jacaranda is not annotated, it was impossible for us to identify the start codon of those sequences. Moreover, as we cannot be sure that our contigs contain the beginning of the real transcript and therefore the start codon, we translated all the transcripts of the library into amino acid sequences according to the three potential open reading frames (ORFs). To convert the transcripts to protein sequences and to perform local alignments, we used Python’s module Seq from the package Bio [6]. For the translation, we used the default table for plant plastids provided by the module. For the alignments, we used a BLOSUM62 transition matrix. We first performed the analysis using parameters with a gap open penalty of 10 and a gap extension penalty of 0.5. We then refined our results with a more stringent gap open penalty of 20. Python script is available here.
To analyze our alignments, we ranked the scores obtained by each contig with each enzyme for one ORF and then compare the results between contigs and ORFs.
Results and comparison the two approaches
We used a set of nine known CLnA enzymes (table 1) to perform local alignments and identified the closest sequences. We then challenged those alignments with five sequences corresponding to Fad2 enzymes (table 2) performing multiple alignment with clustal omega [7] and constructed phylogenetic trees using phylML 3.0 (bootstraping parameter: 100) [8] to ensure that the sequence discovered would correspond to a FadX enzyme. We visualized the multiple alignments with Jalview [9].
| Enzyme | Organism | Uniprot identifier |
|---|---|---|
| Bifunctional fatty acid conjugase/Delta(12)-oleate desaturase | Punica granatum | Q84UB8 |
| Bifunctional fatty acid conjugase/Delta(12)-oleate desaturase | Trichosanthes kirilowii | Q84UC0 |
| Bifunctional desaturase/conjugase FadX | Vernicia fordii | Q8GZC2 |
| Delta13 fatty acid desaturase FadX-1B | Momordica charantia | Q9SP61 |
| Delta(12) acyl-lipid conjugase (11E,13E-forming) | Impatiens balsamina | Q9SP62 |
| Delta13 fatty acid desaturase FadX-1B | Exocarpos cupressiformis | U5LN76 |
| Fatty acid conjugase Fac2 B | Calendula officinalis | Q9FPP7 |
| Fatty acid conjugase Fac2 A | Calendula officinalis | Q9FPP8 |
| Delta(12) fatty acid desaturase DES8.11 | Calendula officinalis | Q9SCG2 |
| Enzyme | Organism | Uniprot identifier |
|---|---|---|
| Delta(12) fatty acid desaturase Fad2 | Calendula officinalis | Q9AT72 |
| Delta(12)-acyl-lipid-desaturase | Punica granatum | Q84VT2 |
| Delta(12)-fatty-acid desaturase | Arabidopsis thaliana | P46313 |
| Delta(12)-fatty-acid desaturase Fad2 | Vernicia fordii | Q8GZC3 |
| (submitted) | Yarrowia lipolytica | Q6CF55 |
Results for the alignment to a reference genome
For comparison, we used the genes predicted for Handroanthus impetiginosus (family Bignoniaceae) and the peptide sequences of the nine CLnA enzymes (table 1). Searching for local alignments with Python’s Bio::seq module (script available here), we identified two gene sequences of Handroanthus impetiginosus (Haimp10002232m and Haimp10040446m) that shared the best similarity with those enzymes (Figure 8). As a preliminary analysis, we performed multiple alignments between the sequences of enzymes in Tables 1 (FadX) and 2 (Fad2), the consensus sequences for the fresh (COS) and germinated (BOS) seeds of Jacaranda (for the three reading frames) and the two peptide sequences of Handroanthus impetiginosus. These alignments showed that the sequence translated with the first reading frame of the fresh seeds was the one sharing most similarity with FadX enzymes while other ORFs as well as the Handroanthus impetiginosus gene sequences would share more similarity with Fad2 enzymes or constitute outliers. We, therefore, decided to focus on this sequence and repeated the alignment with the FadX and Fad2 enzymes constructing a phylogenetic tree to see which enzyme would be closer to our sequence (Figure 8).
The consensus sequence for the first reading frame of the fresh seeds being much longer than the ones of the enzymes, we extracted the sequence where the alignment takes place and repeated the analysis (Figure 9).
We can see that the sequence extracted from the fresh seed is localized between FadX and Fad2 enzymes without really clustering with any specific type.
Results from the de novo transcriptome assembly
Performing local alignment between contigs and FadX enzymes (Table 1), we first filter our results depending on their score. For the fresh seeds, tables 3, 4 and 5 summarize the top-30 contigs for each enzyme according to the three ORFs.
We chose an arbitrary threshold of 450 which corresponds to a minimal level of similarity between the sequences. We selected one contig from the first ORF: DN4193_c0_g2_i1 and four contigs from the third ORF: DN16194_c0_g1_i1, DN17966_c0_g1_i1, DN4193_c0_g1_i1 and DN17966_c0_g1_i2.
For the germinated seeds, we obtained the following rankings for the contigs (tables 6, 7 and 8).
As previously, we used an arbitrary threshold of 450 and selected six contigs for validation, three from the first ORF: TRINITY_DN16261_c0_g2_i1, TRINITY_DN18544_c0_g1_i2 and TRINITY_DN18544_c0_g1_i1; one from the second ORF: TRINITY_DN63614_c0_g1_i1 and two from the third ORF: DN16982_c0_g1_i1 and DN16982_c0_g1_i2 which correspond to two isoforms of the same gene.
We aligned the different contigs selected together with the different FadX and Fad2 enzymes (contigs and enzyme sequences). We then constructed a phylogenetic tree (Figure 11).
The selected contigs seem to cluster preferentially with Fad2 enzymes. We tried to reduce the selection to the top selected contigs (first ranked for each ORF) only and obtain the phylogenetic tree shown in Figure 12.
Discussion of the two approaches
Both approaches, i.e. using a reference genome or performing de novo assembly, led us to the discovery of peptide sequences sharing similarities with the FadX enzymes.
We observe that the consensus sequence (COS_Seq) identified through alignment to a reference genome (Handroanthus impetiginosus) clusters with the contigs identified from the de novo transcriptome assembly which obtained the higher score among all. The reference genome approach seems computationally more efficient but also depends on the quality of the reference chosen. Here, our reference was partly annotated and we had to explore the different ORFs to search for a specific type of enzyme. On the other hand, the de novo assembly approach does not require prior knowledge but the results are harder to interpret. To conclude, using both bioinformatic approaches we were able to identify sequences similar to both FadX and Fad2 enzymes. However, for determining if the sequences identified code for either type, we should perform further analyses, both computational and experimental. Since Fad2 and FadX sequences are very similar as seen in the phylogenetic tree (Figure 6) and the alignment (Figure 7), it is difficult to reach a confident conclusion only considering the results of our analyses. To go further, we could also try to combine both approaches and use the reference genome to facilitate the assembly of the contigs identified in the de novo analysis.
Experimental validations were launched and, as a first step, the Yarrowia lipolytica codon optimized version of a putative bifunctional fatty acid conjugase / desaturase (FadX) from Jacaranda mimosifolia was added to iGEM Registry (BBa_K2983063).
Flux Balance Analysis
In order to choose the best Yarrowia lipolytica strain to produce rare fatty acids such as jacaric acid and punicic acid, we have performed metabolic simulations, also called Flux Balance Analysis (FBA), on several known strains of Y. lipolytica.
First, we have got an improved version of iYali4 [10], a metabolic model of Yarrowia lipolytica W29, from GitHub.
We have used the Python package COBRApy [11] to add a simplified version of the reaction transforming linoleate into jacaric acid. By clicking below, you can download the original model and the code we have created to add the reaction.
Then, we have performed an FBA on this model, which corresponds to the JMY195 [12] strain, with the software OptFlux [13]. This simulation gave us a biomass value of 29.5 and a jacaric acid production of 226.54. These values are specific to the software and cannot be compared to the experimental yields directly. However, these values allow us to compare different strains. We have constructed models of the strains JMY1877 [14], JMY1233 [15], JMY2159 [16], JMY3325 [17] and JMY3820 [18] and performed FBA on them with the OptFlux tool “Gene Under-Over Expression Simulation”. In such models, the normal value for gene expression in set to 1. To knockout genes in the models, we have set this value to 0, and we have set it to 10 to simulate an over-expression. The FBA results are presented in Table 3.
| JMY195 (Po1d) | MATA ura3-302 leu2-270 xpr2-322 | 29.50 | 100.0% | 226.54 | 100.0% |
| JMY1233 | MATA ura3-302 leu2-270 xpr2-322 pox1-6Δ | 29.50 | 100.0% | 226.54 | 100.0% |
| JMY1877 | MATA ura3-302 leu2-270 xpr2-322 dga1Δ dga2Δ lro1Δ are1Δ | 0.00 | 0.0% | 0.00 | 0.0% |
| JMY2159 | MATA ura3-302 leu2-270 xpr2-322 pox1-6Δ dga1Δ dga2Δ lro1Δ fad2Δ | 0.00 | 0.0% | 0.00 | 0.0% |
| JMY3325 | MATA ura3-302 leu2-270 xpr2-322 pox1-6Δ dga1Δ dga2Δ lro1Δ fad2Δ pTEF-FAD2-LEU2 | 0.00 | 0.0% | 0.00 | 0.0% |
| JMY3820 | MATα ura3-302 leu2-270 xpr2-322 pox1-6Δ tgl4Δ + pTEF-DGA2 + pTEF-GPD1 | 29.50 | 100.0% | 386.41 | 170.6% |
For JMY1877, JMY2159 and JMY3325, the results are not relevant since their models were unable to show any growth. This is certainly due to a limitation in the original model, which treats as essential some genes that we have deleted. Indeed, when we have determined critical genes with OptFlux, are1 appeared to be essential, even if it is not in real life. About JMY2159 and JMY3325, it is likely that dga1, dga2 and lro1 are essential when deleted together, which again is not the case as seen experimentally.
For JMY1233, we see that both the biomass and the jacaric acid production have the same value than for the control JMY195.
Finally, the FBA indicates that JMY3820 is able to produce 70% more jacaric acid than JMY195, without changing the biomass production.
These results suggest that JMY3820 is the best strain to produce jacaric acid from linoleate. This would also suggest that it is the best strain to produce punicic acid with linoleate as the starting point, and we kept that in mind during the wet lab experiments. Happily, this has been confirmed by the experiments.
To go further, we have also tried to determine new genes that could be deleted to optimize the production of jacaric acid with Yarrowia lipolytica. To do so, we have used the evolutionary optimization tool of OtpFlux on JMY195 and set the minimum biomass production to 95% of the control. As results, this tool gave us genes or sets of genes the deletion of which increases the jacaric acid production to 1000, which is the maximum potential value of the reaction. These genes are presented in Table 4. The solutions proposed by this tool should always be manually verified since a model may not be refined enough to represent the exact reality and could propose absurd deletions as a result.
| YALI0D00583g | |
| YALI0F15631g | |
| YALI0E33517g | |
| YALI0F24475g | |
| YALI0E25740g | |
| YALI0F01210g | |
| YALI0F26323g | |
| YALI0E18568g | |
| YALI0A05379g | |
| YALI0B09647g | |
| YALI0B22682g | |
| YALI0E33099g | |
| YALI0C05258g | |
| YALI0D04741g | |
| YALI0D10131g | |
| YALI0E31009g | |
| YALI0F14025g | |
| YALI0E31471g | |
| YALI0D24750g | |
| YALI0E24013g | |
| YALI0C11407g |
YALI0D00583g, YALI0F15631g and YALI0F24475g are coding for a proton-transporting protein responsible for the vacuolar acidification; YALI0E33517g is coding for an oxoglutarate dehydrogenase involved in the Krebs cycle; YALI0E25740g is coding for a guanine deaminase involved in the reaction EC:3.5.4.3; YALI0F01210g is coding for a membrane water transporter; etc.
Interestingly, YALI0F14025g is the only gene that was not proposed in a set a gene. The only deletion of this gene is able to increase the production of jacaric acid to its maximum. This gene is coding for a kinase that is involved in the inositol-phosphate biosynthesis.
None of the deletions that have been performed in tested strains, including tgl4 which is deleted in JMY3820, was proposed by OptFlux.
This modelling part has implications well beyond its immediate scope and will definitely help future studies working on optimizing Yarrowia lipolytica for rare fatty acids production.
References
[1]Silva-Junior OB, Grattapaglia D, Novaes E, Collevatti RG. Genome assembly of the Pink Ipê (Handroanthus impetiginosus, Bignoniaceae), a highly valued, ecologically keystone Neotropical timber forest tree. Gigascience (2018) 7, 1-16.
[2]Afgan E, Baker D, Batut B, van den Beek M, Bouvier D, Cech M, Chilton J, Clements D, Coraor N, Grüning BA, Guerler A, Hillman-Jackson J, Hiltemann S, Jalili V, Rasche H, Soranzo N, Goecks J, Taylor J, Nekrutenko A, Blankenberg D. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Res (2018) 46(W1), W537-W544.
[3]Madeira F, Park YM, Lee J, Buso N, Gur T, Madhusoodanan N, Basutkar P, Tivey ARN, Potter SC, Finn RD, Lopez R. The EMBL-EBI search and sequence analysis tools APIs in 2019. Nucleic Acids Res (2019) 47(W1), W636-W641.
[4]Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics (2014) 30, 2114-2120.
[5]Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol (2011) 29, 644-652.
[6]Cock PJ, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, Friedberg I, Hamelryck T, Kauff F, Wilczynski B, de Hoon MJ. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics (2009) 25, 1422-1423.
[7]Sievers F, Wilm A, Dineen D, Gibson TJ, Karplus K, Li W, Lopez R, McWilliam H, Remmert M, Söding J, Thompson JD, Higgins DG. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol (2011) 7, 539.
[8]Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol (2010) 59, 307-321.
[9]Clamp M, Cuff J, Searle SM, Barton GJ. The Jalview Java alignment editor. Bioinformatics (2004) 20, 426-427.
[10]Kerkhoven EJ, Pomraning KR, Baker SE, Nielsen J. Regulation of amino-acid metabolism controls flux to lipid accumulation in Yarrowia lipolytica. NPJ Syst Biol Appl (2016) 2, 16005.
[11]Ebrahim A, Lerman JA, Palsson BO, Hyduke DR. COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Syst Biol (2013) 7, 74.
[12]Barth G, Gaillardin C. Yarrowia lipolytica. In: Wolf K (ed) Non conventional yeasts in biotechnology. Springer, Berlin (1996) 1, 314-388.
[13]Rocha I, Maia P, Evangelista P, Vilaça P, Soares S, Pinto JP, Nielsen J, Patil KR, Ferreira EC, Rocha M. OptFlux: an open-source software platform for in silico metabolic engineering. BMC Syst Biol (2010) 4, 45.
[14]Beopoulos A, Haddouche R, Kabran P, Dulermo T, Chardot T, Nicaud JM. Identification and characterization of DGA2, an acyltransferase of the DGAT1 acyl-CoA:diacylglycerol acyltransferase family in the oleaginous yeast Yarrowia lipolytica. New insights into the storage lipid metabolism of oleaginous yeasts. Appl Microbiol Biotechnol (2012) 93, 1523-1537.
[15]Beopoulos A, Mrozova Z, Thevenieau F, Le Dall MT, Hapala I, Papanikolaou S, Chardot T, Nicaud JM. Control of lipid accumulation in the yeast Yarrowia lipolytica. Appl Environ Microbiol (2008) 74, 7779-7789.
[16]Beopoulos A, Verbeke J, Bordes F, Guicherd M, Bressy M, Marty A, Nicaud JM. Metabolic engineering for ricinoleic acid production in the oleaginous yeast Yarrowia lipolytica. Appl Microbiol Biotechnol (2014) 98, 251-262.
[17]Imatoukene N, Verbeke J, Beopoulos A, Idrissi Taghki A, Thomasset B, Sarde CO, Nonus M, Nicaud JM. A metabolic engineering strategy for producing conjugated linoleic acids using the oleaginous yeast Yarrowia lipolytica. Appl Microbiol Biotechnol (2017) 101, 4605-4616.
[18]Lazar Z, Dulermo T, Neuvéglise C, Crutz-Le Coq AM, Nicaud JM. Hexokinase--A limiting factor in lipid production from fructose in Yarrowia lipolytica. Metab Eng (2014) 26, 89-99.
Our Sponsors
Contact Us
