Comparative transcriptome analysis on candidate genes involved in lipid biosynthesis of developing kernels for three walnut cultivars in Xinjiang

Wenqing Wng, Ho Wen, Qing Jin, Wenjun Yu, Gen Li, Minyu Wu,Hongjin Bi, Lirong Shen,*, Cuiyun Wu,*

a Department of Food Science and Nutrition, Zhejiang Key Laboratory for Agro-Food Processing, Zhejiang University, Hangzhou 310058, China

b College of Horticulture and Forestry, Tarim University, Alar 843300, China

Keywords:

Walnut

Oil accumulation

Triacylglycerol assembly

Transcription factors

Correlation analysis

A B S T R A C T

Walnut (Juglans regia L.) is a good source of lipids and polyunsaturated fatty acids (PUFAs). In order to explore the biosynthesis molecular mechanism of oil accumulation and fatty acid (FA) synthesis in walnut,the samples at different development periods of three walnut cultivars, ‘Zhipi’ (ZP), ‘Xinwu 417’ (W417)and ‘Xinwen 81’ (W81) were collected for transcriptomic analysis. The analysis of oil accumulation and FA profiles showed that the oil content in mature walnut kernel was nearly 70%, and over 90% of FAs were PUFAs. We identified 126 candidate genes including 64 genes for FA de novo synthesis, 45 genes for triacylglycerol assembly, and 17 genes for oil bodies involved in lipid biosynthesis by RNA-sequencing. Ten key enzymes including ACCase, LACS6, LACS8, SAD, FAD2, FAD3, LPAAT1, DGAT2, PDAT2, and PLC encoded by 19 genes were highly associated with lipid biosynthesis. Quantitative PCR analysis further validated 9 important genes, and the results were well consistent with our transcriptomic data. Finally, 5 important transcription factors including WRI1, ABI3, FUS3, PKL and VAL1 were identified, and their main regulatory genes might contain ACCase, KASII, LACS, FAD3 and LPAAT which were determined through correlation analysis of expression levels for 27 walnut samples. These findings will provide a comprehensive understanding and valuable information on the genetic engineering and molecular breeding in walnut.

1. Introduction

Walnut (Juglans regiaL.) is an important nut tree which belongs to the Juglandaceae family. It is widely cultivated in North America,Asia and Europe, especially in China with annual yield accounting for 48% of world production [1]. As one of the most important sources of edible nuts, walnut is rich in proteins, amino acids, lipids,carbohydrates, various trace elements and minerals [2,3]. Particularly,the oil content in walnut kernels is as high as 65% and the polyunsaturated fatty acids (PUFAs) in total fatty acids (FAs) reach up to 90%. Notably, the proportion of the two essential FAs, linoleic acid (LA) andα-linolenic acid (ALA) in walnut is approximately 4:1–6:1 which is considered to be effective in preventing cardiovascular diseases [4], and benefit for patients with type II diabetes [5].

In oil seeds, lipid biosynthesis includes three continuous steps.The initial step involves FA synthesis in plastids, which is mainly catalyzed by acetyl-CoA carboxylase (ACCase), 3-ketoacylacyl-carrier-protein synthase (KAS), 3-ketoacyl-ACP reductase(KAR), 3-hydroxyacyl-ACP dehydratase (HAD), enoyl-ACP reductase (ENR), fatty acyl-ACP thioesterase A/B (FATA/B) [6], and long-chain acyl-CoA synthetase (LACS). The produced acyl-CoAs are then transported to endoplasmic reticulum (ER) for triacylglycerol (TAG)assembly by the catalysis of glycerol-3-phosphate acyltransferase(GPAT), lysophosphatidic acid acyltransferase (LPAAT),phosphatidic acid phosphatase (PAP) and diacylglycerol acyltransferase (DGAT). In addition, phosphatidylcholine (PC)provides a more complex compensation pathway for TAG assembly due to the central intermediate in the flux of FAs [7]. Finally, TAG is stored as oil bodies through the combination with several proteins (i.e.,oleosin, caleosin and steroleosin) [8].

Increasing oil yield has been a vital target for many researchers owing to constant growing demand for oil consumption which is predicted to double by 2030 [7,9]. The metabolic engineering strategies to increase TAG content could be mainly categorized into three aspects including increasing FA supply, TAG assembly activities and decreasing TAG degradation pathways [10-13].Therefore, our knowledge of molecular mechanisms of lipid biosynthesis in seed, especially identification of genes involved in potential rate-determining steps of TAG synthesis is crucial for molecular plant breeding. In the last few years, the in-depth study on FA metabolism inArabidopsisas a model plant provided a theoretical basis for other researches [14]. Numerous candidate genes related to lipid metabolism has been identified in some oil plants by transcriptome sequencing such as oil palm [15,16], jatropha [17],Brassica napus[18,19], hickory [20,21], soybean [22]and tree peony [23,24], which expanded our understanding of molecular mechanisms of seed oil accumulation.

