Combined Transcriptome and Lipidomic Analyses of Lipid Biosynthesis in Macadamia ternifolia Nuts

Macadamia nuts are considered a high-quality oil crop worldwide. To date, the lipid diversity and the genetic factors that mediate storage lipid biosynthesis in Macadamia ternifolia are poorly known. Here, we performed a comprehensive transcriptomic and lipidomic data analysis to understand the mechanism of lipid biosynthesis by using young, medium-aged, and mature fruit kernels. Our lipidomic analysis showed that the M. ternifolia kernel was a rich source of unsaturated fatty acids. Moreover, different species of triacylglycerols, diacylglycerol, ceramides, phosphatidylethanolamine, and phosphatidic acid had altered accumulations during the developmental stages. The transcriptome analysis revealed a large percentage of differently expressed genes during the different stages of macadamia growth. Most of the genes with significant differential expression performed functional activity of oxidoreductase and were enriched in the secondary metabolite pathway. The integration of lipidomic and transcriptomic data allowed for the identification of glycerol-3-phosphate acyltransferase, diacylglycerol kinase, phosphatidylinositols, nonspecific phospholipase C, pyruvate kinase 2, 3-ketoacyl-acyl carrier protein reductase, and linoleate 9S-lipoxygenase as putative candidate genes involved in lipid biosynthesis, storage, and oil quality. Our study found comprehensive datasets of lipidomic and transcriptomic changes in the developing kernel of M. ternifolia. In addition, the identification of candidate genes provides essential prerequisites to understand the molecular mechanism of lipid biosynthesis in the kernel of M. ternifolia.


Introduction
Macadamia is a nut tree with high-quality kernel oil. It belongs to the family Proteaceae, known as Australian or sometimes Hawaiian nuts, and has been grown worldwide in tropical and subtropical regions [1,2]. The genus Macadamia has four types of species that are commercially grown to produce nuts. These species are Macadamia integrifolia, M. tetraphylla, M. ternifolia, and M. jansenii [3]. However, M. integrifolia and M. tetraphylla have more significance due to their high-quality oils [4]. Approximately 44,000 metric tons of macadamia kernels are produced per year in the world and 14,100 metric tons per year are produced alone in Australia [5]. Because of its simple production technology, survival in low-fertility soil, high profitability, and low-temperature resistance compared to other traditional crops of tropical regions [6], the planting area of macadamia promptly increased over the past few years worldwide. China is a world-leading country in terms of planting area and is expected to reach more than 15,000 ha in the southern regions. About

RNA Sequencing, De Novo Assembly, and Annotations
To obtain the genome-wide transcriptional changes in the fruit kernel of M. ternifolia, kernels from three different fruit developmental stages (young (S1-1,1-2,1-3), medium-aged (S2-1,2-2,2-3), and mature (S3-1,3-2,3-3) fruit) were selected to perform RNA sequencing. The mean raw reads in S1 were 54,363,432, whereas an average of 48,083,672 and 52,648,740 raw reads were obtained in the S2 and S3 stages, respectively ( Table 1). The filtering of low-quality reads produced an average of 52,621,607, 45,606,890, and 50,793,072 clean reads in the libraries of S1, S2, and S3, respectively. A total of 7.9 Gb of clean data was obtained from the S1 stage, while 6.8 and 7.6 Gb were acquired from stages S2 and S3, respectively. The Q30% was 93% and the GC content percentage was almost 45 in our sequencing. This result suggested that we had produced high-quality data that was suitable for performing transcriptome downstream analysis. Trinity software was used to assemble the 307,336,019 total bases. The total number of unigenes assembled was 279,972. The average length was 1098 nt, with N50 and N90 length values of 1633 and 432 nt, respectively (Supplementary Materials Table S1). The length distribution range of the unigenes was 200 to 2000 nt ( Figure 1A). Due to the lack of a whole reference genome of M. ternifolia, all unigenes were blasted by using Nr, Swiss-Prot, KEGG, KOG, GO, and Trembl public databases to obtain annotations ( Figure 1B). The majority of unigenes showed functional annotation in the Nr (163,325) and Trembl (162,849) databases. There were 119,590, 109,208, 96,830, 133,334, and 110,274 unigenes that were annotated in KEGG, SwissProt, KOG, GO, and Pfam, respectively. The abundance of each unigene was measured with fragments per kilobase per million reads (FPKM) values. The experiment and sampling suitability was determined with principal component analysis (PCA) and correlation matrixes. The PCA showed that the S1, S2, and S3 fruit kernel stages had significant variations with three distinct clusters ( Figure 2A). However, a weak variation among replicates of each stage was observed. The difference was equally covered by PC1 and PC2, which explained 22.9 and 22.6% of the total variation, respectively. Moreover, a strong correlation (R 2 > 0.90) was observed between the replicates of each stage ( Figure 2B). In contrast, a weaker correlation was determined between different fruit kernel stages. The correlation value between S1 and S2 was around 0.50, whereas S1 and S3 showed a much weaker correlation (R 2 = 0.15). The lower degree of correlation between the three stages indicated strong transcriptional changes.

