AL12 (almost no anthocyanidin in the aboveground organs) and ZL23 (ZL23, anthocyanidin rich in leaves, petioles and hypocotyls) varieties of mung beans were used as plant materials whose seeds were generously provided by the Miscellaneous Grains Research Laboratory of Henan University of Science and Technology. The experiment was carried out from May to September 2022. All experimental procedures on plant materials were conducted at the experimental farm of Henan University of Science and Technology, situated at 34° 41′ N and 112° 27′ E, at an altitude of 280 m above sea level. The mean annual rainfall and temperature in the region is 578.2 mm and 14.8°C, respectively. The soil in the experimental farm is classified as fluvo-aquic soil, and has the following physicochemical characteristics: pH (H₂O) = 7.08, available nitrogen = 80.09 mg·kg^–1, available phosphorus = 3.31 mg·kg^–1, available potassium = 81.32 mg·kg^–1, organic matter = 14.5 g·kg^–1. Cultivation, and cultivation management was kept the same as for the general high-yielding fields. The leaves, petioles, and hypocotyls of both varieties were sampled after the 5-leaves stages. GL, GB, and GP stand for leaves, petioles, and hypocotyls of AL12, respectively; BL, BB, and BP stand for leaves, petioles, and hypocotyls of ZL23, respectively. Samples were snap frozen in liquid nitrogen and stored at -80°C in a refrigerator till further use. In order to ensure the reliability and reproducibility of the experiment, physiological indexes, metabolome analysis and transcriptome analysis both use 3 biological replicates.

The total anthocyanin contents were determined using the hydrochloric acid (HCl)-ethanol method (Wang et al., 2013). Briefly, an accurately weighed 0.1g of sample was taken and ground using a tissue grinder in liquid nitrogen into powder and extracted at 60°C in a water bath for 30 min with 10 ml of 0.1 mol/L HCl-ethanolic solution. The supernatant was carefully separated and transferred to a new 25 ml volumetric flask and extracted again as per the previous procedure. Both the extract fractions were combined, and the final volume was made up to 25 ml using HCl-ethanolic solution. The absorbance of the extract was then recorded at wavelengths of 530 nm, 620nm, and 650nm, respectively. The anthocyanin contents were then estimated using the following equation, and worked out their content by molecular weight of cyanidin-3-O-galactoside.

Where OD is optical density.

Where V is total volume of extract; ε is molar extinction coefficient of anthocyanin; m is sample weight; M is molecular weight of cyanidin-3-O-galactoside (449 g/mol).

Anthocyanin metabolome analysis was completed by Wuhan Metware Biotechnology Co., Ltd. through LC-MS/MS. Lyophilized plant tissue sample was ground into powder by a tissue grinder (MM 400, Retsch). 50 mg sample was dissolved in 500 μL extract solution (50% methanol containing 0.1% HCl). The extract solution was oscillated by vortices for 5 min, and was extracted by ultrasonic extractor for 5 min. Then, the extract solution was centrifuged for 10 min (4°C, 12000 rpm min^-1). The supernatant was transferred to a new centrifuge tube and repeat the extraction process. The extract was filtered with a microporous filter membrane (0.22μm pore size). The extract was then analyzed using a liquid chromatography tandem mass spectrometry system (Applied Biosystems 6500 QTRAP). The liquid chromatographic conditions were as follows: the column was C18 (ACQUITY BEH, 1.7μm, 2.1mm×100mm), and the mobile phase A was ultrapure water (containing 0.1% formic acid), B was methanol (containing 0.1% formic acid), the flow rate was 0.35 mL min^-1, the column temperature was 40°C, and the injection was 2 μL. Mass spectrometry conditions were as follows: electrospray ion source temperature was 550°C, mass spectrometry voltage was 5500V, curtain gas was 35psi. The acquired mass spectrometry data were then qualitatively analyzed by means of a local database. Finally, the anthocyanin components were quantitatively analyzed using the multiple reaction monitoring mode of triple quadrupole mass spectrometry.

