Transcriptome analysis of Chelidonium majus elaiosomes and seeds provide insights into fatty acid biosynthesis

Plant materials and treatments

The sampled plants of C. majus were grown in Beijing Badaling National Forest Park, China, under natural conditions (Fig. 1A). To gain a comprehensive transcriptional profile of its elaiosomes and seeds, sample collection at the optimal developmental stage is crucial. Based on previous observations of C. majus elaiosome and seed development, four different developmental stages were defined (Yang et al., 2015). In the present study, we chose these four developmental stages, namely flower-bud stage (pre-fertilization) (Ch01) (Fig. 1B), flowering stage (post-pollination) (Ch02) (Fig. 1C), young seed stage (Ch03) (Fig. 1E), and mature seed stage (Ch04) (Fig. 1F) for comparative transcriptome analysis to explore gene expression during fatty acid biosynthesis in the elaisomes and seeds of C. majus. Corresponding to the four stages, flower buds (with bracts, petals, and stamens removed), fertilized ovaries, young seeds (testa white to brown), and mature seeds (testa black and hard) were collected as RNA-seq materials. They were frozen in liquid nitrogen and stored at −80 °C until further use. Three samples were collected as biological repeats for each developmental stage.

RNA isolation and sequencing

To obtain transcriptome expression profiles, RNA-seq transcriptome libraries were prepared using the TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA, USA) and five μg of total RNA. First, messenger RNA was isolated through the polyA selection method using oligo (dT) beads and fragmented in a fragmentation buffer. Next, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad CA, USA) and random hexamer primers (Illumina, San Diego, CA, USA). Then, the synthesized cDNA was subjected to end-repair, phosphorylation, and “A” base addition according to Illumina’s library construction protocol. Libraries were size selected for cDNA target fragments of 200–300 bp on a 2% Low Range Ultra Agarose gel followed by polymerase chain reaction (PCR) amplification using Phusion DNA polymerase (NEB, Ipswich, MA, USA) for 15 cycles. After quantification using TBS380, the paired-end RNA-seq sequencing library was sequenced using Illumina HiSeq 4000 (2 × 150 bp read length). RNA samples were named Ch01a/b/c (flower-bud stage), Ch02a/b/c (flowering stage), Ch03a/b/c (young seed stage), and Ch04a/b/c (mature seed stage). A total of 12 cDNA libraries were then constructed and sequenced using the Illumina deep sequencing platform.

De novo assembly, annotation, differential expression analysis