Transcriptome Changes in Different Fruit Kernel Stages
To explore the gene expression changes in developing fruit kernels, we compared the young (S1), medium-aged (S2), and mature (S3) fruit kernel stages. Differentially expressed genes (DEGs) were obtained for two compared samples using p ≤ 0.05 and fold changes (ratio ≥ 1 or ratio ≤ −1) as thresholds ( Figure 3A). The comparison between S1

Transcriptome Changes in Different Fruit Kernel Stages
To explore the gene expression changes in developing fruit kernels, we compared the young (S1), medium-aged (S2), and mature (S3) fruit kernel stages. Differentially expressed genes (DEGs) were obtained for two compared samples using p ≤ 0.05 and fold changes (ratio ≥ 1 or ratio ≤ −1) as thresholds ( Figure 3A). The comparison between S1 and S2 resulted in 12128 downregulated and 12672 upregulated DEGs. A total of 36135 DEGs were determined in the comparison of S1 and S3. From these, 18546 genes were downregulated and 17589 genes were upregulated. When we compared S2 and S3, 14,177 genes were downregulated and 11421 genes were upregulated, with a total of 25598 DEGs. The

Transcriptome Changes in Different Fruit Kernel Stages
To explore the gene expression changes in developing fruit kernels, we compared the young (S1), medium-aged (S2), and mature (S3) fruit kernel stages. Differentially expressed genes (DEGs) were obtained for two compared samples using p ≤ 0.05 and fold changes (ratio ≥ 1 or ratio ≤ −1) as thresholds ( Figure 3A). The comparison between S1 and S2 resulted in 12,128 downregulated and 12,672 upregulated DEGs. A total of 36,135 DEGs were determined in the comparison of S1 and S3. From these, 18,546 genes were downregulated and 17,589 genes were upregulated. When we compared S2 and S3, 14,177 genes were downregulated and 11,421 genes were upregulated, with a total of Life 2021, 11, 1431 5 of 18 25,598 DEGs. The distribution of DEGs revealed that a major portion of the DEGs was overlapped but less stage-specific DEGs were identified among the different fruit kernel stages ( Figure 3B). For example, the number of overlapped DEGs in S1 vs. S2 and S1 vs. S3 was 13,185. In addition, 3316 DEGs were overlapped among the S1, S2, and S3 fruit kernel stages. Globally, these results indicate major gene expression dynamics for fruit kernel development in M. ternifolia.
Life 2021, 11, x FOR PEER REVIEW 5 of 19 and S2 resulted in 12,128 downregulated and 12,672 upregulated DEGs. A total of 36,135 DEGs were determined in the comparison of S1 and S3. From these, 18,546 genes were downregulated and 17589 genes were upregulated. When we compared S2 and S3, 14,177 genes were downregulated and 11,421 genes were upregulated, with a total of 25,598 DEGs. The distribution of DEGs revealed that a major portion of the DEGs was overlapped but less stage-specific DEGs were identified among the different fruit kernel stages ( Figure 3B). For example, the number of overlapped DEGs in S1 vs. S2 and S1 vs. S3 was 13,185. In addition, 3316 DEGs were overlapped among the S1, S2, and S3 fruit kernel stages. Globally, these results indicate major gene expression dynamics for fruit kernel development in M. ternifolia.

