Transcriptome Profiling During Muscadine Berry Development Reveals the Dynamic of Polyphenols Metabolism

Muscadine grapes accumulate higher amounts of bioactive phenolics compared with other grape species. To identify the molecular events associated with polyphenolic accumulation that influence antioxidant capacity, two contrasting muscadine genotypes (C5 and C6) with varied phenolic/flavonoid content and antioxidant activity were investigated via RNA-sequencing during berry development. The results showed that berry development is concomitant with transcriptome profile changes, which was more pronounced at the véraison (V) stage. Despite that the downregulation pattern of gene expression dominated the upregulation through berry development, the C5 genotype maintained higher expression levels. Comparative transcript profiling allowed the identification of 94 differentially expressed genes with potential relevance in regulating fruit secondary metabolism, including 18 transcription factors and 76 structural genes. The genes underlying the critical enzymes in the modification reactions of polyphenolics biosynthetic pathway, including hydroxylation, methylation, and glycosylation were more pronounced during the immature stages of prevéraison (PrV), V, and postvéraison (PoV) in the C5 genotype, resulting in more accumulation of biologically active phenolic/flavonoid derivatives. The results suggested that muscadine grapes, as in bunch grapes (Vitis sp.); possess a similar mechanism that organizes polyphenolics accumulation; however, the set of total flavonoids (TFs) and structural genes coordinating the pathway varies between the two species.


INTRODUCTION
Grapevine is a member of the large genus Vitis that is divided into two subgenera with a varying number of somatic chromosomes (Olien, 1990). The Euvitis (bunch grapes) displays a large group that is classified based on the geographical origin into three groups, including the European (V. vinifera), the East Asian (V. amurensis), and the North American species (V. riparia, V. labrusca, etc.) (Šikuten et al., 2020). However, Muscadinia is a smaller group that comprises three species, including M. rotundifolia, M. munsoniana, and M. popenoei (Olien, 1990). The muscadine grapes (M. rotundifolia), unlike other Muscadinia species, are commercially valuable growing mainly for fresh fruit and wine production but with lesser popularity than the European grapes (Andersen et al., 2021). However, muscadines are receiving growing attention due to their tolerance to several diseases that cause extensive economic losses in bunch grapes (Staudt, 1997;Merdinoglu et al., 2003). Additionally, muscadine berries display enhanced nutraceutical value due to the accumulation of distinctive phytochemical constituents that have a great potential activity against several chronic diseases (Bralley et al., 2007;God et al., 2007;Mellen et al., 2010;Gourineni et al., 2012;Luo et al., 2017;Mendonca et al., 2019). Such nutritional and health merits are not only restricted to muscadine, but muscadine berries accumulate higher amounts of bioactive polyphenolics compared with other grape species (Pastrana-Bonilla et al., 2003;You et al., 2012;Mendonca et al., 2019).
Polyphenolic compounds are a large group of plants' secondary metabolites synthesized in multiple sequential pathways, including the pentose phosphate, shikimate, and phenylpropanoid pathways (Randhir et al., 2004). They have substantial functions in plants, including growth, development, and defense. Polyphenolics also define berry quality and sensory attributes, such as texture, flavor, and pigmentation properties (Gomez-Plaza et al., 2001;Duan et al., 2019). The composition and amount of phenolics are highly coordinated by many factors, including genotypic variations, environmental conditions, and more importantly, developmental stages (Darwish et al., 2021). Grape berry development exhibits a double-sigmoid growth pattern, including (I) berry formation, (II) lag phase, and finally (III) ripening phase (Coombe and McCarthy, 2000). Stages I and III are well-known for their exponential increase in berry size, a missing feature in the lag phase. During immature stages (I and II), the accumulation of organic acids, mainly malic and tartaric acids, raises berry acidity levels (Sweetman et al., 2012). However, the acidity declines when the critical point known as véraison occurs at the end of stage II. Furthermore, the berry begins accumulating sugars (mainly glucose and fructose) and anthocyanins in colored varieties. These physiological and biochemical changes proceed during the ripening stage that exhibits a striking reduction in the organic acid coincided with a considerable accumulation of sugar and pigments (Dai et al., 2013). Among the distinctive characteristics that distinguish each phase, the accumulation of secondary metabolites is of great interest with more than 300 individual phenolic and flavonoid compounds reported being abundant in grapes and wine (Šikuten et al., 2020;Darwish et al., 2021).
The economic value of grapes draws significant attention to manifest all grape-related aspects, emphasizing berry development. Since the genome of grapevine (Vitis vinifera) has been sequenced and annotated (Jaillon et al., 2007), several studies have been conducted to better elucidate the transcriptomic and metabolic reprogramming that govern berry development (Deluc et al., 2007;Nwafor et al., 2014;Palumbo et al., 2014;Degu et al., 2015). In contrast, muscadine received less recognition with relatively few high-throughput studies (Pinu, 2018;Duan et al., 2019). The recent release of muscadine whole genome sequence has facilitated transcriptomic studies associated with several berry quality aspects (Park et al., 2021). In our previous study, we utilized an untargeted metabolomics approach along with multivariate analysis to categorize the metabolome profile throughout the six berry developmental stages of three muscadine grape genotypes: the bronze cultivar Late Fry (LF) and two-colored breeding lines, C5-9-1 (C5) and C6-10-1 (C6) (Darwish et al., 2021). The data illustrated that the differences in total phenolic content (TPC) and total flavonoid content (TFC) levels resulted in alterations in DPPH, FRAP, and ABTS antioxidant activities. However, the véraison (V), postvéraison (PoV), and ripening developmental stages were the most distinguishing stages between genotypes. Given the changes in TPC and TFC levels during development, the following question was raised: what is the crucial gene network that can cause a difference in phenolic/flavonoid accumulation, and hence the antioxidant activity of muscadine grape. To answer this question, two contrasting muscadine genotypes (C5 and C6) with varied phenolic/flavonoid content and antioxidant activity were investigated. The changes in TPC/TFC levels and antioxidant activity during berry development along with RNA-sequencing (RNA-seq) strategy were utilized to identify the molecular events associated with phenolic/flavonoid accumulation. The results showed that berry development is concomitant with modulation of transcriptomic profile changes that were more pronounced at the V stage. Despite that, the downregulation model of gene expression prevailed over the upregulation pattern during the progression in berry development, the expression levels were highly pronounced in C5 during prevéraison (PrV), V, and PoV stages. Of paramount significance, the genes encoded key enzymes involved in the modification reactions of polyphenolics biosynthetic pathway, including hydroxylation, methylation, and glycosylation were more abundant during the V stage of C5 genotype than C6, defining the high accumulation of bioactive phenolic/flavonoid derivatives in C5.