Total RNA was extracted from samples using TaKaRa MiniBEST Plant RNA Extraction Kit (TAKARA 9769), and cDNA libraries were constructed for samples meeting the purity, concentration, and integrity of the RNA. The mRNA was enriched with magnetic beads with Oligo (dT), and randomly interrupted by adding a fragmentation buffer. The first cDNA strand was synthesized with random hexamers using the mRNA as a template, followed by purification through AMPure XP beads. The purified double-stranded cDNA was then end-repaired, A-tailed, and ligated to sequencing junctions, followed by fragment size selection using AMPure XP beads, and finally enriched by PCR to obtain the cDNA library. The cDNA libraries were sequenced by Biomarker Technologies Ltd. (Beijing, China) using the Illumina platform. The adapter sequences and low-quality sequences of the off-machine data were filtered to obtain clean data. Sequence alignment was performed with the reference genome of mung bean, and uniq mapped genes were obtained through HISAT2 and StringTie software. Functional annotation of a single gene was performed using multiple databases including NR (Non-Redundant Protein Sequence Database), Swiss-Prot (Swiss Protein Sequence Database), GO (Gene Gntology Annotation), COG (Clusters of Orthologous Groups), KOG (EuKaryotic Orthologous Groups), Pfam (Protein Family Analysis and Modeling), and KEGG (Kyoto Encyclopedia of Genes and Genomes).

The expression levels of all genes were calculated and normalized to fragments per kilobase of transcript per million fragments mapped (FPKM). Fold Change≥1 and FDR<0.05 were used as the criteria for screening differential expression genes (DEGs). The identified DEGs were further enriched by GO, COG, and KEGG pathway analysis. Gene co-expression network analysis was constructed using the R language software WGCNA package. Since input gene expression was the result of TPM normalization, the expression was first subjected to a further log2 transformation. To ensure that the data conform to the scale-free network distribution, the weighting coefficient β value was screened and the network topology was analyzed using the pickSoftThreshold function in the WGCNA package to obtain the correlation value (R²) corresponding to the alternative soft threshold, which is required to ensure R² is close to 0.8 having a certain degree of gene connectivity. A horizontal line was drawn with R^2 = 0.8 as the threshold, and the first-time power exceeding the threshold was found to be 15. Since power=15 was having a high degree of connectivity, the β=15 was selected to construct the co-expression network and for subsequent analysis. To quantify the co-expression similarity of the whole module, the dynamic tree cut algorithm was used to identify the co-expression modules of the transcriptome expression data, and the minimum number of variables allowed to be included in the module (minModuleSize) was set to 100. To quantify the similarity of co-expression across modules, the characteristic genes of the modules were calculated, with unmerged modules, and other parameters were set by default. The modules are plotted below the tree using the plot dendro and colors function. The soft connectivity function was used to calculate the connectivity between genes, and the top 10 genes with kME values as the core genes of the module (hubgene) were selected. The topological overlap matrix was used to calculate the weights of the different genes in the module, where a higher weight value indicated a higher degree of association between genes. Due to the large network size, a weight value > 0.35 was set as network filtering criteria, and the reciprocal network of 10 hubgene was visualized using Cytoscape v3.6.1. Based on the results of WGCNA, the dark green gene module with the highest correlation to anthocyanins was selected, and the genes within the module were further visualized and analyzed using Cytoscape v3.6.1 for the interactions network. The top 8 genes with the highest degree of linkage were selected as the core genes (hubgene) and were focused owing to key genes affecting anthocyanin biosynthesis.

Total RNA extraction, reverse transcription reactions, and qPCR were performed using the kit method (TaKaRa; 9769, RR047B, and RR430B), and primers were designed ( Table S1 ). A total of 20 µl reaction volume contained 10 µl of TB Green Fast qPCR Mix (2×), 1 µl cDNA, 0.4 µl each of upstream and downstream primers (10 μmol. L^-1), and 8.2 µl doubled distilled water (ddH₂O). The reaction proceeded in a two-step process i.e., pre-denaturation at 95°C for 30 seconds, followed by denaturation at 95°C for 5 seconds, and annealing at 60°C for 30 seconds. A total of 40 cycles were performed, and fluorescence was determined at reversion time. The experiment was conducted in triplicates and relative expression was calculated using the 2^{- ΔΔCT} method with VrActin as the internal reference gene.