Functional Enrichment Analysis of DEGs
To obtain the complete functional information about transcripts that regulate fruit kernel development, GO and KEGG enrichment analyses were performed by using p ≤ 0.05 as the criterion for significant enrichment. In the GO analysis, the majority of DEGs in the S1 and S2 stages were involved in the molecular function of beta-glucosidase and oxidoreductase activity ( Figure 4A). Meanwhile, in pathway enrichment, metabolic and biosynthesis of secondary metabolites had more significance and total genes. A total of 1595 DEGs were enriched in each of these pathways ( Figure 4B). High numbers of oxidoreductase-and peroxidase-activity-associated genes were observed in the DEGs from the comparison of S1 vs. S3 ( Figure 5A). DEGs with annotations regarding metabolism and biosynthesis of secondary metabolites exhibited predominant pathway enrichment ( Figure 5B). Functional enrichment analysis of the S2 vs. S3 stages revealed that the majority of DEGs performed the biological function of cellulose biosynthesis, the molecular function of beta-glucosidase, and the cellular function of the kinesin complex ( Figure 6A). The majority of DEGs in the S2 vs. S3 stages had the same pathway enrichment as the S1 vs. S2 and S1 vs. S3 stages ( Figure 6B). These findings suggested that genes involved in the functional activity of beta-glucosidase and oxidoreductase made significant contributions to the fruit kernel development in M. ternifolia.

Functional Enrichment Analysis of DEGs
To obtain the complete functional information about transcripts that regulate fruit kernel development, GO and KEGG enrichment analyses were performed by using p ≤ 0.05 as the criterion for significant enrichment. In the GO analysis, the majority of DEGs in the S1 and S2 stages were involved in the molecular function of beta-glucosidase and oxidoreductase activity ( Figure 4A). Meanwhile, in pathway enrichment, metabolic and biosynthesis of secondary metabolites had more significance and total genes. A total of 1595 DEGs were enriched in each of these pathways ( Figure 4B). High numbers of oxidoreductase-and peroxidase-activity-associated genes were observed in the DEGs from the comparison of S1 vs. S3 ( Figure 5A). DEGs with annotations regarding metabolism and biosynthesis of secondary metabolites exhibited predominant pathway enrichment ( Figure 5B). Functional enrichment analysis of the S2 vs. S3 stages revealed that the majority of DEGs performed the biological function of cellulose biosynthesis, the molecular function of beta-glucosidase, and the cellular function of the kinesin complex ( Figure 6A). The majority of DEGs in the S2 vs. S3 stages had the same pathway enrichment as the S1 vs. S2 and S1 vs. S3 stages ( Figure 6B). These findings suggested that genes involved in the functional activity of beta-glucosidase and oxidoreductase made significant contributions to the fruit kernel development in M. ternifolia.

Figure 7.
Heatmaps of the lipid content that had the highest differential accumulation among the different stages of the M. ternifolia fruit kernels: (A) in the comparison of the young (S1) and medium-aged (S2) stages, (B) in the comparison of the young (S1) and mature (S3) stages, and (C) in the comparison of the medium-aged (S1) and mature (S3) stages.