Plant Materials
The berry samples were collected from 5-year-old muscadine grapevine (Muscadine rotundifolia (Michx.) small) genotypes C5-9-1 (C5) and C6-10-1 (C6), grown at the experimental vineyard of the Florida A&M University (Tallahassee, FL, United States). These two breeding lines were developed under the grape breeding program of the Center for Viticulture and Small Fruit Research. They were selected according to their diversity in TPC, TFC, and antioxidant activity (Mendonca et al., 2019;Darwish et al., 2021). The berry samples were collected at different developmental stages, including fruit-set (FS), PrV, V, PoV, and ripening ( Figure 1A and Supplementary Table 1). The berries were carefully separated into two different tissues at the ripening stage, designated as ripe skin/flesh (R) and ripe seeds (S) tissues. Five clusters/replicate and three replicates/genotype were randomly collected for all developmental stages, excluding the FS stage. At FS, 70 clusters/replicates were collected due to the small berry size. All samples were immediately frozen in liquid nitrogen and stored at -80 • C for further analysis.

Nucleic Acid Extraction and RNA-seq Library Construction
The total RNA from the muscadine berry tissues was extracted as described previously (Gambino et al., 2008). All RNA extracts were treated with the RNase-Free DNase Set (Qiagen, Valencia, CA, United States), and then cleaned up with the RNeasy Mini Kit (Qiagen, Valencia, CA, United States). A total of 24 RNA-seq libraries (three biological replicates at four stages, FS, V, R, and S) from the C5 and C6 muscadine genotypes were constructed as described previously (Bai et al., 2014)

RNA-seq Data Preprocessing and Identification of Differentially Expressed Genes
The Illumina sequencing of the multiplexed RNA-seq libraries yielded 24 FASTQ files of sequences. Reads quality was checked by FASTQ 1 two times before and after trimming using Trimmomatic v0.39 (Supplementary Table 2; Bolger et al., 2014). The trimmed reads were then subjected to Salmon in non-alignment-based mode to estimate transcript quantification (Patro et al., 2017). All the samples mapped to the novel generated transcriptome of muscadine grape yielded a mapping rate higher than 77.1% (Supplementary Table 2; Park et al., 2021). Differentially expressed genes (DEGs) during berry development were identified between consecutive developmental stages [V-FS, R-V, or seeds -R (S-R)] within the C5 or C6 genotype, and between analogous stages of C5 and C6, using DESeq2 and EdgeR package setting on default parameters (Love et al., 2014;Chen et al., 2016). For convenience, the common and unique DEGs of each comparison generated by DESeq2 and EdgeR pipeline (DEGs, P FDR < 0.05, log2fold change > 1.5 or < -1.5) were considered to be expressed (Supplementary Figures 3-5 and Supplementary Tables 3-5). The web-based tool Venny was used to construct the consensus result (Oliveros, 2007). Finally, K-means clustering was performed with the Hartigan-Wong Algorithm and based on inputs of log2 relative Transcripts Per Million (TPM) values.

Weighted Gene Co-expression Network Analysis
The co-expression network modules were constructed using the variance stabilizing transformation values and the R package WGCNA (v1.69; Langfelder and Horvath, 2008). Before analyzing the data, lowly expressed genes among all sample types were removed by DESeq2, and the remaining 20,886 genes were used in module construction. The co-expression modules were 1 https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ obtained using the default settings, except that the soft threshold power was 10, TOMType was unsigned, minModuleSize was 30, mergeCutHeight was 0.25, and scale-free topology fit index was 0.8 (R 2 = 0.8). A module eigengene (ME) value, which summarizes the expression profile of a given module as the first principal component, was calculated and used to evaluate the association of modules with berries biochemical property [TPC, TFC, and total anthocyanin content (TAC)] and antioxidant activity [DPPH, FRAP, ABTS, NORS, and CUPRAC] of C5 and C6 genotypes at different developmental stages (Darwish et al., 2021). As a result, the final matrix for WGCNA contained a total of 20,886 genes that were assigned to sixteen modules (M1-M16). In addition, the number of non-redundant DEGs from C5 stages , C6 stages , and C5 stage -C6 stage comparisons, which are present in each module, were determined.

Gene Ontology Enrichment and Kyoto Encyclopedia of Genes and Genomes Pathway Analyses
The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were assigned using the g:Profiler website by applying the Benjamini-Hochberg multiple testing correction method with P FDR < 0.05 (Raudvere et al., 2019). However, since the gene ID of muscadine transcriptome is not supported, the Ensembl gene ID of Vitis vinifera was used instead. The Cytoscape plug-in ClueGO was used to visualize the non-redundant BP terms and KEGG pathways for DEGs located in modules ME1, ME5, ME11, and ME14, as well as the 94 genes of interest (Bindea et al., 2009).