To date, a few studies have been conducted at transcription levels to investigate lipid biosynthesis in walnut seed and non-seed tissues.Dozens of genes such as ACCase, LACS and the FAD gene family related to FA synthesis have been determined [25,26]. Similarly,a comparative study by RNA sequencing among roots, stems and leaves of walnut was performed and 16 genes significantly correlated with FA content in different walnut organs were identified [27].However, the genes that play a crucial role on TAG assembly and the transcriptional regulatory genes related to lipid biosynthesis in walnut were still not identified in these studies. In addition, different cultivars vary in oil content and FA composition, and the expression pattern of key genes may be different. Therefore, the molecular mechanism and the gene regulatory network of lipid biosynthesis remain still unclear and need to be further explored. In this study, the comparative transcriptome analysis of walnut kernels from three cultivars over walnut embryo development to identify the key genes and regulatory mechanism associated with oil accumulation were implemented.The results will provide a better understanding and new insight into molecular mechanism of lipid biosynthesis as well as engineering strategies of enhancing oil accumulation, especially proportion of PUFAs in walnut for plant breeders and metabolic engineers.

2. Materials and methods

2.1 Plant materials

Fruits were collected from three 25-year-old walnut cultivars,namely, ‘Zhipi’ (ZP), ‘Xinwu 417’ (W417) and ‘Xinwen 81’ (W81)planted in Wensu Woody Grain and Oil Forest Farm (41°N, 80°E),Xinjiang, China. Walnut samples were taken from 77 days after flowering DAF at early cotyledon stage, and then every 7 days (the sampling time was 98 DAF after 84 DAF) until 147 DAF at maturity,resulting in a total of 10 sampling times (77, 84, 98, 105, 112, 119,126, 133, 140, 147 DAF) (Table S1). For each sampling time, we sampled 5 walnut fruits for each of 5 trees in different directions(i.e., north, south, center, west, and east) for a total of 25 fruits each cultivar. The walnut kernels were picked off, a part of which were stored in a -80 °C freezer after quick freezing in liquid nitrogen for transcriptome sequencing, and the remained samples were directly stored at -40 °C for lipid analysis.

2.2 Determination of lipids and fatty acid content

The walnut kernels were dried to a constant weight at 80 °C, and then total lipids were extracted in a Soxhlet apparatus using petroleum ether as solvent at 60 °C for 6 h. FA methyl esters (FAME) were analyzed based on ISO 5509 method using a GC-2014 gas chromatography(Shimadzu, Kyoto, Japan). According to the retention time and the peak area normalization, the relative content of FAs was determined. The lipid analysis was performed under a completely random experimental design,with 3 biological replicates for each sample.

2.3 RNA extraction, cDNA library construction and sequencing

Samples collected at T1 (84 DAF), T2 (98 DAF) and T3 (119 DAF)were chosen for transcriptome sequencing according to the morphological characteristics and dynamic changes in total oil content during embryo development of three walnut cultivars.Total RNA was extracted using TriZol Reagent according the manufacturer’s instructions (Invitrogen, USA) and genomic DNA was removed by DNase I (TaKaRa, Japan) digestion. RNA quality was determined by Agilent 2100 Bioanalyzer and quantified using the ND-2000 (Nanodrop Technologies, USA), and the RNA integrity number (RIN) was examined using Agilent 2100 Nanodrop. The RIN of all the samples was higher than 6.8, OD260nm/280nmand OD260nm/230nmwere above 2.00 and 1.50, respectively, indicating all the RNA samples could be used to construct sequencing library. The transcriptome library was prepared following TruSeq RNA Sample Preparation Kit from Illumina (San Diego, CA, USA) using 5 μg of total RNA, which was then sequenced with the Illumina HiSeq 4 000 (2 × 150 bp read length) supplied by Shanghai Majorbio Bio-pharm Technology Co.,Ltd. (China).

2.4 Transcriptome data analysis and functional annotation of genes