Lipid-Related Gene Expression Changes during Fruit Kernel Development
To find the key genes associated with fruit kernel development, further mining of metabolic genes was achieved through the selection of DEGs having a fold change value ≥ 5 or ratio ≤ −5 as thresholds. The Venn diagram revealed 282 stage-specific highly expressed DEGs between S1 and S2 (Supplementary Materials Figure S2). The S1 and S3 stages comparison had 246 stage-specific DEGs. It was observed that S2 vs. S3 had 412 stage-specific DEGs with a high fold change ratio. Moreover, 400 overlapped DEGs among all stages were considered as a regulator of the metabolic process during fruit kernel growth. The expression trends of annotated genes associated with lipid biosynthesis are shown in Table 2. Significant low expression levels of four genes related to Figure 7. Heatmaps of the lipid content that had the highest differential accumulation among the different stages of the M. ternifolia fruit kernels: (A) in the comparison of the young (S1) and medium-aged (S2) stages, (B) in the comparison of the young (S1) and mature (S3) stages, and (C) in the comparison of the medium-aged (S1) and mature (S3) stages.
The comparative analysis between the S1 and S3 stages revealed 261 lipid subclasses with differential concentrations (Supplementary Materials Table S4). From these, 187 had a high accumulation in S3 compared to S1 and only 74 subclasses of lipid showed a lower content in S3. Quantification showed that TAGs, including 49:3, 52:2, 52:3, 54:2, 56:2, 56:3, 56:4, 58:1, 60:2, and 62:2, had their highest levels in the S3 stage as compared with S1 ( Figure 7B). Meanwhile, the top ten lipids with the lowest accumulation content in the S3 stage belonged to different subclasses ( Figure 7B) Figure 7C). It was noticed that there were few stage-specific lipids but there were more overlapped or common lipid profiles among the kernel developmental stages (Supplementary Materials Figure S2). Overall, the results of the lipid profiling revealed that during fruit kernel growth, triacylglycerols and diacylglycerol related lipids accumulation were increased. In contrast, phosphatidylethanolamine-associated lipids showed decreased accumulation.

Lipid-Related Gene Expression Changes during Fruit Kernel Development
To find the key genes associated with fruit kernel development, further mining of metabolic genes was achieved through the selection of DEGs having a fold change value ≥ 5 or ratio ≤ −5 as thresholds. The Venn diagram revealed 282 stage-specific highly expressed DEGs between S1 and S2 (Supplementary Materials Figure S2). The S1 and S3 stages comparison had 246 stage-specific DEGs. It was observed that S2 vs. S3 had 412 stage-specific DEGs with a high fold change ratio. Moreover, 400 overlapped DEGs among all stages were considered as a regulator of the metabolic process during fruit kernel growth. The expression trends of annotated genes associated with lipid biosyn-  Table 2. Significant low expression levels of four genes related to 3-ketoacyl-acyl carrier protein reductase (unigene 73227, unigene 211387, unigene 73228, unigene 26385) were observed during the S3 fruit kernel developmental stage compared to the S1 stage. The expression levels of five beta-glucosidase genes, namely, unigene 170150, unigene 81638, unigene 81642, unigene 81652, and unigene 91684, were highest during the S2 stage. Similar results were obtained for unigene 150581, unigene 150580, and unigene 140441, which encode the enoyl-acyl-carrier-protein reductase (NADH) enzyme. In addition, the diacylglycerol-kinase-2-linked unigene 147364 had an unpredictable expression during the fruit kernel development stage. A significantly altered expression was also determined for pyruvate-kinase-2-annotated unigene 63669, unigene 63675, and unigene 86831 during fruit kernel growth.  Triacylglycerols are the key resources of energy in the seeds. Two unigenes encoding glycerol-3-phosphate acyltransferase and denoted as unigene 152617 and unigene 128368 were found with low expressions during the late stage of fruit kernel growth in our study. Moreover, a significant change in the expression trends of some genes related to GDSL esterase lipase was obtained during the fruit kernel developmental stages. For example, unigene 240247, unigene 82774, unigene 183098, unigene 38911, unigene 113245, and unigene 238268 showed increased expression levels during the S1 stage compared with S2 and S3. In contrast, unigene 213504, unigene 62756, unigene 6967, and unigene 103298 had higher levels of expression in the S2 stage than S1 and S3. Globally, the expression of linoleate 13S-lipoxygenase 2-1 gene family members, such as unigene 9629, unigene 66716, unigene 54520, unigene 59928, unigene 213589, unigene 37584, unigene 34267, unigene 13471, unigene 5404, and unigene 139752, were significantly downregulated during the late stages of fruit kernel growth. Various genes involved in phospholipids regulation were determined with altered expression levels during the fruit kernel growth. Among these, unigene 190439, unigene 227869, and unigene 171746 had higher expression values in the S2 stage compared with S1 and S3, whereas unigene 158347 and unigene 122230 showed upregulated expression during the S1 (early fruit kernel growth) stage. Overall, this suggests that the stage-specific regulation of genes related to lipid metabolism and catabolism caused the diversity of lipid compositions and concentrations in M. ternifolia nuts.

