Team:Evry Paris-Saclay/NGS

Title


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).

Figure 1. Ethidium bromide stained agarose gel showing the bands from fresh and germinated Jacaranda seeds RNA extracts.

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.


Figure 2. Galaxy workflow for the alignment to a reference genome approach.

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).


Figure 3. Galaxy workflow for the de novo transcriptome assembly approach.

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].


Table 1. FadX enzymes.
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


Figure 4. Multiple alignments of the FadX enzymes with colors representing the identity percentage.

Table 2. Fad2 enzymes.
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

Figure 5. Multiple alignments of the Fad2 enzymes with colors representing the identity percentage.

Figure 6. Phylogenetic tree of the Fad2 and FadX enzymes.

Figure 7. Multiple alignments of the FadX and Fad2 enzymes with colors representing the identity percentage.

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).

Figure 8. Phylogenetic tree with the sequences from the first ORF of the fresh (COS) seeds in the alignment.

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).


Figure 9. Phylogenetic tree with the sequences from the first ORF of the fresh seeds. Within the first reading frame, we selected the sequence corresponding to the amino acids 91623 to 91980.

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.


Figure 10. Multiple alignments between the selected consensus sequence (amino acids 91623 to 91980) from the first ORF of the fresh (COS) seeds and the FadX and Fad2 enzymes. The colors represent the alignment score for a BLOSUM62 transition matrix.

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.


Table 3. Top-30 results of local alignments between contigs of the fresh seeds (COS) for the first ORF and each FadX enzyme (table 1).

Table 4. Top-30 results of local alignments between contigs of the fresh seeds (COS) for the second ORF and each FadX enzyme (table 1).

Table 5. Top-30 results of local alignments between contigs of the fresh seeds (COS) for the third ORF and each FadX enzyme (table 1).

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).


Table 6. Top-30 results of local alignments between contigs of the germinated seeds (BOS) for the first ORF and each FadX enzyme (table 1).

Table 7. Top-30 results of local alignments between contigs of the germinated seeds (BOS) for the second ORF and each FadX enzyme (table 1).

Table 8. Top-30 results of local alignments between contigs of the germinated seeds (BOS) for the third ORF and each FadX enzyme (table 1).

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).


Figure 11. Phylogenetic tree from the selected contigs and the FadX and Fad2 enzymes.

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.


Figure 12. Phylogenetic tree considering only the top contigs and the FadX and Fad2 enzymes.

Figure 13. Multiple alignments between the top contigs and the FadX and Fad2 enzymes colors representing the alignment score for a BLOSUM62 transition matrix.


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.


Figure 14. Multiple alignments between the selected consensus sequence (amino acids 91623 to 91980) from the first ORF of the fresh (COS) seeds identified through alignment to a reference genome and the selected contigs identified through de novo transcriptome assembly. The colors represent the alignment score for a BLOSUM62 transition matrix.
Figure 15. Phylogenetic tree considering the selected consensus sequence (amino acids 91623 to 91980) from the first ORF of the fresh (COS) seeds identified through alignment to a reference genome and the selected contigs identified through de novo transcriptome assembly.

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).


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.