Validation of Differentially Expressed Genes Subsets by Quantitative PCR
The DNase-treated RNA (5 µg) was reverse transcribed in a reaction of 50 µl using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, United States). The gene-specific primers were designed using Primer Express (v3.0, Applied Biosystems, Foster City, CA, United States) (Supplementary Table 11). The real-time quantitative PCR (qPCR) assays were performed using 20 ng of cDNA and 300 nM of each primer in a 10-µl reaction volume with SsoAdvanced Universal SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA, United States). Three biological and three technical replicates for each reaction were analyzed on a CFX384 Touch Real-Time PCR Detection System instrument (Bio-Rad Laboratories, Hercules, CA, United States) with the first step of 95 • C for 5 min followed by 40 cycles of 95 • C for 10 s, 60 • C for 10 s, and 72 • C for 20 s. Melting curves were generated using the following program: 95 • C for 15 s, 60 • C for 15 s, and 95 • C for 15 s. Transcript abundance was quantified using standard curves for the target and reference genes, generated from serial dilutions of PCR products from corresponding cDNAs. The transcript abundance was normalized to the reference genes MrActin and MrEF1, which showed high stability across the different muscadine genotypes and tissues used. The geometric mean of the selected housekeeping genes was validated as an accurate normalization factor.

Global Changes in Muscadine Transcriptome Through Berry Development
Our previous work built an extensive classification of metabolome profile throughout the six berry developmental stages of contrasting C5 and C6 muscadine genotypes, showing an apparent reduction in the TPC and TFC along with the progression in berry development. However, the differences between the two genotypes were manifested during the developmental stages of V, PoV, and ripening (R). The C5 genotype displayed higher TPC and TFC levels, hence more pronounced DPPH and FRAP antioxidant activities were detected (Darwish et al., 2021). To identify the molecular events associated with metabolome accumulation during berry development, the RNA-Seq strategy was used to generate transcriptome profiles with target berry developmental stages (FS, V, R, and S) from C5 and C6 genotypes cultivated under the same conditions ( Figure 1A). Even though they are complex reproductive organs, the ripe berries' seeds were considered as a fruit tissue/stage. After removing the low-quality reads, 18.8-30.5 million high-quality clean reads per replicate were acquired, which have been deposited in NCBI GenBank (PRJNA775666). The reads from the 24 libraries were mapped against the muscadine transcriptome with a 77.1-80.5% mapping rate (Supplementary Table 1). The principal component analysis (PCA) showed high consistency among transcript profiles and a clear separation of two main components. The first component (PC1) was accountable for 54% of the variance and associated with the berry development procedure. However, the second component (PC2) was responsible for 33% of the variance and mainly associated with genotype disparity, most notably at the V stage ( Figure 1B). Progressive changes to the transcript abundance were more evident when the data were subjected to hierarchical clustering (Supplementary Figure 1).
To identify the DEGs during berry development, the data were subjected to two different statistical packages, DESeq2 and EdgeR (Love et al., 2014;Chen et al., 2016). The DEGs (P FDR < 0.05, log2fold change > 1.5 or < -1.5) of each comparison generated by DESeq2 or EdgeR pipelines were considered (Supplementary Figures 2, 3 and Supplementary Tables 2, 3). Three pairwise transcriptome comparisons between consecutive berry developmental stages (V-FS, R-V, or S-R) within C5 (C5 stages ) and C6 (C6 stages ) genotypes resulted in 10,116 and 10,543 non-redundant DEGs, respectively (Figures 1C-F). In general, a number of 8518 genes were commonly upregulated or downregulated within the stages of both genotypes (Supplementary Figure 4). The data manifested that berry development involves transcriptional reprogramming of a large number of genes. In C5, there was a slight variation in the total and exclusive number of DEGs signified in each comparison ( Figure 1C). However, the number of upregulated and downregulated transcripts considerably varied with the progression in berry development. This was most apparent at the V stage, exhibiting a ∼4-fold higher number of upregulated transcripts than in the R stage (C5 R−V , Figure 1E and Supplementary Table 2). In C6, the total and exclusive numbers of DEGs in the C6 R−V comparison were ∼2-and ∼4-fold lower relative to other comparisons, respectively ( Figure 1D). However, the number of upregulated transcripts during the V stage was ∼3-fold higher than the R stage (C6 R−V ), but roughly 50% of the transcripts during the FS stage (C6 V−FS , Figure 1F and Supplementary Table 3). The data suggested that the downregulation pattern prevailed the upregulation profile of genes expression throughout berry development, particularly within C6 stages, confirming the previously reported results in different Vitis vinifera cultivars (Fasoli et al., 2012;Massonnet et al., 2017). In addition, the V stage harbored the most pronounced transcriptional modulation, especially within C5 stages.

