HPLC-grade acetonitrile and methanol were purchased from Merck (Darmstadt, Germany). MilliQ water (Millipore, Bradford, USA) was used in all experiments. Formic acid and the standards were purchased from Sigma-Aldrich and MCE (MedChemExpress, USA). Four peanut cultivars with different testa color were used as materials in this study. “Zhonghua 9” (ZH9) is an elite peanut cultivar with black testa, and “Kangqibaihong” (KQBH) has white seed coat. The testa color of “HongHong” (HH) cultivar is red, while “Huapi” (HP) shows red with white spot in its seed coat. Peanut materials were planted in the test field of Oil Crop Research Institute (Wuchang), Chinese Academy of Agricultural Sciences, Wuhan, China. Seeds with similar size of the four cultivars were collected at different stages of development (from R4 to R8) (Boote, 1982). For the metabolomic and RNA-Seq (as well as qRT-PCR) experiments, testa was manually peeled from the seeds at the final mature stage (R8), immediately frozen in liquid nitrogen, and stored into -80°C refrigerator until use for analysis.

Total flavonoid content in peanut seed coat was determined by the Plant Flavonoid Detection Kit from Beijing Solarbio Science & Technology Co., Ltd. (Beijing, China). Briefly, testa samples derived from R4 to R8 developmental stages were ground into powders in liquid nitrogen, then about 0.1 g powders were extracted overnight at 4°C with 1.0 ml 70% aqueous methanol. The mixtures were centrifuged at 12,000 rpm for 10 mins, and the absorbance of supernatants was measured at 470 nm. Flavonoid content was qualified by using a calibration curve of rutin with linearity range from 0.04 to 2.5 mg/mL (R2 > 0.99) as described in the manual protocol. Total anthocyanin content in peanut testa was assayed using Plant Anthocyanidin Determination Kit from Shanghai yuanye Bio-Technology Co., Ltd (Beijing, China) according to user’s manual. In brief, testa samples were ground into powders in liquid nitrogen, then about 0.1 g powders were mixed with 2 mL Anthocyanidin Assay Buffer (75% methanol, 25% water, and 0.5% acetic acid). The suspensions were vigorously vortexed for 5 min, sonicated for 30 min, and then placed at 4°C overnight. The mixtures were centrifuged at 12000 rpm for 10 mins. The absorbance of supernatants was read at 530 nm.

The freeze-dried samples were crushed using a mixer mill MM 400 (Retsch, Germany) with a zirconia bead for 1.5 min at 30 Hz. 20 mg powders for each sample (with three biological replicates) were weighted and extracted with 0.5 mL 70% aqueous methanol. The extractions were sonicated for 30 min. After centrifugation at 12, 000 rpm under 4°C for 10 min, the supernatants were collected, and filtrated with a membrane (0.22 μm; Anpel) before UPLC-MS/MS analysis. An AB ExionLC™ Ultra-Performance Liquid Chromatography (UPLC) coupled with a QTRAP^® 6500+ MS System (AB, SCIEX, Framingham, MA, USA) was used for the quantification of target flavonoid metabolites. Separations were carried out using a Waters ACQUITY UPLC HSS T3 C18 column (100 mm×2.1 mm, pore size 1.8 µm) (Waters Corporation, Milford, USA). The mobile phase consisted of eluent A (water, 0.05% formic acid) and eluent B (acetonitrile, 0.05% acetic acid) at a flow rate of 0.35 mL/min with the temperature set at 40°C. The gradient elution program was set as follows: 0-1 min, 10%-20% B; 1-9 min, 20%-70% B; 9-12.5 min, 70%-95% B;12.5-13.5 min, 95% B; 13.5-13.6 min, 95%-10% B;13.6-15 min, 10% B. The injection volume is 2 μL. The effluent was analyzed by an ESI-triple quadrupole-linear ion trap-MS (ESI-Q TRAP-MS/MS). LIT (linear ion trap) and triple quadrupole scans in positive or negative ion mode were acquired on the LC/MS2 system, which was equipped with an ESI Turbo Ion-Spray interface. The optimized ESI parameters were as follows: ion source, turbo spray; temperature, 550°C; and ion spray (IS) voltage, 5500 V (Positive), -4500 V (Negative). Multiple reaction monitoring (MRM) experiments with the collision gas (nitrogen) of 35.0 kPa were performed to acquire QQQ scans. De-clustering potential (DP) and collision energy (CE) used for individual MRM transitions were assessed with further optimization.