Validation of Key Genes Involved in Lipid Regulation Using qRT-PCR Analysis
Based on the transcript downstream data analysis, ten genes putatively involved in the lipid regulation were selected to perform a real-time qRT-PCR analysis in order to confirm the accuracy of the RNA-seq. The selected genes showed differential expressions at the different developmental stages of the fruit kernel in M. ternifolia (Supplementary Materials Figure S3). It was observed that the expression trends of all genes were consistent with RNA-seq, thus confirming the reliability of RNA-seq results in our study.

Discussion
Out of four cultivated species, M. integrifolia and M. tetraphylla are widely grown in commercial orchards [40]. There are many constraints that limit the productivity of M. ternifolia. The major challenge is a high level of cyanogenic glycosides present in the kernel, which are toxic and unsuitable for human consumption without steeping and cooking [4,41]. A detailed understanding of the genetic mechanisms underlying lipid composition in the kernel of M. ternifolia will be useful to promote its commercial production level. Through comprehensive lipidomics and transcriptomic data analyses, this study improved our understanding of the molecular factors that influence oil accumulation in M. ternifolia kernels.

Divergence of Lipid Composition in the Kernels of M. ternifolia
A lipid is a form of stored carbon in plant seeds, yields more energy than carbohydrates, and is an integral part of plant organs [29]. The lipid composition was significantly changed in the developing fruit kernels of M. ternifolia. In our study, a total of 407 different types of lipid species were determined at different fruit kernel growth stages. The majority of these lipids belonged to the category of glycerolipids including triacylglycerols (TAGs). TAGs are the primary components of Macadamia oil, accounting for more than 95% of Life 2021, 11, 1431 12 of 18 total lipids. These influence the nutritional characteristics of nut oil through physical and biochemical changes [42]. TAGs are entirely acylated derivatives of glycerol and the most abundant neutral lipid bodies present in seeds [43], senescence leaves, and pollen grains [44,45]. In the mature fruit kernels of M. ternifolia, the most accumulated species of TAGs were long-chain unsaturated fatty acids, including TAG 52:2;1O, TAG 54:2;1O, TAG 56:3;1O, and TAG 56:2;1O. An increase in unsaturated TAGs was previously observed in the later stages of seed development in Paeonia ostii [46], Glycine max [47], and Brassica napus [48]. The stored level of TAGs has a functional linkage to seed germination [49] and stress management [50]. Moreover, the nuts of Macadamia species are usually consumed in roasted form, which can modify the lipid profiles due to oxidation [51]. Excessive oxidation in unsaturated TAGs produces free fatty acids due to esterases and lipases through the activity of hydrolytic enzymes [52]. Free fatty acids are more vulnerable to oxidation during nut roasting, which causes adverse chemical changes through rancidity and the deterioration process, which ultimately leads to rejection by consumers due to a lower quality of nuts [53]. In a recent study, it was observed that roasting had a non-significant impact on the TAGs composition and it ultimately improved the flavor, aroma, color, texture, and appearance of macadamia nuts [42]. Among other detected glycerolipids in the developing fruit kernels of M. ternifolia, different species of diacylglycerols (DAGs) have a high abundance. In particular, DAGs 40:4 and 42:3 species had a high level of accumulation. DAGs are precursors of glycolipids, are swiftly phosphorylated to produce phosphatidic acid (PA) by diacylglycerol kinase enzyme, and act as a potential secondary messenger in plants [54]. The high contents of TAGs and DAGs in the fruit kernel of M. ternifolia may be useful to improve the oil quality through the biosynthesis of other essential fatty acids and lipids.
Phosphatidylethanolamine (PE) is a polar lipid belonging to the phospholipid category, is dominant in the plasma membrane, and plays an important role in cold stress tolerance [55]. In our study, the quantitative level of many PE species was dramatically reduced at the mature fruit kernel stage. Specifically, PE 35:3, PE 36:5, PE 40:3, PE 42:4, and PA 36:3 had higher levels of accumulation in the young fruit kernel of M. ternifolia. Previous research showed that Pes, together with phosphatidylcholine species, alleviates chilling injury to protect cell membrane damage in cold conditions [56]. PA is a type of glycerophospholipid that mediates the biosynthesis of phospholipids and glycolipids in plants. Moreover, it is involved in secondary signal transduction during abiotic stress, including salinity [57], drought [58], and cold [59]. Our analysis further reported a differential level of ceramides (CERs) in the medium fruit kernel stage, in particular, high levels of Cer 44:1 and Cer 42:2 contents were detected. CERs are sphingolipids that are essential components of the membrane in plants. Functional studies showed that CERs are involved in responses to biotic and abiotic stimuli. However, little is known about their signaling mechanism during the programmed cell death of certain bacteria and fungi pathogens [60]. PEs, PAs, and CERs have a critical role in the defense against various stressful conditions. A similar function in developing the kernel of M. ternifolia might be linked with a higher level of lipid stores in the nut. The data of lipidomic analysis identified the presence of some species of acylated sterol glycosides, phosphatidylmethanol, diacylglyceryl glucuronide, and digalactosyldiacylglycerol in the kernel of M. ternifolia. However, their non-differential trend suggested a minor role for the oil quality index. In contrast, a higher abundance of oil species, such as TAGS, DAGs, PEs, and PCs, in the mature nut is of substantial importance to the improve oil quality in M. ternifolia nut.