The CDS sequences of MYB3 and MYB90 were cloned into pCAMBIA1300-GFP vector without termination. Recombinant plasmid of 35S:MYB3-GFP and 35S:MYB90-GFP were transformed into agrobacterium GV3101 for subcellular localization. Agrobacterium (35S:VrMYB3-GFP and 35S:VrMYB90-GFP) was identified as positive clones which was cultured with liquid LB medium (50 μg mL^-1 of kanamycin and rifampicin each) at 28°C. 1.5 mL of bacterial solution was centrifuged at 5 000 r min^-1 for 3 min to collect bacterial cells. Then, the OD₆₀₀ of agrobacterium suspension, including 1 mol L^-1 MES (pH 5.6), 1 mol L^-1 MgCl₂ and 20 mmol L^-1 acetosyringone, was adjusted to about 0.7. Agrobacterium suspension was standed still at 25°C for 1 h. The bacterial solution was slowly injected into the mesophyll cells on the back of Nicotiana benthamiana leaves, and the injected plants were cultured under light for 48 h. The subcellular localization of fusion protein was observed through LSM 710 NLO confocal microscope (Zeiss, Germany).

All the data are presented as mean ± SD (standard deviation). Microsoft Excel 2016 and SPSS17.0 were used to process the data while comparing the significance between treatments was done using the LSD method.

To investigate the biosynthesis and accumulation of anthocyanins during the growth and development of AL12 and ZL23 mung bean varieties, the phenotype and anthocyanin content of leaves, petioles, and hypocotyls of both varieties were analyzed. All three parts in the common green bean variety were green in color with no anthocyanin accumulation observed during growth and development. In contrast, the leaf veins, petioles, and hypocotyls of the greenish-black bean leaves showed significant anthocyanin accumulation during the same period ( Figures 1A-P ). The anthocyanin content in leaves, petioles, and hypocotyls during the same period were found to be significantly lower in AL12 compared to ZL23. A substantial difference in anthocyanin contents among the two varieties was found in petioles and hypocotyls, where AL12 contained 29.3 mg g^-1 in petioles and 4.5 mg g^-1 in hypocotyls, while the same was found to be 90.2 mg g^-1 and 75.2 mg g^-1 in ZL23, respectively, showing a 3.07-fold and 16.7-fold increased anthocyanin contents in ZL23. Nevertheless, the difference in anthocyanin contents in leaves among the two varieties was relatively small i.e., 10.3 mg g^-1 in AL12, and 17.4 mg g^-1 in ZL23, respectively, where the contents in ZL23 were still 1.7-fold higher compared to AL12 ( Figure 1Q ), which could be attributed to the fact that anthocyanin in leaves is mainly accumulated in the leaf veins.

Principal component analysis (PCA) showed that there were significant differences between the two varieties and three site samples ( Figure 1R ). 72.16% of the differences between samples could be explained by PCA1 (47.79%) and PCA2 (24.37%), indicating that anthocyanins showed a dynamic change pattern during two varieties and three site samples. Hierarchical cluster analysis (HCA) of metabolome data exhibited low difference in biological replicates. The correlation coefficient of anthocyanin components contents was up to 0.93 among biological replicates ( Figure 1S ). This suggests that there was a good correlation among the bio-replicates. The above results demonstrated the results was stable and repeatable, and provided a guarantee for the reliability of the results. Seven categories, consisted of 29 kinds of anthocyanin compounds, in these three stages were identified through LC-MS/MS, including 8 kinds of cyanidin, 7 kinds of delphinidin, 7 kinds of pelargonidin, 2 kinds of malvidin, 2 kinds of peonidin, 2 kinds of petunidin, and 1 kinds of procyanidin ( Figure 2 , Table S2 ). In AL12, cyanidin was the highest proportion of anthocyanins in the leaves and petioles, accounting for 84.1% and 55.9% respectively; pelargonidin was the highest proportion of anthocyanins in the hypocotyl, accounting for 82.8%. In ZL23, cyanidin was the highest proportion of anthocyanins in the leaves, accounting for 86.5%; delphinidin was the highest proportion of anthocyanins in the petioles and hypocotyl, accounting for 95.3% and 64.3%, respectively.