The raw paired-end reads were trimmed and quality controlled by SeqPrep and Sickle with default parameters. Then, clean data were used to do de novo assembly with Trinity (Grabherr et al., 2011). All of the assembled transcripts were searched against the NR, Cluster of Orthologous Groups (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using BLASTX to identify the proteins that had the highest sequence similarity to the given transcripts to retrieve their function annotations; a typical cut-off E-value of less than 1.0 × 10⁻⁵ was set. The BLAST2GO program was used to get Gene Ontology (GO) annotations of uniquely assembled transcripts for describing biological processes, molecular functions, and cellular components (Conesa et al., 2005). COG is an early database that is expected to be a useful platform for the functional annotation of newly sequenced genomes (Tatusov et al., 2003). The KEGG database is the primary public repository devoted specifically to pathways. Metabolic pathway analysis was performed using KEGG (Kanehisa & Goto, 2000), and KEGG classifications facilitate the exploration of categories of genes with specific functions. To identify differentially expressed genes (DEGs) between different samples, the expression level of each transcript was calculated according to the fragments per kilobase of the exon per million mapped reads method. RNA-Seq by expectation-maximization was used to quantify gene and isoform abundances (Li & Dewey, 2011). The R statistical package software, empirical analysis of digital gene expression in R, was utilized for the differential expression analysis (Robinson, McCarthy & Smyth, 2010). In addition, functional-enrichment analysis, including GO and KEGG, were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at a Bonferroni-corrected P-value ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools and KOBAS (Xie et al., 2011).

Quantitative real-time PCR analysis

For each sample, one μg of total RNA was used for cDNA synthesis using the Plant RNA Purification Reagent (Invitrogen, Carlsbad CA, USA). ChamQ SYBR Color qPCR Master Mix (2X) (Vazyme Biotechnology, Nanjing, China) was used for Quantitative real-time PCR (qRT-PCR). qRT-PCR was performed on Line gene 9600 plus Fluorescence quantitative PCR instrument (Hangzhou, China). Each sample contains three independent biological replicates.

Transcriptome sequencing and de novo assembly

In this study, a total of 12 libraries were sequenced. After removing the low-quality reads and adaptor sequences, the remaining clean data comprised of approximately 61 Gb of sequence data (171 million reads from the flower-bud stage (Ch01); 141 million reads from flowering stage (Ch02); 148 million reads from young seed stage (Ch03); and 154 million reads from mature seed stage (Ch04)). The results are shown in Table 1. De novo assembly of the clean reads was performed using Trinity to generate a reliable integrated reference sequence database (https://doi.org/10.6084/m9.figshare.8044838.v1, File 1). Some reads comprise of a contig, and some contigs comprise of a unigene. We obtained a total of 63,064 unigenes (https://doi.org/10.6084/m9.figshare.8044838.v1, File 2) for all studied samples; these unigenes had an average length and N50 value of 1,114.56 bp and 1,565 bp, respectively. After sorting the assembled unigenes from large to small, the N50 value was calculated as half of the total length of the unigenes. The size distribution of the assembled transcripts is shown in Fig. 2. As shown in Table 1 and Fig. 2, the sequencing data were of sufficiently high quality for subsequent analyses. Therefore, gene expression data from C. majus elaisomes and seeds were reliable.

Functional annotation and classification of unigenes

To verify the potential functions of the 63,064 assembled unigenes, we performed sequence similarity BLASTx searches against the NR (NCBI protein nonredundant), COG, GO, and KEGG protein databases with E-values < 1e⁻⁵. After alignment, gene annotation, functional classification, and metabolic pathway assignment, a total of 30,143 (47.80%) of the unigenes were successfully annotated (https://doi.org/10.6084/m9.figshare.8044838.v1, File 3), which indicated that their functions were relatively conserved and diverse.

Characterization of unigene sequencing

Among all the 30,143 annotated unigenes, 29,352 (97.38%) were mapped to the NR database. In addition, the distribution of E-values of the hits in the NR database indicated that 20,378 unigenes, representing more than half of the mapped sequences (69.43%), showed significant homology (E-value < 1.0 × 10⁻³⁰), while 8,973 unigenes (about 30.57%) had E-values ranging from 1.0 × 10⁻³⁰ to 1.0 × 10⁻⁵ (Fig. 3; https://doi.org/10.6084/m9.figshare.8044838.v1, File 3). Smaller E-values distribution (Fig. 3A) indicated a higher reliability of hits. The similarity distribution showed that the majority of the mapped 26,160 unigenes (89.13%) had similarity values >60%, whereas only 3,192 (10.87%) had similarity values <40% (Fig. 3B). To identify species specificity, individual unigenes were annotated based on the highest BLAST scores against the NR database. Among flowering plants, 6,666 (22.71%) unigenes showed significant homology to Vitis vinifera, followed by Theobroma cacao (1,983 unigenes, 6.76%), Phoenix dactylifera (1,230 unigenes, 4.19%), Populus trichocarpa (1,224 unigenes, 4.17%), Prunus persica (1,090 unigenes, 3.71%), and Jatropha curcas (1,065 unigenes, 3.63%) (Fig. 3C). A total of 2,980 unigenes, representing around 10.15% of all unigenes, did not match the NR database, and thus provided novel genomic information on C. majus. Some of these unigenes are likely to be untranslated regions or non-coding RNAs.

Functional classification by gene ontology

By conducting GO functional classification, a total of 258 unigenes were classified into one or more terms and divided into 43 functional groups, belonging to three major categories (Fig. 4). In the biological process category, the majority of the assigned proteins were involved in metabolic processes (177 unigenes, 68.60%), followed by genes involved in cellular processes (169 unigenes, 65.50%) and single-organism processes (156 unigenes, 60.47%). Under the molecular function category, the mapped unigenes were categorized into nine groups, out of which genes were predominantly associated with catalytic activity (165 unigenes, 63.95%) and binding (123 unigenes, 47.67%), which exceeded the total number of the other terms. The cellular component category was divided into 14 groups, in which the largest assignments were found to be in the cell and cell part. There were 164 unigenes in the cell and cell part (63.57%), followed by organelles (115 unigenes, 44.57%), membrane (52 unigenes, 20.16%), and membrane parts (36 unigenes, 13.95%). The above results suggested that the development of elaisomes and seeds in C. majus are closely associated with enzyme-catalyzed reactions, cell structure formation, and metabolic processes.

Functional classification by clusters of orthologous groups

To predict the potential function of genes and evaluate the completeness of the transcriptome libraries, all assembled unigenes were searched against the COG database. A total of 7,323 mapped unique sequences were clustered into 25 functional categories (Fig. 5). Among them, the largest category was “general functions prediction only” (917, 12.52%), followed by “signal transduction mechanisms” (798, 10.90%), “posttranslational modification, protein turnover, and chaperones” (727, 9.93%), and “translation, ribosomal structure, and biogenesis” (596, 8.14%). Thus, we inferred that these genes play an important role in ribosome formation, protein synthesis, and signal transduction during the development of elaisomes and seeds in C. majus. Notably, 322 unigenes (4.40%) were classified under the “lipid transport and metabolism” category, which could provide useful information for further studies on lipid metabolism in C. majus elaisomes and seeds.

Biological pathways identification by KEGG

To identify the biological pathways that were actively involved in lipids accumulation in C. majus elaisomes and seeds, all annotated unigenes were mapped to the reference pathways in the KEGG. The analysis identified a total of 11,540 unigenes that were assigned to five clusters and distributed across 373 KEGG pathways (Fig. 6). Most of the unigenes were grouped under the metabolism cluster (A). Within this cluster, the majority of unigenes were assigned under the “global and overview maps pathway” (2,736, 23.71%), followed by “carbohydrate metabolism” (1,022, 8.85%), “amino acid metabolism” (800, 6.93%), and “energy metabolism” (590, 5.11%) categories. The relative abundance of unigenes involved in “amino acid metabolism” pathways is probably related to the biosynthesis of massive storage protein in the C. majus endosperm. Notably, in the “energy metabolism” cluster, some pathways were found to be closely related to the lipid metabolism. A total of 508 unigenes, including glycerophospholipid metabolism (205, 40.35%), glycerolipid metabolism (108, 21.26%), fatty acid degradation (100, 19.69%), biosynthesis of unsaturated fatty acids (93, 18.31%), fatty acid elongation (75, 14.76%), fatty acid biosynthesis (74, 14.57%), ether lipid metabolism (68, 13.39%), sphingolipid metabolism (54, 10.63%), linoleic acid metabolism (44, 8.66%), steroid biosynthesis (42, 8.27%), arachidonic acid metabolism (21, 4.13%), steroid hormone biosynthesis (16, 3.15%), synthesis and degradation of ketone bodies (12, 2.36%), and primary bile acid biosynthesis (5, 0.98%) pathways were involved in these pathways (Fig. 7).

Analysis of DEGs at four developmental stages

Analysis of the fatty acid compositions of C. majus seeds and elaiosomes showed that the elaiosome-bearing seeds predominantly contained the monounsaturated fatty acid, linoleic acid (C18:2) (77.9%), followed by the monounsaturated fatty acid, oleic acid (C18:1) (16%) (Reifenrath, Becker & Poethke, 2012). Furthermore, saturated fatty acids accounted for only 6% of the composition, and most of the saturated fatty acids were identified as palmitic acid (C18:1) (3.5%). However, in only elaiosomes, the most abundant fatty acid also was linoleic acid (37.2%), followed by palmitic acid (34.8%), and oleic acid (24.5%). Previous developmental and histochemical observations (Yang et al., 2015) indicated that the elaiosome meristem of C. majus was formed during the flower-bud stage (Ch01), while elaiosome cell growth and enlargement primarily occurred from the flowering stage (Ch02) to the young seed stage (Ch03); the maturation of the seed and elaiosome occurred from the young seed stage to the mature seed stage (Ch04). Therefore, to obtain insight into the regulation of the biosynthesis of saturated fatty acids and unsaturated fatty acids in elaiosomes and seeds, we analyzed the DEGs of the lipid metabolism pathways in Ch02 compared to Ch01, Ch03 compared to Ch02, and Ch04 compared to Ch03.

Identification of unigenes related to fatty acid biosynthesis

From the fatty acids biosynthesis pathway (Path: ko00061) and the biosynthesis pathway of unsaturated fatty acids (Path: ko01040) at the four developmental stages of C. majus elaisomes and seeds (File S1), we identified 41 significantly DEGs. These DEGs are shown in File S2 and the corresponding genetic heat map is shown in Fig. 8. Four unigenes encoding the enzymes malonyl CoA-ACP transacylase (fabG) (Toomey & Wakil, 1966), acyl-CoA oxidase (ACOX1) (Oaxaca-Castillo et al., 2007), acyl carrier protein desaturase (DESA1) (Kachroo et al., 2007), and β-ketoacyl–acyl carrier protein synthase II (fabF) (Magnuson, Carey & Cronan, 1995) were found to be upregulated in Ch02 compared to Ch01. By contrast, 24 unigenes encoding the enzymes acetyl-coA carboxylase (Accase) (Bilder et al., 2006), acyl coenzyme A synthetase (ACSL) (Black et al., 1992), DESA1, acyl-carrier-protein S-malonyltransferase (fabD) (Magnuson et al., 1992), fabF, fabG, enoyl-[acyl-carrier-protein] reductase (fabI) (Bergler et al., 1996), β-hydroxyacyl-acyl carrier protein dehydratases (fabZ) (Heath & Rock, 1996), acyl-lipid omega-6 desaturase (Delta-12 desaturase) (FAD2) (Kainou et al., 2006), acyl-ACP thioesterases (FatA/B) (Salas & Ohlrogge, 2002), and stearoyl-CoA desaturase (Delta-9 desaturase) (SCD) (Stukey, McDonough & Martin, 1990) were clearly upregulated in Ch03 compared to Ch02. However, in Ch04 compared to Ch03 identified that out of the 24 previously mentioned unigenes, 13 genes encoding Accase, ACSL, fabD, fabF, fabG, fabI, and SCD were downregulated, while the 11 remaining genes were not differentially expressed. The above results indicated that there are more relevant DEGs in Ch02 vs Ch03 compared to Ch01 vs Ch02 and Ch03 vs Ch04. Based on the observations of C. majus elaiosome development (Yang et al., 2015), the remarkable cell growth of the elaiosome also occurred during the flowering stage to the young seed stage (Ch02 to Ch03). Therefore, we inferred that the biosynthesis of saturated fatty acids and unsaturated fatty acids in C. majus elaisomes might primarily occur during this period.

Identification of unigenes related to the biosynthesis of saturated fatty acid

In the fatty acid biosynthesis metabolic pathway, a comparison of the DEGs in Ch02 compared to Ch01 (Table 2; File S3A) showed that most genes that are primarily involved in fatty acid biosynthesis were not significantly expressed during the flowering stage (Ch02), and only a few genes were upregulated. By contrast, the gene expression analysis in Ch03 compared to Ch02 (Table 2; File S3B) revealed that the majority of genes were significantly expressed during the young seed stage (Ch03). Most of these expressed genes, which encode Accase, ACSL, DESA1, fabD, fabF, fabG, fabI, fabZ, and FatA/B, were significantly upregulated in Ch03. On the other hand, a gene expression analysis in Ch04 compared to Ch03 (Table 2; File S3C) revealed that almost all genes were significantly downregulated during the mature seed stage (Ch04). The above results indicated that the biosynthesis of saturated fatty acid in C. majus elaisomes and seeds also mainly occurred during the flowering stage to the young seed stage. However, some unigenes that encode FabG, exhibited inverse expression patterns. The above-mentioned genes were found to be downregulated in Ch03 compared to Ch02 but were upregulated in Ch02 compared to Ch01 and in Ch04 compared to Ch03 (Table 2). These genes are potentially involved in fatty acid synthesis through negative feedback signaling.

The enzymes FatA/B and DESA1, which are the key enzymes involved in palmitic acid synthesis (see File S3), are particularly noteworthy, because palmitic acid is the most abundant saturated fatty acid in C. majus elaisomes but is present in very low quantities in the corresponding seeds. FatA/B was slightly upregulated in Ch02 compared to Ch01 and significantly upregulated in Ch03 compared to Ch02, but was significantly downregulated in Ch04 compared to Ch03. DESA1 was significantly upregulated in both Ch02 compared to Ch01 and Ch03 compared to Ch02, but significantly downregulated in Ch04 compared to Ch03. The above findings indicate that the biosynthesis of palmitic acid in C. majus elaiosomes occurs from the flower bud stage to the flowering stage (Ch01 to Ch02) but primarily occurs from the flowering stage to the young seed stage. Hughes, Westoby & Jurado (1994) compared the fatty acid compositions of 12 plant species of elaiosomes and seven orders of insects and showed that palmitic acid levels in elaiosomes were particularly similar to those in insects. However, because C. majus was not included in Hughes’s study, whether the high palmitic acid content in C. majus elaiosomes is similar to that in the ants that disperse the C. majus seeds remains unknown and needs further study.

Identification of unigenes related to the biosynthesis of unsaturated fatty acid

Figures 8 and 9 show that significantly upregulated genes were identified in Ch03 compared to Ch02 and were involved in three processes in the biosynthesis of unsaturated fatty acid metabolic pathway. The first process is the conversion of palmitic acid (C16:0) to stearic acid (C18:0), which involves genes encoding the enzyme, fabG. The second process is the conversion of stearic acid (C18:0) to oleic acid (C18:1, Δ9), which involves genes encoding the enzymes, SCD (Delta-9 desaturase) and DESA1. The third process is the conversion of oleic acid (C18:1, Δ9) to linoleic acid (C18:2, Δ9Δ12), which involves genes encoding the enzyme, FAD2 (Delta-12 desaturase). The above findings showed that some palmitic acid was converted into stearic acid (saturated fatty acid), and the latter was subsequently desaturated to form oleic acid and linoleic acid.

Linoleic acid is a poly-unsaturated fatty acid that cannot be synthesized by insects de novo and is thus required to be incorporated in their diet (Dadd, 1973). Previous studies have suggested increasing linoleic acid levels by stimulating the diaspore-carrying behavior of the ants (Gammans, Bullock & Schönrogge, 2005). The high linoleic acid content in both elaisomes and seeds of C. majus indicated that the whole diaspores (seed + elaiosome) were probably used to attract ants that contain high levels of (essential) nutrients. In addition, oleic acid was shown to promote ant removal behavior (Marshall, Beattie & Bollenbacher, 1979; Skidmore & Heithaus, 1988; Brew, O’Dowd & Rae, 1989; Hughes, Westoby & Jurado, 1994). The higher oleic acid content in C. majus elaiosomes and the lower oleic acid content in their corresponding seeds suggested that oleic acid in C. majus might act as the primary stimulating substance to induce the collection of its diaspores by ants, as previously reported in other plant species.

Validation of unigenes by quantitative real-time PCR

To validate the reliability of the RNA-seq results, four unigenes encoding ACSL, DESA1, fabG, and FatB were randomly selected for qRT-PCR validation at the four different developmental stages. These four genes were all related to the fatty acid biosynthesis metabolic pathway and biosynthesis of unsaturated fatty acid metabolic pathway. The qRT-PCR results (File S4) showed that the expression patterns of the four selected genes were very similar to those of RNA sequencing data in Ch03 compared to Ch02 and in Ch04 compared to Ch03 (Fig. 10). In Ch01 compared to Ch02, the expression pattern of DESA1 by qRT-PCR was consistent with that from RNA sequencing (Fig. 10B), and ACSL, fabG, and FatB showed slight differences (Figs. 10A, 10C and 10D). It was acceptable that some differences in direct comparison between qRT-PCR and RNA-seq results would occur due to biases in the library preparation for RNA-seq, different normalization approaches, and other technical biases (Li et al., 2010; Wei, Chung & Zhao, 2011). qRT-PCR results of the selected genes were generally consistent with the RNA sequencing results, indicating that the RNA-seq results were reliable.