Lipid Biosynthesis in Fruit Kernels of M. ternifolia
Lipid biosynthesis in nuts is predominantly altered by the activity of genes that are involved in the de novo biosynthesis of fatty acids, synthesis of triacylglycerol, and formation of oil bodies [35]. The genetic mechanism of fatty acid and lipid synthesis leading to TAG is well known in Arabidopsis [61]. The substrate for all synthesized fatty acid pools is acetyl-CoA. The enzyme pyruvate kinase (PK) in the plastids generates pyruvate, which is catalyzed by the plastidial pyruvate dehydrogenase complex to form acetyl-CoA for the de novo production of fatty acids [62]. The altered expression of PK2 genes in young fruit kernels of M. ternifolia suggests their vital role in the generation of diverse fatty acids. Functional genomics experiments on PK2 genes in rice stated their importance in fatty acid and starch metabolism [63]. Among the fatty acid synthase system, genes associated with 3-ketoacyl-acyl carrier protein reductase (KAR) or synthase and enoyl-acyl-carrier-protein reductase (NADH) (ENR) regulate the synthesis of long-chain unsaturated fatty acids in plants seeds [64]. In the formation of long-chain fatty acids, the KAR enzyme reduces 1,3-ketobutyryl-ACP to 3-hydroxybutyryl-ACP, which extends a hydrocarbon chain of the acyl group up to a 16-or 18-carbon chain through a series of cycle reactions [65]. Similarly, ENR is predicted to be involved in the last step of fatty acid elongation in plants, bacteria, and mammals [66]. The altered expression of genes involved in the functional activity of KAR and ENR enzymes probably contributes to the development of unsaturated fatty acids in M. ternifolia. In this way, it might be interlinked with the unique biochemical and physical properties of nut oil. However, further research evidence is needed to support this statement. In our lipidomic data analysis, TAGs, DAGs, PEs, and PAs are the most synthesized forms of lipids present in the fruit kernel of M. ternifolia. The de novo lipid biosynthesis of M. ternifolia most probably starts with the catalytic activity of glycerol-3phosphate acyltransferase (GPAT) enzyme. The GPAT initiates the preliminary reaction of TAG and phospholipids biosynthesis in the endoplasmic reticulum. It further converts glyerol-3-phosphate into lysophosphatidic acid, which produces phosphatidic acid. The dephosphorylation of phosphatidic acid forms DAG and finally TAG [67]. Furthermore, the oleoyl CoA transferred from plastids is used as a substrate to produce membrane lipids (PCs, PEs, and PIs) or the storage lipid TAG in M. ternifolia. However, further metabolic studies will be essential to prove these findings.
Several studies have reported the role of phosphoglycerolipids in lipid signaling [68,69]. The key genes involved in lipid signal transduction include diacylglycerol kinase (DGK), phosphatidylinositol (PI), nonspecific phospholipase C (NPC), phospholipase C (PLC), phospholipase D (PLD), and PA [70]. Interestingly, we observed substantially altered transcription levels of GPAT, DGK, PI, and NPC in developing fruit kernel M. ternifolia. This indicates their key roles not only in higher TAGs, DAGs, PE, and PA accumulation but also in lipid signaling during the growth of M. ternifolia kernel. It was reported that enzyme GDSL-type esterase/lipase can hydrolyze thioesters, aryl esters, phospholipids, and amino acids. In this way, they altered the mechanism of secondary metabolism [71]. A significant change in expression trends of genes related to the GDSL esterase lipase enzyme probably performed a similar function in developing the M. ternifolia fruit kernel. Furthermore, our research also identified several genes associated with the linoleate 9S-lipoxygenase (LOX) enzyme, particularly LOX6, which had a significantly modified expression in the M. ternifolia fruit kernel. This enzyme utilizes linoleic acid/linolenic acid as a substrate to produce unsaturated fatty acids, which improves the oil content and quality for human health [72]. In addition, LOX genes play a critical role in the enzymatic oxidation of polyunsaturated fatty acids to prevent germination during seed aging [73]. The β-glucosidases (BGlu) enzymes have numerous roles in plants, such as defense, cell wall lignification, and phytohormone activation. These enzymes are also responsible for the release of scented volatile compounds in plants [74,75]. The increased expression of BGlu genes in the developing fruit kernel of M. ternifolia suggests an association with the diversity of functions that directly or indirectly affect oil quality.
Concisely, our comparative lipidomic and transcriptomic data analysis identified key regulatory genes associated with oil content and quality in M. ternifolia kernels. Due to the wide public acceptance of Macadamia nuts as a portion of quality food, these results are vital for promoting the nutritional and commercial value of M. ternifolia kernels. Our results provide a baseline for further transgenic and genetic engineering research with the aim to identify oil candidate genes and their regulators. Hence, further works are needed to boost the application of M. ternifolia nuts in the food, pharmaceutical, and cosmetic industries.