Véraison Signifies the Main Divergent Between Genotypes During Development
To extend our knowledge of how the transcript abundance pattern differentiates in C5 and C6 during berry development, the data were subjected to the same pipelines described above. The different berry developmental stages in C5 were compared against their analogous in C6 (C5 stage -C6 stage ). The four pairwise transcriptome comparisons (C5 FS -C6 FS , C5 V -C6 V , C5 R -C6 R , and C5 S -C6 S ) identified 746, 5647, 2479, and 1996 DEGs, respectively; and 7772 non-redundant transcripts in all comparisons (Figures 2A,B, Supplementary Figure 5, and Supplementary Table 4). In agreement with PCA ( Figure 1B), the number of DEGs in the V stage comparison (C5 V -C6 V ) represented approximately 50% of the total redundant/nonredundant DEGs (Figures 2A,B). Moreover, the number of upregulated transcripts during the V stage in C5 was ∼2fold higher than in C6, while other C5 stage -C6 stage comparisons exhibited a slight variation. Surprisingly, the data showed that 41% of DEGs (5,372 out of 13,000) in C5 stage and C6 stage comparisons were shared within the C5 stages /C6 stages comparisons (Supplementary Figure 6). Only 6.6% of DEGs were exclusively present in C5 stage -C6 stage or within C5 stages analysis comparing to 9.4% of DEGs within C6 stages . The C5 stage -C6 stage data agreed with C5 stages and C6 stages comparison results, confirming that the dramatic transcriptomic changes that differentiate C5 from C6 occurred during the V stage.