The PCA results of transcriptome data showed a clear separation between the two varieties and three site samples, relatively. A negligible variability in transcriptome expression was found in GB and GP of AL12, and BB and BP of ZL23, while the same was higher in the leaves (GL, and BL). The 47.4% of the variation between samples was explained by PCA1 (31.58%) and PCA2 (15.82%), indicating that anthocyanin accumulation showed significant differences between the two varieties and different parts ( Figure 3A ). Correlation coefficient analysis between biological replicates as shown in Figure 3B indicated that the correlation coefficient of gene expression levels among biological replicates of all samples was >0.91, indicating an excellent reproducibility among biological replicates and that the data could be further used to determine DEGs.

To further clarify the molecular mechanism of anthocyanin biosynthesis during the growth and development of stated organs in ZL23, 18 cDNA libraries comprising two varieties of three organs in triplicates were prepared for transcriptome analysis. As shown in Table S3 , a total of 116.64 Gb of clean data was obtained after filtering, with 19, 236, 210 to 38, 429, and 484 clean reads for each library. The percentage of Q30 (sequences with a sequencing error rate<0.1%) for each library was over 92%, while GC contents of the samples ranged between 44.67 and 46.20%. The clean reads were mapped between 83.58 and 93.87% to the reference genome of each library and which was then used for further analysis.

A total of 8211 DEGs were screened in the three comparison combinations as shown in Figure 3C , among which, a total of 4921 were up-regulated and 3284 were down-regulated. Additionally, 2803 DEGs were screened in the GB-1_GB-2_GB_3_vs_BB-1_BB_2_BB_3 comparison sample, among which 1882 were up-regulated and 981 were down-regulated. The GL-1_GL-2_GL_3_vs_BL-1_BL_2_BL_3 comparison sample screened a total of 3548 DEGs, with 2007 up-regulated and 1541 down-regulated genes, while GP-1_GP-2_GP_3_vs_BP-1_BP_2_BP_3 comparison sample screened a total of 1860 DEGs, where 1098 were up-regulated and 762 were down-regulated. These results suggested that the DEGs in AL12 and ZL23 might be involved in the synthesis of anthocyanins in mung beans. As shown in the Wayne diagram ( Figure 3D ), the number of DEGs specific to leaves, petioles, and hypocotyls were 2053, 1167, and 535, respectively, where the number of common DEGs between leaves and petioles was 567, and between leaves and hypocotyls was 256. Similarly, the number of common DEGs between petioles and hypocotyls was 397, while 672 were among leaves, petioles, and hypocotyls. Further, 833 transcription factors were found to be differentially expressed in the two varieties ( Figures 3E-F , Table S4 ), belonging to 15 transcription factor families, among which 109 and 92 were significantly up- and down-regulated, respectively.

The GO enrichment analysis results are shown in Figure S1 and Table S5 , where 52 functional components were grouped into 3 categories, including ‘biological processes’, ‘molecular functions’ and ‘cellular components’. The highest number of annotations was shown for ‘metabolic processes’ within ‘biological processes’, followed by ‘cellular processes’. In ‘cellular components’, the most numerous genes were annotated as ‘cellular’, ‘cellular components’, and ‘membrane components’. Similarly, in ‘molecular function’, the most genes were annotated as ‘binding’ and ‘catalytic activity’. These results indicated that many types of enzymatic pathways were active in mung beans. To further identify the metabolic pathways of DEGs, KEGG enrichment analysis was performed. As shown in Figure S2 ; Table S6 , the top 20 KEGG pathways with the lowest Q values were mainly concentrated in ‘phytopathogen interaction’, ‘phenylpropane biosynthesis’, ‘photosynthesis’, ‘flavonoid biosynthesis’, and ‘flavonoid and flavonol biosynthesis’ pathways. The ‘flavonoid biosynthesis’ pathway was significantly enriched in all three comparative combinations, and results suggested that anthocyanin biosynthesis belonged to the flavonoid pathway and was also associated with the phenylalanine metabolic pathway. The large enrichment of the flavonoid biosynthesis pathway correlated with the large amount of anthocyanin synthesis in ZL23.