Plant Materials, RNA Extraction, and Sequencing
The plant material used in this study was collected from M. ternifolia plants cultivated in the Yunnan province of China. The developing fruit kernels comprised of young (S1), medium-aged (S2), and mature (S3) fruiting stages were picked in three biological repeats. The harvested samples were immediately placed in liquid nitrogen and stored at −80 • C before further processing. Total RNA was isolated from all nine samples with a Trizol reagent (Invitrogen, San Diego, CA, USA) following the manufacturer's standard protocol. RNase-free DNase I (TaKaRa, Kyoto, Japan) was used to purify the total RNA. The RNA quality and quantity were determined by 1% agarose gel and Nanodrop spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA), respectively. Sequencing libraries were prepared with the Illumina TruSeq Standard mRNA library preparation kit following the manufacturer's standard protocol. After the quality tests, the Illumina HiSeq Ten X platform (Illumina Inc., San Diego, CA, USA) was used according to the recommended protocol to perform de novo pair-end RNA sequencing.

Transcript Assembly, Annotation, and Analysis of DEGs
By removing the adaptors, ambiguous bases, and low-quality reads, the original sequencing data in the form of raw reads were filtered to obtain high-quality clean reads. Because of the non-availability of the reference genome of M. ternifolia, high-quality clean reads were assembled into unigenes for subsequent analysis with Trinity 2.6.6 software [76]. The BLAST software was used to compare the unigenes sequence with the Kyoto Encyclopedia of Genes and Genomes database (KEGG), NCBI non-redundant protein sequences (NR), Gene Ontology (GO), euKaryotic Ortholog Groups (KOG), Protein family (Pfam), Trembl, and Swiss-Prot databases to retrieve the unigenes' functional annotations [77]. The expression levels of the mapped transcript were determined in the form of fragments per kilobase of transcript per million fragments mapped (FPKM) values. The R package DESeq2 version 1.22.2 [78] was used to obtain the total number of differential genes and the number of upregulated and downregulated genes among two different samples. The significant criteria for DEGs included log2 fold change ≥1 and ≤−1 and a false discovery rate < 0.05. The KEGG for linking genomes to life and the environment platform was used to perform the pathway enrichment analysis [79].