Flavonoids were identified and annotated using a self-built database MWDB (MetWare biological science and Technology Co., Ltd. Wuhan, China) as well as the publicly available metabolite databases (Chen et al., 2013) according to their retention time (RT), MS2 spectrums obtaining from the MRM mode of QQQ mass spectrometry. The stock solutions of standards were prepared at the concentration of 5 μmol/L in 70% MeOH, and then diluted to a series of working solutions with gradient concentrations. The standard curves were established between the concentrations and their corresponding chromatographic peak areas obtaining from specific precursor ions in the MS2 spectra. The absolute contents of metabolites were calculated based on linear regression equations obtaining from standard curves ( Table S1 ). Since a part of anthocyanin compounds were not commercially available so far, a calibration curve of delphinidin-3,5-O-diglucoside was used to calculate the concentrations of these anthocyanins. The experimental data were obtained from three independent biological replicates. For multivariate statistical analysis, Partial Least Squares-Discriminant Analysis (PLS-DA) was used to identified the alternation of metabolites among different groups. An absolute value of log2 (fold change) ≥ 1 and variable importance in the projection (VIP) > 1.0 were introduced to screen significant differential metabolites (SDMs).

Total RNA was extracted from peanut seed coats collected at the final mature stage (R8). RNA concentrations and purities were measured using NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE). 1 μg of RNA per sample was used for preparing transcriptome library. The cDNA libraries were sequenced and constructed using the Illumina sequencing platform. Clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N and low-quality reads from raw data. The clean paired-end reads were then mapped to the peanut reference genome (Arachis hypogaea cv. Tifrunner, version 2, https://www.peanutbase.org/peanut_genome) using HISAT2 software (version 2.0.5) (Kim et al., 2019). The multi-mapped reads were filtered out, and only uniquely mapped reads were used for quantification. The gene expression levels were estimated by FPKM values (fragments per kilobase of transcript per million fragments mapped). Differentially expressed genes (DEGs) were screened using the DESeq2 R package (1.20.0) with the criteria of |log2FC (fold change) | > 1 and an adjusted P-value < 0.01 (Love et al., 2014). KEGG enrichment analysis of DEGs was performed by KOBAS software (Mao et al., 2005).

Total RNA was isolated from peanut seed coats using an EASYspin plant RNA extraction kit (Aidlab Biotechlogies, Co., Ltd, Beijing, China) according to the manufacturer’s instructions. Genomic DNA was removed from RNA samples by use of DNase I, then about 1 μg RNA were reverse transcribed into first strand cDNA with MMLV reverse transcriptase kit (Thermofisher Scientific, USA). qRT-PCR reactions were analyzed on Bio-Rad CFX96 RT-PCR Detection system (Bio-Rad, Hercules, CA, USA) using Hieff qPCR SYBR Green Master Mix (YEASEN, Shanghai, China). The relative transcript levels were quantified by the 2^-ΔΔCt method (Livak and Schmittgen, 2001), using A. hypogaea Actin gene (accession number: Aradu.W2Y55) as an internal standard gene for normalization. qRT-PCR primers were listed in Table S2 .

A weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA R package (v1.68) (Langfelder and Horvath, 2008) based on FPKM values of the filtered DEGs. Network construction and module identification were conducted using the topological overlap measure (TOM) with following parameters: “weighted network” = unsigned, “soft thresholding power” =15, “minimum module size” = 30, “minimum height for merging modules” = 0.25. Then, the calculated module eigengenes and Pearson’s correlation coefficient values were used to determine the association of modules with the major metabolite contents for the 12 samples. The number of edges in a node represented the hubness of the gene. The gene networks and top 20 hub genes within a module were visualized by Cytoscape software (Shannon et al., 2003).

In general, peanut seed development can be divided into nine reproductive stages (R1-R9) (Boote, 1982). Four cultivar seeds with typical testa color were collected from full pod stage (R4) to harvest maturity stage (R8), respectively ( Figure 1A ). Testa of ZH9 is black, and KQBH is white without pigment accumulation during seed maturation. The seed coat color of Honghong (HH) and Huapi (HP) are both light pink at the full pod stage (R4), and gradually get deeper during seed development. At the harvest maturity stage (R8), the seed coat of HH (Red) completely turns to red, while that of HP (Spot) is dark red with some white spots.

ZH9 with black testa had a large amount of total anthocyanin as early as the full pod stage (R4), whose contents were significantly higher than those of red (HH) and/or spot (HP) cultivars during the whole developmental stages ( Figure 1B ). In contrast, a trace amount of anthocyanin was detected in white seed coat (KQBH), which is in accordance with its phenotype. Moreover, TFCs in the four cultivars increased from full pod stage (R4), and reached to their maximum levels at the stage of R6 or R7. At harvest maturity stage (R8), TFCs in red (HH) and/or spot (HP) peanut skin were much higher than those in black (ZH9) and/or white (KQBH) cultivars ( Figure 1C ). These results suggest that peanut testa color is associated with anthocyanin as well as other flavonoid levels, which conforms with previous studies (Shem-Tov et al., 2012; Attree et al., 2015).

Early investigations have shown type-A procyanidins predominated in peanut skins, while type-B were commonly found in important dietary sources such as apples and grapes (Zhang et al., 2013; Dudek et al., 2017; Rue et al., 2018). In current study, the amounts of type-A procyanidins were less than those of type-B in the red and spot cultivars, but type-A overwhelmed type-B in black peanut. These results indicated that the regulation mechanisms of procyanidin biosynthesis were distinct among different cultivars, and the differential accumulation pattern of procyanidins was possibly related to their color formation. Comparing the results for the different analyzed peanut skin samples (except for white testa), it is notable that the content ratio of procyanidin oligomers to their monomers (flavan-3-ols) displayed more or less the same levels although the absolute concentrations differed ( Table 1 and Figure S3 ). This aspect has also been observed in some other foods and nuts (De Pascual-Teresa et al, 2000; Rzeppa et al., 2011). Four beans belonging to the same genus in legumes, such as black bean and kidney bean, had different concentrations for the procyanidin and monomeric flavan-3-ol, as well as for the total amounts, but their relative compositions displayed a similar pattern (Bittner et al., 2013). The result suggests that the content ratio of procyanidins to monomeric flavan-3-ols might be characteristic for the different sample types and can be used for authenticity control, while their absolute concentrations varied possibly due to the origin of samples or stage of development.

In our transcriptome data, several MYB-like transcription factors (eg. MYB30, MYB80, MYB86, MYB102) and an ANR (arahy.IK60LM) were found to be hub genes involved in anthocyanin and proanthocyanin metabolism by WGCNA. Previous study showed that R2R3-MYB transcription factor and ANR were important regulatory genes that confer pigment accumulation in plant tissues (Yang et al., 2022). A major gene (AhTc1) controlling peanut purple testa was characterized to be encoding a R2R3-MYB transcription factor using QTL-seq method together with transcriptome analysis and gene expression profiling. Overexpression of AhTc1 (arahy.J3K16L) in tobacco led to different degrees of purple color in leaves, flowers, fruits and testa (Zhao et al., 2020). Recently, a major location (AhRt2) controlling peanut testa color was fine-mapped to a 0.5 Mb genomic region on chromosome 12 using the BSA−seq and linkage mapping approaches, and the ANR (arahy.IK60LM) was suggested to be the possible candidate gene within this region (Zhang et al., 2022), which was in accordance with our WGCNA analysis.

In the tailing process, anthocyanidins are converted into colored pigments via glycosylation catalyzed by UGTs. Compared with other cultivars, black peanut had substantially higher level of anthocyanidin 3-O-glycoside, which is a main stable colored pigment in several plant tissues. Based on the transcriptome data, four putative UGTs and three ANRs were identified to be involved in testa pigmentation. The UGT transcripts have relative higher expression level in black peanut. It is proposed that UGTs and ANR compete for the substrate cyanidin ( Figure 4C ), and the prevalence of UGTs activities over ANR one will determine the color pattern of peanut testa. Moreover, cyanidin 3-O-sophoroside and cyanidin 3-O-sambubioside are the two major anthocyanins in black peanut skin, indicating glycosyl extension from anthocyanidin 3-O-glucosides plays a vital role in testa pigment formation (Zhao et al., 2017; Huang et al., 2019). Two specific UGTs (Ib3GGT and At3GGT) were characterized in Arabidopsis and purple sweet potato. The two enzymes were responsible for catalyzing the conversion of anthocyanidin 3-O-glucosides to anthocyanidin 3-O-sophorosides or anthocyanidin 3-O-sambubioside using UDP-glucose or UDP-xyloside as a sugar donor, respectively (Yonekura-Sakakibara et al., 2012; Wang et al., 2018b). Five peanut UGTs were found to have high sequence similarities with Ib3GGT and At3GGT, and putatively annotated as anthocyanidin 3-O-glucoside 2’’-O-glycosyltransferase ( Figure S5A ). Phylogenetic analysis showed that AhUGT236 (arahy.BNFJ8I) exhibited the most closed relationship with Ib3GGT, and its transcript level was higher in black testa than red and white samples ( Figures S5B, C ). These results indicated AhUGT236 was possibly a candidate gene involved in the biosynthesis of anthocyanidin 3-O-sophoroside, and might play vital roles in testa pigmentation. However, the gene function and its regulatory mechanism warrant further investigation in future studies.