The published literature suggests that anthocyanins are mainly synthesized by 12 structural genes i.e., CHI, CHS, C4H, CHS, 4CL, DFR, F3’H, F3’5’H, F3H, LDOX, PAL, UFGT, and three major transcriptional factors namely MYB, bHLH, and WD40, which were made into HCA ( Figures 4 , 5 ). As shown by HCA results, the reliability of the anthocyanin content and PCA results was further illustrated by the fact that the genes related to anthocyanin synthesis were significantly different between the two varieties. The HCA of structural genes is shown in Figure 5A , which were at significantly higher expression levels in the petioles and hypocotyls of ZL23 compared to AL12. However, the differences in expression in leaves were relatively small. Moreover, the three regulatory genes i.e., MYB, WD40, and bHLH showed similar results as shown in Figures 5B-D and Table S4 .

To further characterize the expression of genes related to anthocyanin biosynthesis, the relative expression of 12 structural genes for anthocyanin biosynthesis and MYB, bHLH, and WD40 were determined by RT-qPCR ( Figure 6 ). The expression profile results indicated a high correlation between RNA-seq and RT-qPCR, thereby further cementing the reliability of the RNA-seq results.

To identify co-expressed gene modules and explore the relationship between gene networks and anthocyanin synthesis, the WGCNA was further performed, including 6415 DEGs. The correspondence between correlation coefficients and average connectivity at different thresholds (power values = 1 to 30) was measured. As shown in Figure 7 , correlation coefficients, and average connectivity were guaranteed when the power value was 30. Finally, 14 modules were obtained from the clustering dendrogram ( Figure 7D ), where each module contained a different number of genes, with the ‘BROWN’ module containing the most genes at 1511 and the ‘BLUE’ module containing the least genes at 49, with an average of 458 genes per module. Pearson correlation analysis was performed between the genes and traits in each module according to the amount of anthocyanins in each site, and results are shown in Figure 7 . The results demonstrated that seven of the 14 modules were positively correlated and seven were negatively correlated with anthocyanin synthesis. The green and dark-grey modules had the highest positive and negative correlation with anthocyanin content, respectively. The top 8 genes in the green module were selected as the core genes, including Vigna_radiata_var._radiata_newGene_3976, Vigna_radiata_var._radiata_newGene_4853, gene10693, gene25494, gene6451, gene12756, gene22866, and gene21475. These genes and their associated genes were mapped in a visual gene interaction network ( Figure 8 ). Results showed that gene22866 and gene21475, which encode LDOX and F3’5’H respectively, were structural genes for anthocyanin synthesis, while Vigna_radiata_var._radiata_newGene_3976, Vigna_radiata_var._radiata_newGene_4853, gene10693, gene25494, gene6451, gene12756 belonged to MYB transcription factors. Due to the low expression of Vigna_radiata_var._radiata_newGene_3976, Vigna_radiata_var._radiata_newGene_4853, gene10693, gene25494 (FPKM ≤10), they will be ignored. Based on the combined analysis of the above results, gene12756 (VrMYB3) and gene6451 (VrMYB90) were found to play important roles in the regulation of anthocyanin synthesis in mung beans.

To confirm the location of VrMYB3 and VrMYB90 proteins in cells, the transiently expressed fusion proteins of VrMYB3 and VrMYB90 in tobacco leaves were detected by confocal laser microscopy. Figure 9 showed that there were GFP signals on both cell nucleus and cell membrane after transiently expressed 35S:GFP vector in tobacco epidermal cells. However, transiently expressed of 35S:VrMYB3-GFP and 35S:VrMYB30-GFP vectors in tobacco epidermal cells, GFP signals were present on the nucleus. It suggests that VrMYB3 and VrMYB90 have potential function of transcriptional regulatory.