Construction of Gene Co-expression Networks
The WGCNA system biology approach was utilized to construct co-expression networks based on pairwise correlations of nonlowly expressed genes of the 24 samples. Sixteen modules were identified and labeled in distinct colors shown in a hierarchical clustering dendrogram and network heatmap ( Supplementary  Figures 7, 8). Next, we performed the analysis of moduletrait correlations between the 16 modules and the biochemical property data (TPC, TFC, and TAC) or antioxidant activity data (DPPH, FRAP, ABTS, NORS, and CUPRAC) from both genotypes over berry developmental stages (ME1-ME16, FIGURE 1 | Overall experimental procedure of muscadine berries. (A) The layout of the study used the two muscadine genotypes, C5 (C5-9-1) and C6 (C6-10-1), exhibiting differential phenolic/flavonoid levels and antioxidant activity profiles (Darwish et al., 2021). The developmental stages highlighted in yellow represent the samples selected for RNA-seq, including fruit-set (FS), véraison (V), ripe skin/flesh (R), and ripe seeds (S). However, pre-and postvéraison stages (PrV and PoV) were only included in the quantitative PCR (qPCR) assay. (B) The principal components analysis (PCA) of the samples' replicates, where normalized counts were transformed to the variance stabilizing transformation (VST) using DESeq2. Each unique combination of a developmental stage was given a distinctive color. (C,D) By using DESeq2 and EdgeR pipeline, each stage was compared to its preceding one within the same genotype. As a result, 10,116 and 10,543 non-redundant differentially expressed genes (DEGs) were identified in C5 and C6, respectively, with fold2change > ± 1.5. (E,F) Bar plot of DEGs in C5 and C6 when each stage was compared to its former one. Figure 2C). Out of the 16 modules, eight modules showed significant correlations with the trait data. To link the WGCNA modules with the previously analyzed data (C5 stages , C6 stages , and C5 stage -C6 stage ), the acquired DEGs from each analysis were assigned to each module ( Figure 2C and Supplementary  Figure 9). Hence, our focus was restricted to the 11,595 nonredundant DEGs located in modules potentially associated with the phenolic/flavonoid accumulation and the antioxidant activity, FIGURE 2 | Differentially expressed genes among C5 stages and their corresponding stages in C6. (A,B) Venn diagrams and bar plots of DEGs resulted from comparing equivalent stages in C5 and C6 genotypes, using DESeq2 or EdgeR pipelines. As a result, 7772 non-redundant DEGs were identified with log2fold change > 1.5 or < -1.5. (C) Module-trait associations between RNA-seq data and evaluated biochemical-related traits, including TPC, TFC, total anthocyanin content (TAC), and different types of antioxidant assays (DPPH, FRAP, ABTS, NORS, and CUPRAC) from the C5 and C6 genotypes at different berry developmental stages. The color of the cell indicates the correlation coefficient between a given module and biochemical-related traits at the row-column intersection. Each row corresponds to a module (M1-M16). The left panel shows the assigned number of genes to each module, whether the number of total input genes or DEGs from C5 stages , C6 stages , and C5 stage -C6 stages comparisons. In addition, modules of interest were shown in bold line and selected for further analysis (indicated as in focus). Green and red are the color key that represents r 2 values from -1 to +1. FS, fruit-set; V, véraison; R, ripe skin/flesh; S, ripe seed.
To provide a broad overview of the types of DEGs located in the modules of interest, GO terms and KEGG enrichment analyses were performed based on V. vinifera Ensembl GeneID (Raudvere et al., 2019; Supplementary Table 6). Among the eight modules of interest, four modules (ME1, ME5, ME11, and ME14) exhibited high enrichment in GO terms in the molecular function (MF) for transferase and oxidoreductase activities. Moreover, ME1 and ME11 were highly enriched in GO terms belonging to the biological process (BP) for flavonoid metabolic process (GO:0009812) and phenylpropanoid metabolic process (GO:0009698), respectively (Supplementary  Figures 10, 12). Similarly, ME1, ME5, ME11, and ME14 were highly enriched in the KEGG pathway for isoflavonoid (vvi00943), diterpenoid (vvi00904), phenylpropanoid (vvi00940), and terpenoid (vvi00900) backbone biosynthesis, respectively (Supplementary Figures 10-13 and Supplementary Table 6). Moreover, ME1 and ME5 exhibited significant enrichment for flavonoid biosynthesis (vvi00941). To identify the major transcriptional dynamics associated with berry development, K-means clustering was performed to partition DEGs within modules of interest into specific behaviors. The 11,595 DEGs were classified into 15 K-means clusters (K1-K15) that were comparable but with different kinetics in both genotypes (Supplementary Figures 14, 15 and Supplementary  Tables 7, 8). Interestingly, C5 showed a complete K7 cluster that has no comparable pattern in C6. Temporal expression trends within each module were characterized by building modules/K-clusters Circos. Although the modules were represented by the same number of DEGs, they displayed distinctive profiles between C5 and C6 ( Supplementary  Figures 16, 17).

Identification of Genes in the Polyphenolics Pathway
To identify the key genes involved in phenolic/flavonoid biosynthesis, annotation information of DEGs within modules of interest was extracted from the de novo muscadine reference genome annotation (Park et al., 2021). A number of 94 genes were identified based on their predicted functions (Supplementary Table 9). These hub genes were located in six modules (ME1, ME5, ME6, ME11, ME13, and ME14) and different K-means clusters (Supplementary Figure 18). The GO term and KEGG analyses among the candidate genes showed significant enrichment of biological pathways for shikimate metabolic process, glutathione metabolism, drug transmembrane transport, and drug transport, along with flavonoid, flavone, flavonol, phenylalanine, phenylpropanoid, anthocyanin, stilbenoid, and lignin biosynthetic process (Figure 3 and Supplementary Table 10).
To validate the expression patterns and identify the crucial candidate genes underlying the diversity in phenolic/flavonoid accumulation and antioxidant activity in both genotypes, we performed qPCR to the hub genes at different developmental stages. Since the variance associated with the genotype (PC2) was remarkably evident at the V stage (Figures 1A,B), we decided to include the PrV and PoV stages in the experiment. The qPCR results of selected genes at FS, V, R, and S stages were significantly correlated (r 2 ≥ 0.933, P < 9.8 × 10 −7 ) with their TPM values, confirming the RNA-seq expression data s (Supplementary Table 11). For convenience, the 94 genes were classified into six groups based on their location in the phenolic/flavonoid biosynthesis (Figures 4-8). Generally, all genes showed comparable expression patterns but with different kinetics between genotypes.
For genes involved in the shikimic and amino acids pathway (group-1), 3 out of 12 genes showed distinct accumulation patterns, including the nicotinamidase (NIC), the arogenate dehydratase/prephenate dehydratase (ADT), and the UDPphenylpropanoid glucosyltransferases (TOGT) (Figure 4 and Suppplementary Figure 19). Relative to C6, NIC transcripts steadily increased to a peak at the PrV stage. The expression level of NIC was ∼2-and ∼3-fold higher in C5 at the PrV and V stages. Similarly, the abundance of ADT was visibly higher in C5 (∼9.5fold) at the same stages. Interestingly, a strong peak signal of TOGT was detected at the PrV stage of C5, exhibiting ∼175-fold higher transcript abundance and sharply decreased afterward. In general, all genes related to the shikimic and amino acids pathway were more abundant in C5 than in the C6 genotype, particularly at the PrV and V stages, excluding the chorismate mutase2 gene (CM2).
The comparable accumulation pattern was also noted for gene-encoded enzymes involved in the hydroxycinnamic acids pathway (group-2). Only two out of the 14 genes coordinating the pathway exhibited differential expression profiles between the two genotypes ( Figure 5 and Supplementary Figure 20). The omega-hydroxypalmitate O-feruloyltransferase2 (HHT2) expression displayed a shifted pattern in C6 genotype toward induction during the late stages of PoV and R. However, the expression pattern of the caffeoyl-CoA O-methyltransferase (CCoAOMT) was indistinguishable between the two genotypes before it peaked by ∼6.5-fold at the R stage of C6.
Regarding the lignin and pterostilbene pathway (group-3), the cinnamyl alcohol dehydrogenase 1 (CAD1) and the vestitone reductase (VR) genes displayed high accumulation levels that occurred mainly at the PrV stage, reaching ∼7.5and ∼3.5-fold higher in C5 genotype, respectively (Figure 6 and Supplementary Figure 21). However, this pattern was inversed in the case of stilbene synthase genes (STS1-3) that manifested higher levels at the V, PoV, and R stages of C6, suggesting the intense activation of the pterostilbene pathway during the late developmental stages of the C6 genotype.
In the flavonoids/anthocyanins pathway (group-4), several genes exhibited similar accumulation profiles between the two genotypes, considering that both genotypes are colored, albeit higher transcription levels were detected in C5 during PrV, V, and PoV stages such as chalcone isomerase (CHI), FIGURE 3 | A network view for the predefined biological processes Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEEG) pathways for the 94 hub genes. BP GO terms and KEEG pathways for the selected 94 genes (p-adjusted < 0.05), extracted by g:Profiler website with Benjamini-Hochberg false discovery rate (FDR) multiple testing correction method. The default ClueGO settings were applied, and the terms are functionally grouped based on shared genes (kappa score). The size of the nodes indicates the number of mapped genes ranged from 0-5, 5-10, 10-20, 20-30, and ≥30 genes. The size of the nodes indicates the number of mapped genes, while the color indicates the degree of significance (0.1 < p-value < 0.0005). The most significant term defines the name of the group.
In group-5, different TF-related genes either acting as negative (MYB86, MYB17, and SPLs) or positive (TCPs, AGLs, MYB90, ANL2, bHLH2/3, and WD40) regulators of the phenylpropanoid/flavonoid pathway generally displayed a similar accumulation pattern between the two genotypes. However, higher mRNA levels were detected in C5 during PrV, V, and PoV stages albeit a few exceptions. They showed a rapid decline in their transcription levels along with berry development. Nevertheless, three TFs displayed strong mRNA signals detected at one or more developmental stages of the C5 or C6 genotype (Figure 8 and Supplementary Figure 23). The major differences between C5 and C6 genotypes were evident during the V, PoV, and R stages, particularly for the MYB-transcriptional activators (MYBF1 and MYB1L), the homeobox-leucine zipper protein anthocyaninless (ANL1), and bHLH1. The MYBF1 sharply peaked at the PoV, reaching ∼75fold higher in C5 compared with C6 that retained at basal low levels after the PrV (Figure 8; Li and Zachgo, 2013). The expression of MYB1L was initially high at the FS stage and rapidly declined in both genotypes; however, it exhibited a strong induction, reaching a peak at the PoV stage with ∼9-fold higher in C5. Similarly, the bHLH1 displayed more abundance at the V stage that continues to a peak at the PoV stage of C5. In C5, bHLH1 exhibited ∼6.7-and ∼4.8-fold higher than in C6 at corresponding stages, respectively. Interestingly, ANL1 showed more transcript abundance at the V, PoV, and R stages of C6. Overall, these results highlight genotype-specific trends in the accumulation of the TFs, coordinating the expression of downstream genes associated with the phenolic/flavonoid pathway by which they are abundant mainly during immature developmental stages of both genotypes.
Finally, the genes in group-6 exhibited a similar expression pattern. The peroxidase proteins (PODs) are involved in the turnover of vacuolar phenolic metabolites (Barceló et al., 2003). The methanol O-anthraniloyltransferase (AMAT), the flavin-containing monooxygenase (FMO), and the UDPrhamnose: rhamnosyltransferase (GRT4s) are enzymes involved in the production of O-methyl anthranilate, formation of the epicatechin 3 -O-glucoside, and metabolism of phenolic compounds, respectively (Wang and De Luca, 2005;Nakamura et al., 2012;Hsu et al., 2017). The PODs, AMAT, GRT4s, and FMO mRNAs declined rapidly in both genotypes, although pronounced levels were retained in C5 at one or more developmental stages (Supplementary Figure 24). By contrast, one version of the genes encoded the heavy metal transport/detoxification protein (HMA2) showed expression induction at the V stage of C5, while the HMA1 was more abundant in the seeds of both genotypes.

