Phenotypic variation and flavonoid metabolic features
A white seed coat peanut mutant (wsc) was identified from a mutation population derived from an elite pink seed coat variety Zhonghua 16 (cultivated by Oil Crops Reasearch Institute of Chinese Academy of Agricultural Sciences) once treated with 60Co. The WT accumulated pink anthocyanins and no obvious color change was observed in wsc during development (Fig. 1A). A histochemical analysis was therefore performed to determine the nature of polyphenol compounds in developing testae between wsc and WT. Toluidine blue O (TBO) staining of developing testae transverse sections revealed the distribution of polymeric phenolic compounds and showed no significant differences between wsc and WT in the three seed development stages assessed (Fig. 1B). A series of flavonoid metabolome profiles were therefore conducted between wsc and WT. We detected a total of 199 flavonoids in testae samples including 19 anthocyanins, 23 flavanones, 63 flavones, 32 flavone C-glycosides, 38 flavonols, three flavonolignans, one hydroxycinnamoyl derivative, 15 isoflavones, and five proanthocyanidins (Additional file 1). Data revealed that the total flavonoid content in wsc was 1.67 times higher than in WT at DAF20, declining to just 1.43 times and 0.97 times that much at DAF40 and DAF60, respectively (Fig. 1C). It is clear that both wsc and WT accumulated flavonoids while developing; flavonoid content peaked at DAF40 in wsc compared with WT samples at DAF60 (Fig. 1C).Anthocyanins, proanthocyanidins, and flavonols, all the dominate flavonoids, were enhanced throughout development except flavonolignan. The anthocyanins detected in this study included cyanidin, delphinidin, malvidin, pelargonidin, peonidin, rosinidin, and derivatives, data showed that peonidin, delphinidin,cyanidin O-diacetyl-hexoside-O-glyceric acid, cyanidin 3,5-O-diglucoside (cyanin), pelargonin, cyanidin 3-O-rutinoside (keracyanin),andcyanidin levels all increased throughout the development of testae in WT while remained at very low stable levels in wsc. Levels of the four glycosylated anthocyanins (i.e., malvidin 3,5-diglucoside (Malvin), cyanidin 3-O-glucoside (Kuromanin), cyanidin 3-O-malonylhexoside, and rosinidin O-hexoside) all increased in wsc during the two early developmental stages and then decreased to equal amounts in the WT at DAF60 (Fig. 1D), only malvin keeps the increasing trend in wsc. At the same time, procyanidin A1, procyanidin A2, procyanidin A3, procyanidin B1, and procyanidin B2 all occurred at higher levels in WT samples, especially at DAF40 and DAF60, while procyanidin contents in wsc remained three orders of magnitude lower (Fig. 1E). The content of flavonols was also redirected over the course of this experiment between wsc and WT, and eight of these were up-regulated while nine were down-regulated in wsc (Fig. 1F). Indeed, both the majority and total amounts of isoflavones were enhanced in wsc (Fig. 1G), while flavonolignan volumes did not change significantly between WT and wsc (Fig. 1H). Flavones were also redirected between wsc and WT; the total contents of these compounds were surprisingly enhanced in wsc (Figs. 1C, 1I). A total of 13 flavone C glycosides were up-regulated in this analysis while eight were down-regulated in wsc although total contents obviously increased overall (Fig. 1C and 1J). Considerable increases in four flavanones were also seen during the two early developmental stages, three were significantly strengthen at DAF60, while two of the other four differentiated flavanones had declined at DAF60 as well as one at DAF20 and DAF 40, respectively (Fig. 1K).
An untargeted lipid metabolome analysis was also performed as part of this study identifying 6,424 differential known metabolites (Additional file 2), KEGG enrichment analysis of differentiated metabolites (Additional file 2) validated the enrichment of flavonols and also showed that isoflavonoid biosynthesis between wsc and WT remained quite consistent across the flavonoid metabolome.Differential over-representation was also revealed in several metabolic pathways, including brassinosteroid, cutin, suberin, and wax biosynthesis, as well as arginine, proline, starch, sucrose, and galactose metabolism (Additional file 2).
A seed metabolome analysis was performed to investigate whether, or not, wsc influence embryo nutrient quality. A total of 439 metabolics were identified at differential levels within embryos of wsc and WT (Additional file 3) with the top five enriched pathways recorded as metabolic pathways, biosynthesis of secondary metabolites, isoflavonoid biosynthesis, flavonoid biosynthesis, and flavone and flavonol biosynthesis (Additional file 4A). Particularly, all differential metabolites from isoflavonoid biosynthesis were up-regulated during this analysis (Additional file 4A). At the same time, five of the six metabolites from flavone and flavonol biosynthesis were enhanced, and the contents of 3,7-Di-O-methylquercetin, quercetin 3-O-[beta-D-xylosyl-(1->2)-beta-D-glucoside], and kaempferide in wsc rose to more than five times those seen in WT (Additional file 3). All metabolites involved in phenylpropanoid biosynthesis and phenylalanine metabolism declined at different levels (Additional file 3), while both the two involved in cutin, suberin, and wax biosynthesis were down-regulated in wsc (Additional file 4A). The results of this analysis are consistent with those from testae flavonoid metabolomes. Testae play crucial role in transporting, as the main energy storage form, fatty acid contents in wsc and WT were tested. Results reveal that almost all fatty acid (FA) components together with the total FA contents decreased over the course of this study between 5.47%–19.41% (Additional file 4B).
Global transcriptome analysis reveals the involvement of multiple biological processes during testa development
In order to identify the genes that underlie the regulation of testa pigment biosynthesis, metabolic pathways, and related hormone signaling, RNA-seq were performed on both wsc and WT samples collected at DAF20, DAF40, and DAF60. A stringent FDR value ≤ 0.001 and a fold change ≥ 2 were used as thresholds to identify DEGs. Results revealed a total of 17,428 DEGs between wsc and WT (Fig. 2A; Additional file 5). A total of 6,397DEGs, 10,471 DEGs, and 9,059 DEGs were identified at DAF20, DAF40, and DAF60, respectively (Fig. 2A). Clearly more up-regulated genes (3,725) were present at DAF20 than down-regulated ones (2,672), while the number of up-regulated genes (7,036) at DAF40 was much more than those down-regulated (3,395). Results revealed 7,238 up-regulated genes and 1,821 down-regulated ones at DAF60 (Fig. 2A).
In order to investigate altered functional categories between wsc and WT lines, a series of GO classifications were applied for further DEG analysis. The results of these comparisons (Fig. 2B) revealed that the GO-term distributions of DEGs correspond with biological processes, molecular functions, and cellular components (Additional file 6). Data revealed that metabolic process, cellular process, single-organism process, biological regulation, and localization were all dominant biological process categories (Fig. 2B). In terms of GO terms relating to cellular component, most DEGs were correlated with five major biological processes, including cell as well as cell part, membrane, membrane part, and organelle (Fig. 2B). Catalytic activity, binding, transporter activity, structural molecule activity, and nucleic acid binding transcription factor activity were the top categories annotated for the molecular function (Fig. 2B).
The subsequent use of KEGG pathway analysis assigned DEGs to 274, 249, and 269 metabolic pathways, respectively, associated with the three different developmental stages of the two lines. A total list of the metabolic pathways identified in this study is presented in Additional file 7, alongside the top 16 common metabolic and biological pathways identified in the wsc compared with WT (Fig. 2C). These results revealed that metabolic pathways were the most enriched, followed by secondary metabolite biosynthesis, plant hormone signal transduction, carbon metabolism, endocytosis, starch and sucrose metabolism, and amino acid biosynthesis (Additional file 7).
qPCR verification of wsc and WT DEGs
In order to evaluate the accuracy of the RNA-seq data generated in this study, biologically independent quantitative reverse-transcription PCR (qRT-PCR) was utilized to verify DEGs. A total of 17 genes with Fragments per Kilobase Million (FPKM) values ≥ 2 were selected for confirmation (the gene-specific primers used for this analysis are listed in Additional file 8 (Additional file 9A). The results of a linear regression analysis revealed an overall correlation coefficient (R) of 0.81, indicative of a strong relationship between the transcription profiles revealed by RNA-Seq and abundances assayed using qRT-PCR (Additional file 9B).
WSC mutation influences flavonoid metabolic pathways
Flavonoid components determine the color of the peanut testa [21–23]. Flavonoid biosynthesis pathway genes were searched based on their KO identifications within the KEGG database as well as synonyms identified in combined with functional annotations, this enabled the identification of 433 DEGs in flavonoid metabolic pathways (Table 1). The genes involved in the three secondary metabolic pathways (i.e., flavonoid, anthocyanin, and flavones/flavonol biosynthesis) related to pigmentation were analyzed in this study using testa transcripts, the core genes in these pathways were studied in detail, and results demonstrated that most exhibited significant changes in expression levels over the course of this analysis. Indeed, regardless of whether these were EBGs (e.g., CHI) or LBGs (e.g., ANS, UFGT),, all exhibited higher transcript abundances within WT compared to wsc with the exception of 4-coumarate:CoA ligase genes and their flavonol synthase counterparts (Fig. 3). The dihydroflavonols represent a branch point in flavonoid biosynthesis, being the intermediates in the production of both the colourless flavonols, through FLS, and the colored anthocyanins, through DFR. Anthocyanins were down-regulated in wsc, which corresponds to lower level expression levels of DFR genes, while myricetin and kaempferols contents were enhanced in wsc, consistent with FLS expression (Fig. 3). Combining the transcriptome and metabolome information, it could be inferred that DFR might be the target gene for the loss of pink color in white wsc. Competition between FLS and DFR for common dihydroflavonols substrates might partially block the synthesis of anthocyanins and cause the increased production of flavonols such as myricetin and kaempferol, shifting the flavonol:anthocyanin ratio in wsc.
Suberin in wsc testa epidermis
Suberin and its associated waxes play important roles at these tissue and plant–environment interfaces by serving as a barrier controlling the movements of water and solutes [24]. As flavonoids in the testa of wsc changed greatly over the course of this analysis, KEGG pathway involved in cutin, suberin, and wax biosynthesis were further analyzed in this study. Data indicated that most cutin and suberin synthesis pathway genes showed significantly higher levels of gene transcripts in wsc (Fig. 4A and 4B) while the genes involved in wax synthesis exhibited no obvious regularity (Fig. 4C). Five of the six FAR gens involved in the synthesis of aromatic suberin monomers for the aliphatic domain [25] were significantly depressed at either DAF20 or DAF60, eight of 11 HHT1 DEGs were up-regulated, and all CYP86B1 genes exhibited a higher level of expression in wsc (Fig. 4A). Four out of five CYP86A4S genes were classified as DEGs while all four HTH genes exhibited enhanced expression levels compare with WTs (Fig. 4A). Suberin staining revealed the presence of almost no red staining in all three developmental stages of wsc compared with the deeper red staining in WT during development (Fig. 4D).
Primary metabolic differences
Transcriptome analysis revealed that the glycolysis and gluconeogenesis pathway was differentially regulated between the wsc and WT. The expression of coding genes for the 13 enzymes catalyze the 14-step reactions in glycolysis was also enhanced to varying degrees (Additional file 10, 11). All nine DEGs in hexokinase which function as one of the three key enzymes (i.e., hexokinase, phosphofructokinase, and pyruvate kinase) within the glycolysis pathway markedly increased in expression, especially in the two later stages (Additional file 10A, 11A), while the six DEGs in phosphofructokinase exhibited higher expression levels in wsc via different patterns (Additional file 10A, 11B). Two of the 12 pyruvate kinase genes were down-regulated in the wsc, especially during the early stage at DAF20, including the highest expressed gene, Araip.VA90H. The remaining ten pyruvate kinase genes also exhibited enhanced expression in the mutant (Additional file 10A, 11C).
The products of the glycolysis pathway were then catalyzed by the pyruvate dehydrogenase complex and entered into the TCA cycle. As data revealed that 19 DEGs, 26 DEGs, and 37 DEGs were enriched in this cycle at DAF20, DAF40, and DAF60, respectively, the ten enzyme encoding genes that control this process were analyzed in detail and six were found to be differentially expressed in the wsc. Four pyruvate dehydrogenase genes were up-regulated in all three developmental stages in the wsc, while the other three were only enhanced at DAF40 (Additional file 10B, 11D). At the same time, two citrate synthase genes increased and one deceased in the wsc (Additional file 10B, 11E), while all aconitase genes were characterized by activated expression in this mutant, especially the three high expressed genes (Aradu.S0KU9, Aradu.83N8C, and Araip.SB6JF) (Additional file 10B, 11F). Isocitrate dehydrogenase expression was also universally higher in the wsc, reaching its highest level at DAF40 (Additional file 10B, 11G), while both DEGs in succinatedehydrogenase differed at DAF40 with one activated and the other one depressed (Additional file 10B, 11H). Five of the six malate dehydrogenase genes were up-regulated at either DAF40 or DAF 60 (or both), while Araip.M2TZZ was down-regulated in all the three wsc developmental stages (Additional file 10B, 11I).
DEG investigations within amino acid metabolism pathways revealed that glutamate amino acid synthesis occurred at significantly higher expression levels within wsc (Fig. 5A). Five of P5CS (one of the two key enzyme coding gene in proline synthesis) genes were all characterized by stronger levels of expression at all three developmental stages while the other three had higher expression levels in the wsc at the two early stages. Two P5CR (another key enzyme coding gene in proline synthesis) gene were detected in differential expression; Aradu.A85Y5, was characterized by a significantly higher expression level at all three stages while its counterpart, Aradu.B47VE, was only activated at DAF20 and DAF40 (Fig. 5A). Metabolome analysis of testae revealed that the proline content in wsc was three times that of WT while the N-Carbamoylsarcosine content reached 1.8 times that level when compared with WT (Fig. 5B).
Sugars supply plants with energy and functioned in flavonoids accumulation [26, 27]. Three gene families encoding sucrose transporters have been implicated in plant sugar accumulation, including the H+/sucrose antiporters located on tonoplast (TSTs) [28], H+/sucrose symporters (SUTs) [29–31], and clade III of SWEET [32–34]. We compared the expression patterns of these gene families between wsc and WT (Additional file 12) and found seven SUTs to occur at significantly higher levels within wsc than WT even though expression patterns were different; Aradu.M7BX0, Aradu.QA3D9, and Aradu.SQ8X7 were all owned highest levels at the early developmental stage and adjusted to a moderate level; in contrast, Araip.53MMF exhibited a higher expression level within wsc at DAF60 while the other three SUTs were characterized by consistently higher expression levels in mutants at the all three stages (Additional file 12). Nine of the 13 DEGs in SWEETs were enhanced and four decreased in expression (Additional file 12). In the case of the monosaccharide transporter, expression of the two STPs in wsc was dramatically higher than in WT (Additional file 12). The metabolome profile data elucidated that wsc testae included more than ten-fold increases in maltose and levanbiose contents compare with WT (Additional file 12C).
The vital roles of phytohormones and TFs in the white testa phenotype of wsc
Hormone and transcription factors have both been reported to perform functions in flavonoid accumulation [15, 35]. Genes involved in hormone synthesis and signal transductions were compared in this study. Results revealed the presence of five DEGs involved in the brassinosteroid (BR) synthesis pathway, of which three BR6OX1 were up-regulated either at DAF20 or at DAF40 and DAF60; in contrast, the two BAS1 genes which reduce the level of active BRs were significantly down-regulated in wsc at DAF60 (Additional file 13A). Analysis of the testae metabolome uncovered 16 metabolites within the BR biosynthesis pathway, and all of them were up-regulated (Additional file 2). At the same time, 16 DEGs involved the GA synthesis pathway including KAO (3), GA2ox (3), GA3ox (1), and GA20ox (9) were identified in this analysis; the majority of these genes were up-regulated in wsc with the exception of Aradu.BUI3V,Araip.KVM2C, and Aradu.U5UGY that all exhibit opposite expression patterns (Additional file 13B). Further, testae metabolome analysis also supported the enhancement of GA synthesis by identifying three pathway products that increased in content within wsc (Additional file 2). While DEGs involved in JA biosynthesis included two PLA1 (i.e., BGI_novel_G004010 and Aradu.Y096Z),, one AOC (Aradu.H2RVW),, four OPR (i.e., Aradu.FV6YV,Araip.F1ZZD,Aradu.0C0YF,andAraip.22RGE),,and three MFP (i.e., Araip.UPC07,Araip.M6X0I, and Araip.60KBF) were all enhanced in wsc (Additional file 13C). These expression and metabolome data taken together therefore supported the enhancement of BR, GA, and JA synthesis pathways.
A series of DEGs involved in seven hormone signal transduction pathways were also identified in this study, including nine in ABA, 12 in AUX, five in BR, ten in cytokinin (CTK), seven in GA, ten in JA, and four in salicylic acid (SA) (Additional file 14). Seven ABA signaling DEGs were enhanced, including SRK2 gene Aradu.SKS95 that exhibited the highest FPKM value (Additional file 14A). The vast majority of AUX signaling pathway genes increased in expression in wsc; only two deviated from this pattern, Araip.TU273 and Araip.60Q0E, both of which exhibited obvious decrease in mutants at DAF20 (Additional file 14B). Further, four out of five DEGs involved in BR signaling were highlighted in wsc (Additional file 14C); and nine of the ten DEGs involved in JA signal pathway were up-regulated and all of them were MYC2 factors (Additional file 14D). Expression patterns and FPKM values of DEGs in hormone signal transduction pathways provided collective evidence of ABA, AUX, JA, and BR signaling enhancements (Additional file 14).
Flavonoid production is transcriptionally regulated by MYB factors and the MBW complex. A total of 24 R2R3-MYB (nine up and 15 down), 26 bHLH (17 up and nine down), and 21 WD40 (19 up and two down) genes were identified in this analysis (Additional file 15) including homologs (Aradu.CA8XJ and Araip.MHR6K) of AtTT8 which enables strong, cell-specific accumulation of flavonoids in Arabidopsis thaliana [36] and homologs of AtMYB5 (i.e., Aradu.JK51Z and Aradu.WZF00) control outer seed coat differentiation alongside TTG1 and TT2 [37].
White testa phenotype candidate genes revealed by multi-omics analysis
In order to identify candidate genes controlling the white testa mutant phenotype in peanut, we analyzed DEGs in common between WT and wsc at the three different developmental stages. This analysis resulted in the identification of 1,646 unigenes (Fig. 2 and Additional file 16). Observations across whole growth stages and tissues revealed no obvious differences between wsc and WT with the exception of the seed phenotype (Additional file 17), therefore WSC was identified as a seed-specific gene. Previously published data [38] was interrogated to identify seed specific expressed (SPE) common DEGs between wsc and WT which resulted in 86 candidates genes for the testa phenotype (Additional file 18). The wsc was identified in the M3 generation of the mutant population while the segregation ratio of the mutant line between WT and wsc was nearly 3:1 (156:49). Combining SPE common DEGs with the results of the genetic analysis revealed the putative candidate genes Araip.M7RY3 (CSN1), Aradu.R8PMF (MYB) and Araip.MHR6K (bHLH) (Additional file 18). The FPKM value of the Araip.M7RY3 gene in the WT seed coat decreased from 2.15 (DAF20) to 1.15 (DAF40) before falling further to 0.19, while the value for this gene in wsc increased from 7.55 (DAF20) to 12.21 (DAF40) and then rose further to 12.84 (DAF60). It is annotated that Araip.M7RY3 is homolog of a COP9 signalosome complex subunit 1(CSN1) encoding gene. Previous studies have revealed that several subunits of the COP9 signalosome complex are involved in regulating flavonoids and that phenylalanine metabolism further regulates proanthocyanidin biosynthesis [39–41]. Aradu.R8PMFand Araip.MHR6Kare components of MBW complex which are widely reported function in flavonoids metabolism (Li, 2014). The information that these genes cause a white seed color phenotype require further confirmation in future functional genomics studies.