Real-Time qRT-PCR Analysis
To perform the qRT-PCR analysis for ten selected genes, the cDNAs were produced with iScript cDNA Synthesis Kit (BIO-RAD). Primer Premier 5.0 was used to design the gene-specific reverse and forward primers (Table S6). The reaction mixture in three technical repeats was prepared with iTaq Universal SYBR Green Supermix 50 mL (BIO-RAD) kit. The quantitative PCR analysis was conducted in triplicate biological replicates on a BIO-RAD iCycler iQ5 real-time PCR system. The Actin2 gene was used as an internal control. The running and data analysis protocol of the qRT-PCR was followed as detailed in a previous report [80].

Lipid Extraction and Analysis Using UPLC-MS/MS
Different fruit kernels consisting of young (S1), medium-aged (S2), and mature (S3) fruiting stages were utilized to extract lipids. In brief, a 15 mg sample was taken to make a fine powder. First, 2 mL of precooled (−20 • C) methanol was added, followed by 4 mL dichloromethane, and the mixture was vortexed for 1 h. After 10 min of incubation at 4 • C and sonication for 10 min in an ultra-sonication bath, 1.6 mL double-distilled water was added. Samples were then centrifuged for 5 min. Two phases were produced after the centrifugation. The upper clean and residual-free lipophilic phase was carefully removed. The lipophilic phase was then dried under a vacuum. Later, 1 mL isopropyl alcohol was added in the lipophilic phase and used for further lipidomic analysis. The lipid profiling was executed on a 6600 plus AccurateMass Q-TOF (AB Sciex TripleTOF ® 6600, Framingham, MA, USA) mass spectrometer system. The process of liquid chromatography was completed using UPLC with an autosampler LC-30A (Shimadzu Corporation) column (100 × 2.1 mm, 2.6 m). The lipid measurement protocol was followed as in the description of Wan et al. [81]. The MSDIAL ver. 4.0, PeakView2.1, and MultiQuant 3.0 were used to perform the quantitative and qualitative data analyses of the lipids [82].

Conclusions
Our results revealed lipid and transcript profiles at various developing fruit kernel stages of M. ternifolia nuts. Moreover, the key lipids and their putative candidate genes were identified. In particular, GPAT, DGK, PI, NPC, PK, KAR, and LOX enzymes encoding genes were most probably rate-limiting for the biosynthesis of glycerolipids and phospholipid in M. ternifolia nuts. These results help to understand lipid diversity and provide first insights into the genetic mechanism of lipid accumulation. Gene functional characterization and further transcript data from different genotypes are mandatory to explore the regulatory network of lipid biosynthesis in M. ternifolia nuts.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/life11121431/s1, Table S1: Statistics of assembly for the different fruit kernel stages of M. ternifolia. Table S2: Total diversity of lipids and their content level in the different fruit kernel stages of M. ternifolia. Table S3: Total contents of lipids differentially accumulation in the comparison of the S1 vs. S2 stages. Table S4: Total contents of lipids differentially accumulation in the comparison of the S1 vs. S3 stages. Table S5: Total contents of lipids differentially accumulation in the comparison of the S2 vs. S3 stages. Table S6: Primer sequences used in the qRT-PCR experiment. Figure S1: Distribution of the overlapped and unique number of lipid species among the different fruit kernel stages of M. ternifolia. Figure S2: Distribution of the overlapped and unique highly expressed metabolic DEGs among the different fruit kernel stages of M. ternifolia. Figure S3: qRT-PCR validation of selected differentially expressed genes. The Actin2 gene was used as the internal control and the transcript levels in S2 and S3 were expressed as compared with S1. Young (S1), medium-aged (S2), and mature (S3) fruiting stages.