DISCUSSION
Despite the nutraceutical qualities of muscadine berries, no previous investigation has been conducted to elucidate the molecular basis of potent berry bioactivity. The current investigation provides the first extensive record for mRNA expression profiling in parallel with the changes in TPC/TFC accumulation and antioxidant activity throughout muscadine berry development. The two contrasting muscadine genotypes, C5/C6, were selected based on their diversity in TPC/TFC levels that influence antioxidant capacity (Darwish et al., 2021). The differences between these genotypes exhibiting similar berry color and grown under the same cultivation circumstances were sufficient to identify the critical DEGs coordinating antioxidant bioactivity throughout berry development in a genotypedependent manner. The data indicated that transcriptional modulation of a large number of genes occurs during berry development. The PCA of the entire transcript profiles showed a clear separation among berry developmental stages. Likewise, the berry transcriptome of different white and red berry grape varieties exhibited a similar extensive shift along with development (Massonnet et al., 2017). Most importantly, the PC2 pointed out the V stage as a decisive event when the separation between genotypes occurred in terms of the levels of TPC/TFC and antioxidant activity. Indeed, the number of upregulated transcripts during C5 V was higher than in C6 V . Furthermore, comparisons among stages showed that the FIGURE 6 | Lignin and pterostilbene pathway. The colored boxes show the accumulation pattern of genes involved in lignin and pterostilbene pathway that were assessed by qPCR at different berry developmental stages, including fruit-set (1), prevéraison (2), véraison (3), postvéraison (4), ripe-skin/flesh (5), and ripe-seeds (6). Red boxes indicate higher levels of expression, and green boxes indicate lower expression levels. The color brightness is directly proportional to the expression ratio, according to the color scale. UCC, uclacyanin-3; MVC, mavicyanin; ENODL, early nodulin; CAD, cinnamyl alcohol dehydrogenase; CCR, cinnamoyl-CoA reductase; VR, vestitone reductase; LAC, laccase; STS, stilbene synthase; ROMT, trans-resveratrol di-O-methyltransferase. Underlined genes, including PAL, C4H, 4CL, and CCoAOMT are shown in Figure 5. number of upregulated transcripts during the V stage was higher than through the R stage. The data demonstrated the predominance of downregulation mode of expression during berry development, but not for the seed (S) tissue, which is consistent with previous results using different grape cultivars (Fasoli et al., 2012;Massonnet et al., 2017). This downregulation profile was more pronounced within C6 stages . Moreover, the most transcriptional modulation between C5 and C6 arises during the V stage, affecting mainly phenolic/flavonoid accumulation. Therefore, the piling up of phenolic/flavonoid and antioxidant activity appeared to be directly associated with the transcriptomic modulation of berry development, initiated at the V stage of C5. By contrast, the seeds did not display such considerable differences between genotypes.
In general, berry development is accompanied by alterations in secondary metabolites accumulation, considering it as one of the main factors influencing berry quality. However, the relationship between the phenolic/flavonoid content and muscadine grapes has not been explored yet at a molecular level. The transcriptome profile of the muscadine genotypes during berry development emphasized the remarkable transcriptional shift at the V stage. The identified DEGs output of C5 stages , C6 stages , and C5 stage -C6 stage comparisons are highly informative about phenolic/flavonoid related genes that become strongly expressed and co-regulated during berry development. Therefore, in parallel with K-means clustering, we engaged the WGCNA approach to provide significant insight into the spatiotemporal dynamics underpinning differences in gene expression associated with contrasting phenolic/flavonoid content and antioxidant activity (Langfelder and Horvath, 2008). We identified eight out of 16 modules that were highly correlated with the examined traits. Interestingly, not all modules exhibited a positive association. The negatively correlated modules (i.e., ME1, ME5, and ME10) were of great interest since the reduction tendency of gene expression and phenolics accumulation prevailed during berry development (Fasoli et al., 2012;Massonnet et al., 2017;Darwish et al., 2021). To demonstrate the phenotype-genotype variation between C5 and C6, our focus was restricted to the 11,595 non-redundant DEGs located in these modules of interest. The GO and KEGG analyses confirmed the successful application of the WGCNA in analyzing the multivariate data as reported in different biological contexts (Li et al., 2019;García-Gómez et al., 2020;Wang et al., 2020). Only four modules (ME1, ME5, ME11, and ME14) were highly enriched in GO term in the BP and KEGG pathways for phenylpropanoid, flavonoid, isoflavonoid, and terpenoid backbone biosynthesis. Furthermore, we identified 94 phenolic/flavonoid biosynthetic and signaling genes located in ME1, ME5, ME6, ME11, ME13, and ME14, showing different K-means clusters. As expected, these hub genes were strongly enriched in flavone, flavonol, phenylalanine, anthocyanin, stilbenoid, lignin biosynthesis, and shikimate metabolic process, in addition to the previously stated overrepresented GO and KEGG pathways within the modules of interest. The GO and KEGG analyses supported the narrow-down strategy, offering biological knowledge about the hub genes (Consortium, 2004). Some of these overrepresented GO and KEGG pathways were reported during fruit development (Berdeja et al., 2015;García-Gómez et al., 2020). Despite the decline profile of gene expression, the C5 genotype held higher transcription levels of most of these genes mainly during the V stage, justifying the higher TPC/TFC content and antioxidant capacity.
The transcriptome survey identified both common and unique molecular events that characterize the action of muscadine berry development. One of the most distinctive molecular incidents between C5 and C6 was the massive activation of phenolic/flavonoid-related transcripts. Generally, both genotypes exhibited comparable expression patterns but with different kinetics, excluding seed tissue. Of paramount importance, the genes with a contrasting accumulation profile between genotypes expose the contribution of the genotypic diversity factor. For instance, the shikimic and amino acids pathway in C5 was emphasized by the expression of NIC, ADT, and TOGT mRNAs. The NIC encodes a nicotinamidase enzyme that converts nicotinamide to nicotinic acid, representing the most similar homolog of the 2,3-dihydro-2,3-dihydroxybenzoate FIGURE 8 | Transcription factors regulating the phenylpropanoid/flavonoid pathway. The colored boxes show the accumulation pattern of genes involved in the regulation of flavonoid and lignin pathway that were assessed by qPCR at different berry developmental stages, including fruit-set (1), prevéraison (2), véraison (3), postvéraison (4), ripe-skin/flesh (5), and ripe-seeds (6). Red boxes indicate higher levels of expression, and green boxes indicate lower expression levels. The color brightness is directly proportional to the expression ratio, according to the color scale. MYB, myb-related regulatory gene; TCP, teosinte branched1/ Cincinnata/proliferating cell factor; SPL, squamosa promoter-binding-like protein; AGL, agamous MADS-box protein; ANL, homeobox-leucine zipper protein anthocyaninless; bHLH, basic helix-loop-helix; WD40, tryptophan-aspartic acid repeat. EBGs and LBGs for early and late biosynthetic genes, respectively. synthase (EntB) (Hunt et al., 2007). The EntB is involved in a sequential step, leading to the conversion of the isochorismate to the 2,3-dihydroxybenzoic acid (2,3-DHBA) and shikimic acid in plants (Widhalm and Dudareva, 2015). TOGT encodes a phenylpropanoid glucosyltransferase functions mainly on hydroxycoumarin glycosylation (Fraissinet-Tachet et al., 1998). The TOGT-suppressed plants displayed a significant reduction in scopoletin glycosylation, the essential process for accumulating phenylpropanoids (Chong et al., 2002). Furthermore, the ADT and TYRAATs catalyze the final steps in phenylalanine and tyrosine biosynthesis, respectively (Rippert et al., 2009). The conversion of phenylalanine to cinnamic acid is the admission point of the phenylpropanoid pathway, where the latter serves as the substrate for most hydroxycinnamic acid derivatives (Alam et al., 2016). Certainly, the expression level of general phenylpropanoid-related biosynthetic genes, including PAL and 4CL was upregulated in C5. Similarly, the CAD1 and VR that are required for monolignol and isoflavonoids biosynthesis were more abundant in C5. By contrast, the three STS genes were steadily abundant in C6 (Dixon, 1999;Boudet et al., 2003). The higher transcription level of NIC, ADT, TOGT, CAD1, and VR during the PrV stage, followed by the abundance of PAL and 4CL during V and PoV stages in C5 may reflect their triggering rules that facilitate the stimulated accumulation of the phenylpropanoid derivatives.
Consistently, the C5 genotype retained higher amounts of the TPC/TFC and antioxidant capacity than in C6, albeit the decline behavior along with development (Darwish et al., 2021). The expression profiles of flavonoid/anthocyanin biosynthetic genes simulated these results. Of special interest, those genes showed a genotype-specific pattern, such as CHS, HID3, F3H2, F3 5 H, LAR2, and FOMT. The CHS catalyzes the first step in the flavonoid biosynthesis, while HID enzyme activity is a critical limiting factor in the isoflavonoids biosynthesis (Winkel-Shirley, 2001;Akashi et al., 2005;Du et al., 2010). The F3 H and F3 5 H enzymes are necessary for the biosynthesis of flavonols (quercetin and myricetin), flavan-3-ols (procyanidin and prodelphinidin), and anthocyanins-dependent compounds (cyanidin and delphinidin) (Jeong et al., 2006;Wang et al., 2014). Proanthocyanidins (PAs; also referred to procyanidins or condensed tannins) are made of oligomers and polymers of (epi)catechin units, exhibiting strong antioxidant activity (Prasain and Barnes, 2014). The LAR and ANR are major enzymes that catalyze the formation of PAs (Bogs et al., 2005;Gagné et al., 2009). Moreover, C5 exhibited upregulation of the OMT gene (FMOT) involved in the flavonoids methylation. The alteration in the enzymes underlying the modification reactions of the flavonoid pathway, such as methylation and glycosylation reactions can modify and enhance the bioactivity of the flavonoids and their derivatives (Kim et al., 2010). For instance, the O-methylation of hydroxyl groups in flavonoids reduces their reactivity and increases their antimicrobial capacity (Ibrahim et al., 1998). Finally, the flavonoids glycosylation process via UFGT results in color-intensive pigments (Boss et al., 1996). In both genotypes, the abundance of UFGT was increased toward the R stage, although higher levels were depicted in C5.
The flavonoid biosynthetic genes can be divided into early biosynthetic genes (EBGs), which catalyze the production of dihydroflavonols, and late biosynthetic genes (LBGs) that lead to the biosynthesis of proanthocyanidins and anthocyanins (Ferreyra et al., 2012;Xu et al., 2015). Flavonoid/anthocyanin biosynthesis is coordinated by a transcription complex composed of two transcription factors that belong to the R2R3-MYB and the bHLH-MYC families, and a WD40 cofactor. The three proteins cofunction by establishing the MYB-bHLH-WD40 (MBW) complex to activate the expression of a downstream cascade of flavonoid/anthocyanin structural genes (Yan et al., 2021). Despite the critical contribution of bHLH and WD40 in the complex, the MYB protein is the key component in providing specificity and determining the rate of the subsets of biosynthetic genes (Zimmermann et al., 2004). Evaluation of several mutants provided strong evidence regarding the autonomous role of MYB in mediating the transcription of EBGs that lead to the production of the colorless dihydroflavonols. However, the activation of LBGs that leads to the production of anthocyanins (color pigmentation), requires the MBW complex (Baudry et al., 2004;Stracke et al., 2007;Gonzalez et al., 2008). The MBW complex is regulated by several factors (Li, 2014). The MYB genes that act as negative regulators of the pathway suppress the formation of the MBW complex by competing with R2R3-MYBs for interactions with the bHLH component of the MBW complex (Dubos et al., 2008;Matsui et al., 2008). The squamosa promoter-binding-like protein (SPL) disrupts the MBW complex by competitively binding to the R2R3-MYB subunit of the MBW complex (Gou et al., 2011). By contrast, the positive fine-tuners, including TCP3, organize the formation of the MBW complex by synergistically associating with R2R3-MYBs proteins (Li and Zachgo, 2013). Furthermore, the lignin inhibitor agamous-like protein (AGL) stimulates the flavonoid pathway by maintaining the balance of flavonoid/lignin pathways, redirecting the metabolic flux between the two pathways (Giménez et al., 2015). Finally, the anthocyaninless, a homeodomain protein, controls the accumulation of anthocyanins in subepidermal tissue (Kubo et al., 1999).
The transcription activators MYBF1, MYB90, and MYB1L were strongly accumulated in C5. The MYBF1 acts downstream TCP to stimulate the EBGs pathway (Li and Zachgo, 2013). Consistently, expression of the EBGs, including CHS, CHI, F3H, F3 5 H, and FLS were more abundant in C5. The MYBF1 is a key flavonol regulator activating the expression of flavonol synthase, and hence accumulates flavonol at early stages of development (Czemmel et al., 2009). However, the earlier and higher levels of FS in C6 suggested undetected version(s) of FS in muscadine transcriptome. In V. vinifera, the VvMYBPA1, VvMYBPA2, and VvMYB5a are required for triggering the EBGs (Czemmel et al., 2012). The VvMYBPA1 is considered a key regulator of PAs synthesis, controlling the general flavonoid pathway and the PAs branch genes ANR and LAR Gagné et al., 2009). Interestingly, none of these TFs was captured in the transcriptome. However, the transcription activators of LBGs, MYB90, and MYB1L showed strong induction in C5. The MYB90 can activate the promoters of UFGT, anthocyanidin 3-O-glucoside 6 -O-acyltransferase (3AT), and F3 5 H (Matus et al., 2017). In both genotypes, the abundance of MYB90 was increased concomitantly with the abundance of UFGT1, although higher levels were detected in C5. It is documented that the MYB90/MYB1L acts downstream the AGL/ANL activators and the SPL repressor (Kobayashi et al., 2002;Walker et al., 2007;Xu et al., 2015). Concordantly, the expression levels of SPL sharply decreased during development, showing slight differences between genotypes. Interestingly, AGL2 and WD40 displayed higher levels during C5 berry development with the strongest signal detected in the seed of both genotypes. However, AGL1 and ANL2 contributed more during immature berry stages. One of the most interesting findings was the abundance of TFs encoded flavonoids/anthocyanin repressors, including the MYB and the SLP (Qian et al., 2017;Song et al., 2020). It could be speculated that the plenty of suppressor TFs in highly flavonoid accumulation genotypes will have a strong effect in adjusting the expression of downstream structural genes involved in the flavonoid pathway. Thus, such TF proteins can be present in plants and function as a backup system, avoiding the unnecessary excessive accumulation of flavonoids. Taken together, these results suggest that the modification reactions of flavonoids such as hydroxylation (via F3 H and F3 5 H), methylation (via FOMT), and glycosylation (via UFGT) are highly activated in C5, resulting in the production of more biologically active flavonoid derivatives.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA775666.