The raw paired end reads were trimmed and quality was controlled by SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. Then clean reads were separately aligned to the reference genome ofJuglans regia(https://www.ncbi.nlm.nih.gov/genome/?term=txid51240 [orgn]) with orientation mode using TopHat (http://tophat.cbcb.umd.edu/) [28].The detected genes were annotated against the following six protein databases, including EggNOG (http://eggnogdb.embl.de/#/app/home),Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg), Gene Ontology (GO, http://www.geneontology.org), Nonredundant (NR, http://www.ncbi.nlm.nih.gov), Swiss-Prot(http://www.expasy.ch/sprot), and Pfam (http://pfam.xfam.org) using BLAST software (v2.9.0).

2.5 Analysis of differentially expressed genes

The gene expression levels were quantified as transcripts per million reads (TPM) by RSEM (http://deweylab.github.io/RSEM/) [29].DESeq2R package (v1.24.0) was used to identify differential expression genes (DEGs) defined asPadj< 0.05 and |log2(Fold change)| ≥ 1 between two different samples based on TPM values [30]. Clustering analysis of DEGs was performed by R statistical package, and KEGG pathway enrichment analysis of DEGs was carried out at BH(Benjamini and Hochberg)-correctedP≤ 0.05 compared with the whole-transcriptome background by KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [31]. The correlation analysis of gene expression based on TPM values over walnut embryo development in the three walnut cultivars was measured by Person’s correlation coefficient,and the significance was calculated usingPcorrected by BH.

2.6 Quantitative RT-PCR analysis of lipid-related genes

Total RNAs of walnut kernels from three developmental stages (T1, T2 and T3) in the 3 walnut cultivars were reverse-transcribed to synthesize cDNAs using MonScript RT All-in-One Mix (Monad,Zhuhai, China). Primers (Table S2) for 9 target genes involved in lipid biosynthesis and β-actin which was chosen as the internal reference gene for normalization were designed on Primer-Blast (https://blast.ncbi.nlm. nih.gov/Blast.cgi). qRT-PCR was performed on ViiA 7 (Applied Biosystems) using MonAmp SYBR Green qPCR Mix with Low ROX (Monad) in three biological replicates and three technical repetitions. The relative expression level was calculated with 2–ΔΔCtmethod [32].

3. Results

3.1 Oil content and fatty acid composition in developing kernels of three walnut cultivars

Three walnut cultivars, ZP, W417 and W81 showed an almost consistent change trend of morphological characteristics as well as dynamic oil accumulation which was identified as a “fast followed by slow” type (Fig. 1). To be more specific, during the early developmental stage at 77 DAF) the oil content was relatively low (i.e., ～4%)in accordance with gelatinous and watery embryos. At 84 DAF, the walnut embryos began to solidify and were filled with white solid cotyledon. Simultaneously the oil content increased dramatically and reached up to 35.19%, 20.36% and 29.33% in ZP, W417 and W81,respectively. This result indicated that the oil content varied among the different cultivars at 84 DAF. Subsequently, it remained high oil accumulation rate at 98 DAF and slowed down a little at later stage until 119 DAF from which the oil accumulation nearly stopped with slight separation between walnut kernel and cotyledon. Finally, the walnuts were completely mature with green peels off, and the oil content of the three walnut cultivars was almost the same and up to ～69%.

Fig. 1 Walnut profiles and total oil content in developing walnut kernels. (A) Morphological characteristics of walnut profiles during embryo development at 84 DAF, 98 DAF, and 119 DAF for 3 walnut cultivars. (B) Dynamic changes of total oil content during embryo development for 3 walnut cultivars (mean ± SD, n = 3).

Then, the FA composition of walnut oil of different walnut cultivars was analyzed. There were 6 dominant FAs in walnut, 2 kinds of saturated FAs (i.e., palmitic acid (16:0), and stearic acid(18:0)) and 4 kinds of unsaturated FAs (oleic acid (18:1), LA(18:2), ALA (18:3), and eicosenoic acid (cis-11-eicosenoic acid,20:1)) (Fig. 2). The 3 walnut cultivars revealed similar FA composition and change trends. The relative content of LA was always the highest at the whole developmental stages, and dropped from ～70% at 84 DAF to ～50% at 119 DAF except W81, and subsequently it increased to ～60% at 126 DAF, after which it remained almost the same (56%). Compared to LA, the content of oleic acid in ZP and W417 showed opposite trends, increasing from about 6% to over 30% and then performing as small fluctuation. The proportions of ALA and palmitic acid changed a little and fluctuated by～1%. The remaining two FAs accounted for below 2%. In general, the composition of PUFAs reached up 90%. Those results demonstrated that the oil content increased rapidly at 84 DAF (T1) and 98 DAF (T2),and remained relatively stable at 119 DAF (T3). Therefore, the kernel samples from these three developmental stages for the 3 walnut cultivars were used for subsequent transcriptome analysis.

Fig. 2 Fatty acid composition during embryo development for 3 walnut cultivars (mean ± SD, n = 3). (A) ZP, (B) W417, (C) W81.

3.2 Transcriptome sequencing and alignment to reference genome

Transcriptome sequencing of walnut kernels involved in different developmental stages and cultivars was performed to explore the potential molecular mechanisms in lipid biosynthesis. A total number of 67.66, 72.21 and 70.64 Gb raw bases were obtained in W417, W81 and ZP, respectively. After removing adaptor sequences, low-quality bases and undetermined bases, over 97% high-quality clean bases were produced. The error rate was below 0.03, and the Q20 and Q30 value were above 97.70% and 93.34%, respectively, indicating the sequencing data was of high quality (Table 1).

Table 1Summary of the sequencing data of walnuts.

All clean reads were aligned to the reference genome by Hisat2 (v2.1.0) software to acquire mapped data for subsequent transcript assembly, expression calculation, etc. As a result, over 95% of the clean reads were mapped (Table 1), indicating the selected reference genome was quite appropriate and there was no pollution in the experiments. Furthermore, mapped reads in each sample were spliced and finally merged together using Stringtie (v1.3.3b).100 907 transcripts were produced and the length distribution showed length between 400 bp and 1 800 bp was evenly distributed and accounted for 7%-8%. Notably, 33 738 (33%)transcripts were longer than 1 800 bp, indicating long sequences dominated the transcripts of walnut (Fig. 3).

Fig. 3 Length distribution of walnut transcripts.

3.3 Functional annotation of expressed genes

All the 34 112 expressed genes with non-zero expression levels in at least one sample were annotated against 6 public databases,NR, Swiss-Prot, Pfam, COG, GO and KEGG. A total number of 30 320 (88.88%) expressed genes had the most significant BLAST matches with known proteins in the 6 databases and 8 538 genes were found to be annotated in all the 6 databases (Table 2 and Fig. S1). There were the most annotated genes accounting for 87.54% (29 862) and 88.86% (30 312) in COG and NR among the databases, respectively, whereas GO and KEGG possessed the lowest accounted rate, 44.79% (15 279) and 38.09% (12 992),respectively. To some extent, it’s very necessary to choose a reference genome for transcriptome analysis of walnuts due to low accounted rate in GO and KEGG.

Table 2Functional annotation of walnut genes and transcripts in 6 databases.

3.4 Analysis of DEGs in developing kernels of 3 walnut cultivars

Principal component analysis (PCA) of the 27 walnut samples (3 typical times each cultivar, 3 duplicates each typical time) showed good clustering among all the biological replicates, especially at T2 and T3 in the three walnut cultivars, suggesting that the samples were collected in the right developmental stages. PC1 contributed to the most variance (53.32%) in the datasets and distinguished samples at T1 from T2 and T3. However, PC2 captured much lower variance (13.77%), accounting for the differences among samples from T2 and T3, and the difference between ZPT1 and W417T1 presented by PCA analysis might be associated with the major difference in oil content between ZP and W417 at T1 stage mentioned before (Fig. 4A).

Fig. 4 Principal component analysis and Venn diagram representation of DEGs identified in 3 walnut cultivars. (A) Principal component analysis of 27 walnut samples. (B) Number of intersections of DEGS between different stages in W417, W81, and ZP. (C) Number of unions of DEGs between W417, W81, and ZP at each of the 3 developmental stages. (D) Combined number of DEGs from union of B and intersection of C.

In order to find out the dynamic changes of gene expression over the kernel development, analysis of differentially expressed genes (DEGs) was performed using transcripts per million reads (TPM)data from 27 samples. Three pairwise comparison groups including samples from different developmental stages and cultivars were analyzed. 1 800, 1 145, and 1 352 DEGs were contained in all three comparison groups between developmental stages (i.e. T1vs. T2,T2vs. T3, and T1vs. T3) within W417, W81, and ZP, respectively, and this analysis yielded a total of 3 375 non-redundant DEGs (Fig. 4B).A different statistical approach by taking the union of DEGs first and intersection next was adopted when it came to the identification of DEGs between cultivars due to only a few differentially co-expressed genes found in comparisons between 3 cultivars at three stages.There were 12 094, 8 133, and 7 402 non-redundant DEGs existing in ZP vs. W417, ZP vs. W81, and W417 vs. W81 at 3 stages, respectively.In this way, 2 578 DEGs contained in all three comparisons were determined (Fig. 4C). Overall, 5 040 non-redundant DEGs redundant genes were identified, and 913 DEGs were uncovered both in comparisons from different stages and cultivars (Fig. 4D).

To better understand the expression pattern of the DEGs among the comparison groups, clustering analysis has been searched. A hierarchical cluster analysis of 5 040 DEGs co-expressed among different stages and different cultivars based on the TPM values was conducted. As a result, T1 and T2 stages of 3 cultivars showed good clustering (Fig. 5A).All the DEGs were clustered into 5 groups. Subcluster I consisted of 1 118 genes which represented a decrease in expression levels from T1 to T3 stage in all cultivars, especially nearly linear decrease in ZP and W417, indicating negative associations with oil content in walnut kernels. 235 genes in Subcluster II were down-regulated from T1 to T2 stage in W417 and W81 and up-regulated thereafter,which was different from the constant increase from T1 to T3 stage in ZP. In contrast, 412 genes in Subcluster III displayed a relatively big rise in gene expression from T1 to T2 stage and an abrupt decline soon afterwards. Subcluster IV contained 1 289 genes which exhibited increasing expression pattern from T1 to T2 stage and remained unchanged or a slight decline in gene expression thereafter. Similarly, 1 986 genes in Subcluster V were positively associated with the oil content due to the constant up-regulated expression levels over the fruit development. The clustering analysis suggested the genes in Subcluster IV and Subcluster V might play a crucial role in lipid biosynthesis and kernel development (Fig. 5B).

Fig. 5 Cluster analysis of identified DEGs in developing walnut kernels. (A) Hierarchical cluster analysis of DEGs, the color bar showed expression levels of DEGs standardized as lg TPM, the red represented DEGs with high expression levels and the blue represented DEGs with high expression levels. (B) The 5 subclusters of different expressional patterns. (B1-B5) Subcluster I (1 118 genes), subcluster II (235 genes), subcluster III (412 genes), subcluster IV (1 289 genes),subcluster V (1 986 genes), respectively.

3.5 GO annotation and KEGG enrichment analyses

To determine specific functions of these DEGs, GO analysis was carried out. 1 348 DEGs between cultivars (group 1) and 1 795 DEGs between developmental stages (group 2) were annotated in GO database, and both divided into three main functional categories including biological process (BP), cellular component (CC) and molecular function (MF) with 36 subcategories which contained more than 10 genes in group 1 or group 2 (Fig. 6). In detail, the BP category was assigned into 13 subcategories, in which “metabolic process” and“cellular process” were the two predominant subcategories including 681 (50.52%) and 584 (43.32%) DEGs, respectively in group 1,and 915 (50.97%) and 789 (43.96%) DEGs, respectively in group 2.The CC category was further classified into 12 subcategories,and the largest subcategory was both “cell part” which contained 453 (33.61%) and 636 (35.43%) DEGs, respectively in group 1 and group 2, followed by “membrane part” with 346 (25.67%) and 522(29.08%) DEGs. The MF category was mapped into 11 GO terms,and the most two abundant subcategories were “catalytic activity” and“binding”, which contained 732 (54.30%) and 586 (43.47%) DEGs in group 1, and 960 (53.48%) and 797 (44.40%) DEGs in group 2.These results suggested strong metabolic activities occurred during the growth and development of different walnut cultivars.

Fig. 6 GO classification of DEGs between cultivars and between stages in walnut.

Furthermore, 3 375 DEGs between developmental stages and 2 578 DEGs between cultivars were analyzed to explore the metabolic pathways involved by these DEGs through KEGG enrichment analysis. For DEGs between stages, the top enriched KEGG pathways were “starch and sucrose metabolism” (60), “protein processing in endoplasmic reticulum” (47), “amino sugar and nucleotide sugar metabolism” (40), “glycolysis/gluconeogenesis” (39).Notably, there were several highly enriched pathways involved in lipid metabolism, including “FA degradation”(19), “FA biosynthesis” (17), “steroid biosynthesis” (15),“biosynthesis of unsaturated fatty acids” (12) (Fig. 7A). For DEGs between cultivars, “protein processing in endoplasmic reticulum” (41)and “glycolysis/gluconeogenesis” (33) were two of the most enriched KEGG terms, followed by “pyruvate metabolism” (19)and circadian rhythm-plant (19), whereas only one pathway, “FA biosynthesis” (15), was enriched in lipid metabolism (Fig. 7B).

Fig. 7 KEGG enrichment analysis of the DEGs. (A) KEGG enrichment pathways of DEGs between developmental stages in 3 walnut cultivars. (B) KEGG enrichment pathways of DEGs between 3 walnut cultivars during kernel development.

3.6 Identification of genes related to lipid biosynthesis during kernel development in 3 walnut cultivars

DEGs involved in lipid biosynthesis were crucial to explore the molecular mechanisms leading to the increase of oil content in walnut.Therefore, we constructed the lipid biosynthesis pathway visualized by the schematic diagram and heatmap to display the dynamic changes of genes at expression levels during embryo development in three walnut cultivars (Fig. 8). Three classes of genes involved in FAde novosynthesis (64 genes) in the plastid, TAG assembly (45 genes)on the ER and oil bodies (17 genes) were uncovered (Table S3). In order to find out all the highly lipid-related genes which had similar expression pattern during developmental stage in ZP, W417 and W81,we assumed genes up-regulated or down-regulated simultaneously in at least two comparison groups as co-expressed DEGs. In this way, ACCase, LACS, FAD3, LPAAT, PAP, PLC involved in 11 co-expressed DEGs altogether were found up-regulated during the kernel development. In contrast, 17 co-expressed DEGs including KASIII, HAD, KASI, SAD, LACS, FAD2, GPAT, PLC, PLD, PDAT,PDCT and oleosin in total were found down-regulated. It should be noted that some enzymes such as SAD, LACS, PLC and PLD possessed several copies which exhibited the opposite expression patterns. Obviously, those up-regulated DEGs from T1 to T3 for the 3 cultivars with high expression levels were more relevant to the oil accumulation pattern. In addition, the genes which showed much higher expression levels at T2 than T1 should also be noticed due to the rapid oil accumulation from T1 to T2.

Fig. 8 Transcriptional visualization of differentially expressed lipid-related genes between developmental stages in three cultivars. Names of the genes were exhibited on the right of the heatmaps and the key enzymes were marked in red. The color of the heatmaps represents log2(fold change) in which the red showed up-regulated genes and the blue showed up-regulated genes. Abbreviations: MCAAT, malonyl- CoA-ACP transacylase; FAD, fatty acid desaturase; LPCAT,1-acylglycerol -3-phosphocholine acyltransferase; PDCT, phosphatidylcholine: diacylglycerol choline phosphotransferase; PLA2, phospholipase A2; PLC,phospholipase C; PLD, phospholipase D; PDAT, phospholipid: diacylglycerol acyltransferase; SAD, stearoyl-ACP desaturase.

3.6.1 DEGs related to fatty acid synthesis

Plastidial ACCase which catalyzes the first committed step in FA synthesis is a multi-subunit complex comprising four proteins: biotin carboxylase (BC), biotin carboxyl carrier protein (BCCP), carboxyl transferase subunit alpha (α-CT) and carboxyl transferase subunit beta (β-CT). 15 genes in total encoding the four subunits except forβ-CT were detected (Fig. 8 and Table S3). It was noteworthy that two BCCP copies (Jr32143andJr35731) and oneα-CT copy (Jr4589)were relatively abundant and showed much higher expression levels at T2 than T1 in three cultivars which may be associated with the rapid oil accumulation from T1 to T2 (Fig. 9A and 1B). Interestingly,Jr25895(BCCP) andJr20278(α-CT) were quite inactive at T1 but were activated at T2 in ZP and W417 and at T3 in W81 (Fig. 9A),suggesting coordinated expression of the two subunits dominated the first rate-limiting step in FA synthesis.

Fig. 9 Temporal changes at transcriptional levels based on TPM values of (A) ACCase, (B) LACS, (C) SAD, (D) FAD2, (E) FAD3, (F) LPAAT1, DGAT2,PDAT2, (G) PLC, and (H) caleosin, steroleosin in the 3 walnut cultivars.

Fig. 9 (Continued)

Eight members of LACS family consisting of 11 genes named LACS1, LACS4, LACS9 possessing two copies for each and LAC, LAC2,LACS6, LACS7, LACS8 with one copy for each in total were obtained in walnut (Table S3). Only LACS6 (i.e.,Jr9141) and LACS8 (i.e.,Jr14571)were relatively abundant and exhibited a major upward trend at expression levels from T1 to T3 in the three cultivars (Fig. 9B), emphasizing the significance of these two genes in FA transport to ER.

3.6.2 DEGs related to synthesis of unsaturated fatty acids

The silence construct of SAD which converts stearic acid into oleic acid was proved to double the stearic acid content of triacylglycerol inChlamydomonas reinhardtii[33], highlighting the essential role in oleic acid content. Among 4 SAD copies uncovered in our study, it was noteworthy thatJr14834were up-regulated from T2 to T3, and exhibited an abrupt increase at expression levels with 16-fold change in the three cultivars, which was strongly associated with the rapid rise in oleic acid content at later stage. In addition,Jr29135andJr8426were more abundantly expressed at T2 and T3 than T1 in ZP, whereas no obvious changes of expression levels were observed in W417 and W81 (Fig. 9C and Table S3), which may explain the more significant increase in oleic acid content from 6.02% to 14.34% in ZP than W417 and W81 with increase by 2%-4% (Fig. 2).The results indicated the coordinating expression pattern of the three SAD copies in oleic acid accumulation in walnut.

Enzymes associated with synthesis of unsaturated FAs also include FAD2 and FAD3, FAD7/FAD8, which can further desaturate the 18:1-PC in PC pool converted from the 18:1-ACP generated by SAD in plastid into 18:2-PC and 18:3-PC, respectively. FAD2 was encoded by two genes, in whichJr14671was down-regulated both in T1vsT2 and T1vsT3 with over 8-fold changes of expression levels in the three cultivars (Fig. 9D and Table S3), indicating the close connection with the decline by about 20% in LA content from T1 to T3 (Fig. 2). Furthermore, we detected 3 genes encoding FAD3,two of which displayed quite low expression levels at the whole developmental stages. Strikingly,Jr595was transcribed remarkably with TPM values from 1 987 to 9 765 at T2 and T3, and up-regulated both in T1vsT2 and T1vsT3 in the three cultivars (Fig. 9E and Table S3), denoting the dominant role in ALA synthesis in walnut.Interestingly, three chloroplasticω-3 fatty acid desaturase genes including FAD7/FAD8 kept low expression levels all along.

3.6.3 DEGs related to TAG assembly and oil bodies

The conventional Kennedy pathway is considered to play a great part in the de novo assembly of TAG in most organisms [9]. There are 3 important acyltransferases including GPAT, LPAAT and DGAT and one PA phosphatase (PAP), in which GPAT and PAP encoded by 5 genes and 3 genes, respectively, both exhibited low expression levels (Table S3). Five genes have been annotated for LPAAT in our study including 3 LPAAT1 and 2 LPAAT2 copies, in which only 1 copy (i.e.,Jr28986) encoding LPAAT1 was abundantly expressed and showed much higher expression levels at T2 and T3 than T1.We detected 2 genes and 1 gene encoding DGAT1 and DGAT2,respectively, and DGAT2 was more highly expressed than DGAT1,emphasizing the importance of DGAT2 for TAG synthesis. PDAT which utilized PC and diacylglycerol (DAG) as substrate for TAG synthesis is considered to partially compensate for TAG assembly and made a major contribution to oil accumulation in absence of DGAT [34,35]. PDAT2 encoded by one gene (i.e.,Jr4463) was more likely to lead a decisive role in compensation for DGAT on TAG assembly with much higher expression levels at T2 than T1, instead of 3 PDAT1 copies which were all lowly expressed during embryo development (Fig. 9F and Table S3).

PLC and PLD can hydrolyze the PC to PA which is then converted to DAG by PAP. We examined nine genes encoding PLC,most of which were lowly expressed. Notably, three genes includingJr12203,Jr20089andJr38162were relatively abundant and more active with over two-fold higher expression levels at T1 than T2 (Fig. 9G and Table S3), emphasizing the important compensatory effect at early stage on TAG assembly. This might explain why the quite low expression levels of GPAT and PAP in Kennedy pathway only had slight impact on oil accumulation in walnut. Almost all the eleven PLD copies were lowly expressed (Table S3), implying PLC may dominate the alternative pathway of DAG production.

Oil bodies is a common form of oil storage in mature seeds due to their stability provided by steric hindrance and electronegative repulsion of proteins [8]. There are three types of oil bodies, namely,oleosin, caleosin and steroleosin encoded by 9, 3 and 5 genes,respectively. The most abundant proteins in oil bodies are oleosins,in accordance with our discovery that almost all the oleosin copies displayed extremely high expression levels (Table S3). Furthermore,they had higher TPM values at T2 and T3 than T1 in ZP and W417.However, there was only one caleosin (i.e.,Jr18750) and one steroleosin (i.e.,Jr28664) copy which were transcribed much more highly than the other copies and showed higher abundance at T2 and T3 than T1 in agreement with the oil accumulation trends (Fig. 9H and Table S3).

3.7 Identification of transcription factors involved in lipid biosynthesis

The carbon metabolism in seeds is strictly regulated by a series of transcription factors such as WRI1, ABI3, FUS3, PKL and VAL1,which contribute to oil accumulation, embryo development and seed maturation [13]. We detected 3 genes (i.e.,Jr12895,Jr23096andJr28183) for WRI1, in which onlyJr23096was almost silenced all along, whileJr12895andJr28183were gradually activated during the embryo development, especially from T1 to T2 in ZP and W81,coinciding with the rapid increase in oil content at the early periods.Similar to the expression pattern of WRI1s, ABI3s and VAL1s encoded by 2 genes and PKL encoded by 1 gene displayed higher expression levels at T2 and T3 than T1, which was contrary to FUS3 with down-regulated expression pattern both in T2vsT1 and T3vsT1 in W417 and W81 (Table S4).

To further investigate the correlations between several TFs mentioned above and genes involved in lipid biosynthesis at transcriptional levels, the correlation coefficients calculated by Person’s correlation algorithm was visualized by heatmap (Fig. 10). It was notable that one WRI1 repeat,Jr28183,was positively and significantly correlated with 4 ACCase copies(i.e.,Jr40555,Jr35731,Jr8139andJr4589), all the 3 KASII copies and all the five LPAAT copies. In particular, there were strong positive correlations betweenJr28183andJr8139(BCCP),Jr11691(KAR),Jr34592(KASII),Jr11465(LACS),Jr595(FAD3),Jr13816andJr28986(LPAAT) at expression levels with correlation coefficient above 0.8, indicating the high positive regulation on these genes. Similarly, 2 ABI3 and 2 VAL1 copies were all positively correlated with 2 or 3 KASII copies, 1 FAD3 copy, 1 SAD copy, LACS8 and 4 LPAAT copies. Strikingly, only few genes displayed negative correlations, including ABIs, PKL,VAL1s with HAD, and PKL, VAL1 with SAD. Our data indicated the regulation of TFs was almost all positive and the main regulated targets contained ACCase, KASII, LACS, FAD3 and LPAAT.

Fig. 10 Correlation analysis between 5 transcription factors and lipid-related genes. The color of heatmap represented the correlation coefficient with positive correlation in red and negative correlation in blue. The dark green signified that the absolute value of correlation coefficient of these genes were less than 0.5 with Padj > 0.05.

3.8 Validation and analysis of candidate genes involved in lipid biosynthesis

Nine important genes, ACC, LACS6, LACS8, SAD, FAD2,FAD3, DGAT2, PDAT2 and PLC related to lipid biosynthesis were selected for validation based on relative expression level by qPCR to assess the accuracy of TPM values in RNA sequencing in ZP, W417 and W81, respectively (Figs. 11 and S2). The increasing relative expression level of all selected genes in three cultivars from T1 to T2 which were considered as important periods for oil accumulation was well consistent with our transcriptomic data.Especially,α-CT (Jr4589), LACS8 (Jr14571), FAD2 (Jr14671),FAD3 (Jr595), DGAT2 (Jr14467) and PLC (Jr12203) validated in ZP showed a high accordance in change trend between two results during three stages. Notably, the constant decline in expression level of FAD2 may be mostly responsible for the dramatic decrease from 73.22% (T1) to 49.73% (T3) in LA content of ZP, whereas extremely high expression level of FAD3 in three cultivars could explain why walnut was rich in ALA accounting for ～10% of FAs to a large extent. Our results indicated the significance of these genes in lipid biosynthesis.

Fig. 11 Expression patterns of nine candidate genes in ZP walnut kernels at 3 developmental stages (T1-T3). (A) ACC/α-CT (Jr4589), (B) LACS6 (Jr9141),(C) LACS8 (Jr14571), (D) SAD (Jr29135); (E) FAD2 (Jr14671), (F) FAD3 (Jr595), (G) DGAT2 (Jr14467), (H) PDAT2 (Jr4463), (I) PLC (Jr12203).

4. Discussion

Mature walnut kernel possesses high oil content up to ～70% and the PUFAs account for over 90% of FAs. In this work, we profiled oil content and FA composition of three widely grown walnut cultivars in Xinjiang (China), and investigated the underlying genetic basis.ZP, W417, and W81 presented almost the same developmental pattern and the kernels began to solidify at 84 DAF while the date in‘Lvling’ walnut cultivar was 50 DAF [26], which may be attributed to special topography and climate conditions in Xinjiang. We examined 6 dominant FAs in which only oleic acid and LA content changed dynamically accounting for a total proportion of nearly 80% in total FAs over kernel development, indicating a possible conversion from oleic acid to LA. These changing trends of FA profiles were similar to the storage component accumulation in other wood plants, such as oil palm [15], olive [36], and tung tree [37].

Genes related to oil accumulation were essential to investigate the molecular mechanism of lipid biosynthesis. As the first committed enzyme in FA synthesis, ACCase has received abroad attention in most studies [38]. It was reported that all the four subunits were coordinately expressed to promote FA synthesis [39], and lack ofβ-CT might limit the activity of ACCase [40]. In our study, 3 BCCP genes and 2α-CT genes were much more abundant at T2 than T1, implying that they might play an important role at early developmental stages in walnut. However, noβ-CT gene was obtained in the 3 walnut cultivars,which was similar to the discovery inCarya cathayensis[20]andC. illinoinensis[21], indicating a lack of importance ofβ-CT in nuts.LACS is generally presumed to convert FAs to acyl-CoAs which are transferred to ER for TAG assembly on the outer plastid envelope [7].Disruption or overexpression of LACS8 and double mutant of LACS8 and LACS9 hardly compromised FA content inArabidopsisseed, and LACS8 was not involved in acyl-CoA formation and FA transport from the plastid due to the location to the ER, unlike LAC9 located to the plastid [41]. LACS1, LACS7 and LACS9 in the “Lvling”walnut cultivar were supposed to matter in FA synthesis [26], while LACS6 and LACS8 were determined to be more associated with oil accumulation in this study, suggesting key genes involved in lipid biosynthesis may vary among walnut cultivars.

FAD2 and FAD3 located in ER catalyzed the formation of LA and ALA, respectively [42,43]. The overexpression of PrFAD2 and PrFAD3 inArabidopsisincreased LA and ALA content, and modulated the final proportion of the 2 FAs [23]. In our study,decreased FAD2 abundance was likely to result in the decline of LA proportion in PUFAs from T1 to T3 and high content of ALA may be attributed to high activity and expression of FAD3. However,not FAD3 but FAD7 genes were observed in pecan [21]and the“Lvling” walnut cultivar [26]. Interestingly, another study indicated a significant positive correlation of 0.991 between the expression ofJrFAD3and ALA content in walnut [25], and FAD3 contributed to ALA content in olive oil in the seed but FAD7 was mostly responsible for ALA content in the mesocarp [44]. These studies indicated FAD3 was more likely to be responsible for ALA accumulation.

DGATs, including DGAT1 and DGAT2, which acted as the rate limiting acyltransferase in determination of TAG assembly in Kennedy pathwayhave been recognized to play a crucial role on oil accumulation in plant tissues [45-47].ArabidopsisDGAT1 expressed in tobacco contributed to a surprising 20-fold increase in TAG content in tobacco leaves [48], whereas high expression of DGAT2 which showed ambiguous function inArabidopsiswas believed to take responsibility for oil and oleic acid accumulation inRicinus communis[49],Zea maysL. [50], andVerniciafordii[51]. It was DGAT2 that played an important role in TAG assembly in this study, revealing the similarity that DGAT2 rather than DGAT1 made a major contribution to TAG synthesis to pecan [21], and hickory [20]. Except for Kennedy pathway, alternative pathways or combinations of pathways were also utilized to produce TAG in different plants [9]. Dynamic changes of PLC at transcriptional levels were considered to play an important compensatory role in TAG assembly in the case of low expression of the key genes, such as GPAT and PAP in Kennedy pathway. Actually,a route alternate to Kennedy pathway was also observed to enhance TAG accumulation and ALA content inP. rockii[23].

Embryo development is subject to an array of transcription factors (TFs), such as WRI1, LEC1 and LEC2, PKL, FUS3 and ABI3. Mutants of WRI1 led to an 80% decrease of FA content in seed oil inArabidopsis, which may be explained by the dramatic reduction in carbon flux that provided pyruvate from sucrose [7,52],and overexpression of WRI1 remarkably increased seed oil content by 30% in maize [53]. However, WRI1 was considered to be lack of importance over pecan embryo development [21], differing from the high consistency between expression level of WRI1 and oil accumulation uncovered in this study. Moreover, the correlation analysis indicated WRI1 strongly associate with the expression of ACCase, KASII and LPAAT. It has been documented that WRI1 regulated the expression of at least 15 enzymes including ACCase,EAR and KASI involved in FA synthesis [54,55], and regulate positively in oil accumulation [56,57], which were confirmed in our study to some extent. While inTorreya grandiskernels, most of WRI1 transcripts were negatively correlated with the genes related to FA synthesis, suggesting that WRI1 may had opposite regulatory effect on these genes in various species. LEC2, FUS3 and ABI3 (AFLs) could not only directly interact with some genes involved in lipid synthesis, but also participate in oil accumulation in seeds in a more complicated way [58,59]. Besides, AFLs can activate the maturation program with combined action of LEC1 [13]. LEC1 and LEC2 may regulate the expression of WRI1, ABI3 and FUS3 in turn [56,60], whereas these two TFs were not detected in our study,indicating LEC1 and LEC2 might not associate with oil accumulation and embryo development in walnut. Considering that there are few studies on the regulation of these TFs in lipid biosynthesis, more research about molecular biology experiment such as gene knockout will enhance our knowledge on their interaction in enzyme activity.The results revealed the potential molecular mechanism of oil accumulation and lays the foundation for metabolic engineering and molecular breeding of walnut to enhance oil accumulation and FA composition.

5. Conclusion

In summary, 19 candidate genes encoding 10 key enzymes including ACCase, LACS6, LACS8, SAD, FAD2, FAD3, LPAAT1,DGAT2, PDAT2, and PLC were identified to be highly associated with oil accumulation. In particular, decreased FAD2 abundance is likely to result in the decline of LA proportion in PUFAs from T1 to T3 and high content of ALA may be attributed to high activity and expression of FAD3. Dynamic changes of PLC in expression level is considered to play an important compensatory role in TAG assembly in the case of low expression of the key genes in Kennedy pathway. The qPCR results have also well validated the transcriptome sequencing data. Furthermore, we explored the regulation of transcriptional factors including WRI1, ABI3, FUS3,PKL, and VAL1 on lipid-related genes through correlation analysis of their expression level in 27 walnut samples. ACCase, KASII,LACS, FAD3 and LPAAT might be highly regulated by these TFs.Considering that there are few studies on the regulation of these TFs in lipid biosynthesis, more research about molecular biology experiment such as gene knockout will enhance our understanding on their interaction in enzyme activity. Our study reveals the potential molecular mechanism of oil accumulation and lays the foundation for metabolic engineering and molecular breeding of walnut to enhance oil accumulation and FA composition.

Declaration of competing interest

The authors declare that they have no competing interests.

Acknowledgments

This work was supported by Major scientific and technological projects of Xinjiang Production and Construction Corps (2017DB006 and 2020KWZ-012). We would like to thank Prof. Liang Liu from Department of Statistics, The University of Georgia, United States for his help on writing and revising the manuscript, Prof. Qikang Gao of Center of Analysis and Measurement, Faculty of Agriculture, Life and Environment Science, Zhejiang University for his assistance in experiment of quantitative RT-PCR analysis.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://doi.org/10.1016/j.fshw.2022.04.020.