Lignin and Quercetin Synthesis Underlies Berry Russeting in ‘Sunshine Muscat’ Grape

In order to further explore the mechanism of ‘sunshine muscat’ grape russet formation, transcriptomic and metabolomic analyses were performed on ‘sunshine muscat’ grape peels with and without russet. A total of 1491 differentially expressed genes (DEGs) were discovered based on these analyses. The phenylpropane synthesis pathway was the key metabolic pathway identified, and 28 DEGs related to phenylpropane synthesis pathway were screened, of which 16 were related to lignin synthesis. In addition, 60 differential metabolites were screened. There were 29 phenolic substances among the differential metabolites, which were all up-regulated and 10 were quercetin-related glycosides. Our results indicate that phenols likely play a dominant role in the formation of ‘sunshine muscat’ grape russet, and the synthesis of lignin and quercetin may be the key factors underlying russet formation.


Introduction
Russet refers to the yellow-brown spots appearing on the surfaces of some fruits during their development. It is a physiological abnormality that occurs as fruits mature and is essentially a layer of secondary protective tissue that serves the function of resistance to adverse environmental conditions and protection of the fruit itself [1]. The problem of russet is prominent in the cultivation of fruit trees, in which it seriously affects the smoothness of fruit surfaces and impairs the quality and the commodity value of fruits. Therefore, this significantly decreases their commercial value [2,3]. In recent years, there have been a number of studies that have focused on apple and pear russet. However, little research has been focused on the grape russet. Researchers have found that cinnamyl-CoA reductase (CCR), cinnamyl alcohol dehydrogenase (CAD), and peroxidase (POD) genes involved in lignin biosynthesis affect russet formation in sand pears and apples [4]. Thus, lignin accumulates to promote the formation of russet by enhancing the expression of genes related to lignin synthesis [5,6]. The formation of apple russet is caused by the rupture of the stratum corneum and the formation of a cork layer. The MdSHN3 transcription factor is a positive regulator of the formation of cuticles in apple fruits and inhibits the formation of the apple russet [7]. The formation of 'Dangshansu pear' russet may be related to phenylpropane metabolism, ethylene metabolism, and secondary metabolism, while the formation of the peel russet may be regulated by lignin synthesis, polyamines, and H 2 O 2 signals [8]. Studies have shown that the formation of polyphenolic compounds is the main cause of russet formation in 'sunshine muscat' grape berries during ripening [9,10]. The synthetic pathway of phenolic compounds is quite complex. First, phosphoenolpyruvate(PEP) and erythrose phosphate (E4P) form phenylalanine through the shikimic acid pathway. Then, under the action of various enzymes, phenylalanine passes through the phenylpropanoid metabolic pathway and flavonoid pathway to form phenolic compounds [11]. Thus far, studies on 'sunshine muscat' grape russet formation are very limited, and, until now, the molecular and metabolic pathways engaged in grape russet remain elusive.
With the development of multiple omics technologies, a single omics technique can no longer meet the demands of scientific research, especially research related to plants that have more secondary metabolites. It is difficult to clarify the mechanisms through which environmental changes and other factors impact fruit quality using a single set of omics technologies. In recent years, an increasing number of researchers have applied multiple omics techniques to address research objectives, and the integrative application of multiple omics technologies is becoming more extensive. Genes and metabolites involved in the same biological process will likely exhibit similar changes. Thus, association analysis of transcriptomes and metabolomes is an effective way to find key metabolic pathways and key genes and reveal the underlying molecular mechanisms. A recent study based on metabolomics and transcriptomics revealed the mechanism of peel coloring during jujube maturation, which found that underlying flavonoid metabolites have changed and that changes in structural genes or their regulators (i.e., transcription factors) involved in flavonoid biosynthetic pathways may be the mechanism delaying the accumulation of red anthocyanins in the peels of jujube fruits [12]. This integrated metabolomics and transcriptomics approach has revealed significantly regulated metabolites and biological pathways in citrus fruits, providing new insights into the mechanism of fruit quality deterioration and the induction of resistance against Penicillium digitatum in Citrus fruit [13]. Such an analysis has been conducted on the genes and metabolites encoding stilbene synthase (a key enzyme involved in resveratrol synthesis) in grape peel after UV-C irradiation using full transcriptomics sequencing and metabolomics techniques. The anabolism network of astragalus in fruit after UV-A irradiation has been conducted for the first time [14].
In order to explore the mechanism of grape russet formation, 'sunshine muscat' grapes were examined in this study. Russet and non-russet grapes peel samples were collected at five different stages and mixed uniformly. Transcriptomic sequencing, metabolomic detection, and association of transcriptomic and metabolomics analyses were conducted to screen for the key metabolic pathways and identify differentially expressed genes (DEGs) and metabolites related to russet formation. Knowledge about the physiological and molecular mechanisms underlying 'sunshine muscat' grape russet formation will establish a theoretical foundation for the screening of 'sunshine muscat' grape russet control measures and high-quality fruit formation.

Plant Materials
'Sunshine muscat' grapes were harvested from a vineyard located at the modern agriculture research and development base of Sichuan Agricultural University (30 • 33 46" N, 103 • 39 36" E) in July 2019 (the critical period before russet formation). The row spacing of the vineyard was 1.5 m × 3.0 m, and the four-year-old vines were characterized by the uniformity of their vigor and soil and water management in the field. Sampling was performed 70, 80, 90, 100, and 110 days after flowering, and 60-80 berries were collected each time. To study the russet formation mechanism, russet and non-russet grape peels were separated from collected grapes (Figure 1), frozen in liquid nitrogen, and stored at −80 • C in an ultra-low temperature refrigerator until subsequent use. The russet and non-russet grapes peels, which were collected across five stages, were separately and uniformly mixed, and the mixed russet peels (CKY) and mixed non-russet peels (CKN) were used for transcriptomic sequencing and metabolomic assays. Biomolecules 2020, 10, x FOR PEER REVIEW 3 of 19 Figure 1. 'Sunshine muscat' grapes with different degrees of the russet phenotype: (A) non-russet grape, (B) slightly russet grape, and (C) extremely russet grape.

Total RNA Extraction and Sequencing
Grape peels with and without russet were separately mixed at equal quantities from each stage, and RNA was extracted. Each sample consisted of three replicates for RNA-seq sequencing. Grape peels were ground with liquid nitrogen for RNA extraction. The RNAprep Pure Polysaccharide polyphenol plant total RNA extraction kit (Qiagen, Hilden, Germany) was used for RNA extraction, and the constructed sequencing library (prepared with the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® was used to construct the library) was sequenced on the Illumina Hi SeqTM4000 platform (Illumina, San Diego, CA, USA).

Data Quality Control
All raw sequencing reads were first processed with in-house Perl scripts (Novogene Technologies Co., Ltd., Beijing, China). In order to ensure the quality and reliability of data analysis, the original data were filtered through the removal of reads containing adapter sequence, uncalled 'N' bases, and too many low-quality bases (i.e., reads for which Qphred ≤ 20 bases accounted for more than 50% of the total read length). At the same time, Q20 (the percentage of bases with phred > 20, phred = −10log10(e)), Q30 (the percentage of bases with phred > 30, phred = −10log10(e)), and GC (the percentage of G and C among bases called in clean reads) content calculations were performed on the clean reads. All subsequent analyses were high-quality analyses based on clean data.

Differentially Expressed Gene Analysis and Enrichment Analysis
RSEM software [15] was used to quantitatively analyze the gene expression levels of each sample. Unigene expression was expressed as fragments per kilobase of exon per million fragments mapped (FPKM) values. The method of Benjamini and Hochberg was used to correct the p-values of differentially expressed genes (DEGs) for multiple corrections. The thresholds of differentially expressed unigenes were an absolute fold change >1 and p < 0.05 [16], which were used as standards to identify genes that were significantly differentially expressed (where fold change was the FPKM ratio of each unigene between the russet and non-russet grape peel libraries). Cluster Profiler (3.4.4) software was applied to conduct Gene Ontology (GO) enrichment analysis of DEGs and statistical enrichment of DEGs in various Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

qRT-PCR Analysis
Sixteen genes involved in lignin biosynthesis were selected to validate the results of RNA-seq by qRT-PCR. The primers used in qRT-PCR are shown in Table 1. The TB Green Premix Ex Taq II (Tli RNaseH Plus, Takara, Beijing, China) was used to perform qRT-PCR. The thermal cycling program for qPCR was 95 °C for 30 s, which was followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. VvGAPDH was used as the reference gene for normalization. The relative expression levels of the genes were analyzed by the comparative threshold cycle(CT)method (2 − ΔΔ CT method) [17].

Total RNA Extraction and Sequencing
Grape peels with and without russet were separately mixed at equal quantities from each stage, and RNA was extracted. Each sample consisted of three replicates for RNA-seq sequencing. Grape peels were ground with liquid nitrogen for RNA extraction. The RNAprep Pure Polysaccharide polyphenol plant total RNA extraction kit (Qiagen, Hilden, Germany) was used for RNA extraction, and the constructed sequencing library (prepared with the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® was used to construct the library) was sequenced on the Illumina Hi SeqTM4000 platform (Illumina, San Diego, CA, USA).

Data Quality Control
All raw sequencing reads were first processed with in-house Perl scripts (Novogene Technologies Co., Ltd., Beijing, China). In order to ensure the quality and reliability of data analysis, the original data were filtered through the removal of reads containing adapter sequence, uncalled 'N' bases, and too many low-quality bases (i.e., reads for which Q phred ≤ 20 bases accounted for more than 50% of the total read length). At the same time, Q20 (the percentage of bases with phred > 20, phred = −10log 10 (e)), Q30 (the percentage of bases with phred > 30, phred = −10log 10 (e)), and GC (the percentage of G and C among bases called in clean reads) content calculations were performed on the clean reads. All subsequent analyses were high-quality analyses based on clean data.

Differentially Expressed Gene Analysis and Enrichment Analysis
RSEM software [15] was used to quantitatively analyze the gene expression levels of each sample. Unigene expression was expressed as fragments per kilobase of exon per million fragments mapped (FPKM) values. The method of Benjamini and Hochberg was used to correct the p-values of differentially expressed genes (DEGs) for multiple corrections. The thresholds of differentially expressed unigenes were an absolute fold change >1 and p < 0.05 [16], which were used as standards to identify genes that were significantly differentially expressed (where fold change was the FPKM ratio of each unigene between the russet and non-russet grape peel libraries). Cluster Profiler (3.4.4) software was applied to conduct Gene Ontology (GO) enrichment analysis of DEGs and statistical enrichment of DEGs in various Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

qRT-PCR Analysis
Sixteen genes involved in lignin biosynthesis were selected to validate the results of RNA-seq by qRT-PCR. The primers used in qRT-PCR are shown in Table 1. The TB Green Premix Ex Taq II (Tli RNaseH Plus, Takara, Beijing, China) was used to perform qRT-PCR. The thermal cycling program for qPCR was 95 • C for 30 s, which was followed by 40 cycles of 95 • C for 5 s and 60 • C for 30 s. VvGAPDH was used as the reference gene for normalization. The relative expression levels of the genes were analyzed by the comparative threshold cycle(CT)method (2 −∆∆CT method) [17].

Measurement and Analysis of Metabolomes
Liquid chromatography-mass spectrometry (LC-MS) technology [18,19] based on the highly sensitive SCIEX QTRAP ® 6500+ mass spectrometry platform (SCIEX, Framingham, MA, USA) was used to conduct a quasi-targeted metabolomics assay. We collected 100-mg tissue samples, which were frozen in liquid nitrogen in an eppendorf(EP) tube, to which 500 µL of 80% methanol aqueous solution containing 0.1% formic acid was added. Samples were then whirled, shocked, immersed in an ice bath for 5 min, and centrifuged at 4 • C and 15,000 rpm for 10 min. A volume of supernatant was added to a one-half volume of mass spectrometry water diluted to 53% methanol content. Supernatant samples were placed in a centrifuge tube and centrifuged at 4 • C and 15,000× g for 20 min. Then, the supernatant was collected for LC-MS analysis. An equal volume of the sample was taken from each experimental sample and mixed for use as a quality control (QC) sample. Based on the Novogene database, the multi-response monitoring model (MRM) was used to assess experimental samples, and the data were analyzed by principal component analysis (PCA), partial least squares-discriminant analysis (PLS-DA), and other multivariate statistical analyses, with the aim of elucidating differences in metabolites between the russet and non-russet grapes.

Transcriptome Sequencing Yield Statistics
The constructed grape peel transcriptome libraries were sequenced using the high-throughput Illumina platform, obtaining 87,831,088 and 84,645,916 raw reads for the russet CKY and non-russet CKN libraries, respectively. After removing low-quality sequences, 85,425,282 and 82,194,414 clean reads were obtained for the russet CKY and non-russet CKN libraries, respectively. The clean bases totaled 6.45 G and 6.17 G after filtering the CKY and CKN libraries, respectively. Q20 and Q30 sequences comprised more than 95% (sequencing error rate <1%) and 90% of all clean reads, respectively, and the GC content exceeded 45% (Table 2), which indicates that the sequencing data volume and quality were both satisfactory for subsequent sequence assembly and analysis.

Analysis of Differentially Expressed Genes and Gene Enrichment in CKY versus CKN Libraries
After screening DEGs according to adjusted p-values (P adjusted < 0.05), a total of 1491 significantly DEGs were identified between the CKY and CKN transcriptomes, of which 574 and 917 genes were down-regulated and up-regulated, respectively, in CKY ( Figure 2). was taken from each experimental sample and mixed for use as a quality control (QC) sample. Based on the Novogene database, the multi-response monitoring model (MRM) was used to assess experimental samples, and the data were analyzed by principal component analysis (PCA), partial least squares-discriminant analysis (PLS-DA), and other multivariate statistical analyses, with the aim of elucidating differences in metabolites between the russet and non-russet grapes.

Transcriptome Sequencing Yield Statistics
The constructed grape peel transcriptome libraries were sequenced using the high-throughput Illumina platform, obtaining 87,831,088 and 84,645,916 raw reads for the russet CKY and non-russet CKN libraries, respectively. After removing low-quality sequences, 85,425,282 and 82,194,414 clean reads were obtained for the russet CKY and non-russet CKN libraries, respectively. The clean bases totaled 6.45 G and 6.17 G after filtering the CKY and CKN libraries, respectively. Q20 and Q30 sequences comprised more than 95% (sequencing error rate <1%) and 90% of all clean reads, respectively, and the GC content exceeded 45% (Table 2), which indicates that the sequencing data volume and quality were both satisfactory for subsequent sequence assembly and analysis.

Analysis of Differentially Expressed Genes and Gene Enrichment in CKY Versus CKN Libraries
After screening DEGs according to adjusted p-values (Padjusted < 0.05), a total of 1491 significantly DEGs were identified between the CKY and CKN transcriptomes, of which 574 and 917 genes were down-regulated and up-regulated, respectively, in CKY ( Figure 2). Based on the Gene Ontology (GO) database, GO functional clustering was performed on the DEGs obtained by sequencing (Padjusted < 0.05). GO classifications are divided into three categories: cell components (CC), biological processes (BP), and molecular function (MF). The identified DEGs provided insights into the molecular mechanisms related to russet formation and were subjected to GO analysis to identify their enrichment in various roles. A total of 35 enriched GO terms were found, including 2 cell components, 13 biological processes, and 35 molecular functions. The CC Based on the Gene Ontology (GO) database, GO functional clustering was performed on the DEGs obtained by sequencing (P adjusted < 0.05). GO classifications are divided into three categories: cell components (CC), biological processes (BP), and molecular function (MF). The identified DEGs provided insights into the molecular mechanisms related to russet formation and were subjected to GO analysis to identify their enrichment in various roles. A total of 35 enriched GO terms were found, including 2 cell components, 13 biological processes, and 35 molecular functions. The CC classification mainly included the overall components of membranes and cell wall morphology. The BP classification mainly included biosynthetic processes and biological stress responses, and the MF classification mainly included oxidoreductase activity, hydrolase activity, antioxidant activity, peroxidase activity, and carbohydrate binding ( Figure 3). classification mainly included the overall components of membranes and cell wall morphology. The BP classification mainly included biosynthetic processes and biological stress responses, and the MF classification mainly included oxidoreductase activity, hydrolase activity, antioxidant activity, peroxidase activity, and carbohydrate binding ( Figure 3). Based on the KEGG database, the candidate metabolic pathways related to grape russet formation were identified to further study the biological functions of these genes. As seen in Table 3, the metabolic pathways related to grape russet formation were mainly phenylpropane synthesis, plant hormone signal transduction, and glutathione metabolism. The number of up-regulated DEGs was greater than the number of down-regulated DEGs generally (Padjusted < 0.05). The number of genes up-regulated in secondary metabolic pathways, such as phenylpropane synthesis, plant hormone signal transduction, and glutathione metabolism, was significantly greater than that of downregulated genes, which indicates that the formation of 'sunshine muscat' grape russet might be related to these up-regulated genes.  Based on the KEGG database, the candidate metabolic pathways related to grape russet formation were identified to further study the biological functions of these genes. As seen in Table 3, the metabolic pathways related to grape russet formation were mainly phenylpropane synthesis, plant hormone signal transduction, and glutathione metabolism. The number of up-regulated DEGs was greater than the number of down-regulated DEGs generally (P adjusted < 0.05). The number of genes up-regulated in secondary metabolic pathways, such as phenylpropane synthesis, plant hormone signal transduction, and glutathione metabolism, was significantly greater than that of down-regulated genes, which indicates that the formation of 'sunshine muscat' grape russet might be related to these up-regulated genes.

Differentially Expressed Genes Involved in Phenylpropane Synthesis
Differentially expressed unigenes were assessed at thresholds of a corrected absolute log 2 (Fold Change) > 1 and p < 0.05 in order to screen genes that were significantly differentially expressed. Analyzing and screening DEGs in the phenylpropane synthesis pathways between CKY and CKN libraries yielded 28 related DEGs of which 17 and 11 were up-regulated and down-regulated, respectively (Table 4). Based on the number of up-regulated and down-regulated genes and the involved pathways, it appears that grape peel russet is connected to the up-regulation of caffeic acid 3-O-methyltransferase (100854172), cinnamyl alcohol dehydrogenase (100262421), hydroxycinnamoyl transfer enzyme (100262421), caffeoyl coenzyme AO methyltransferase (100233087), peroxidase (100854817, 100249955, 100242338, and 100262575), 4-coumarin-CoA ligase (100254698), and cinnamyl-CoA ligase (100251623) genes, which are involved in lignin synthesis. Thus, synthesis of lignin in peels might underlie the formation of the grape russet.

Reliability Validation of RNA-seq by qRT-PCR
In order to further confirm the reliability and accuracy of Illumina RNA-seq analysis results, 16 key genes involved in lignin biosynthesis were selected, and their expression was evaluated by real-time fluorescent quantitative PCR(qRT-PCR). The expression patterns of these genes obtained from qRT-PCR were highly consistent with those in the RNA-seq data (Figure 4), which indicates that the RNA-seq data were reliable.

Metabolite PCA and PLS-DA
The multidimensional statistical analysis method principal component analysis (PCA) was conducted to classify the samples. Because no external factors were considered, the obtained PCA model reflected the overall differences in the metabolome data, which clarified the metabolites present in the skins of grapes with and without russet. The distance between the CKY and CKN samples was very great, which indicated a substantial difference between the metabolites in the CKY and CKN treatments ( Figure 5A).
Partial least squares discrimination analysis (PLS-DA) utilized partial least squares regression [20] to establish a relationship model between metabolite expression and sample type, which aimed to predict sample types based on the metabolites observed. The PLS-DA models of each comparison group were established, and the model evaluation parameters (R 2 , Q 2 ) were obtained through 7-fold cross-validation. R 2 and Q 2 values closer to 1 indicated a model was more stable and reliable. The results of the PLS-DA score map ( Figure 5B) and displacement test for CKY and CKN samples ( Figure  5C) were similar to the PCA analysis results. The model coefficients of CKY and CKN samples were Q 2 Y = 0.94 and R 2 Y = 1.00, which indicated that the model had high predictive ability and goodness-

Metabolite PCA and PLS-DA
The multidimensional statistical analysis method principal component analysis (PCA) was conducted to classify the samples. Because no external factors were considered, the obtained PCA model reflected the overall differences in the metabolome data, which clarified the metabolites present in the skins of grapes with and without russet. The distance between the CKY and CKN samples was very great, which indicated a substantial difference between the metabolites in the CKY and CKN treatments ( Figure 5A).
Partial least squares discrimination analysis (PLS-DA) utilized partial least squares regression [20] to establish a relationship model between metabolite expression and sample type, which aimed to predict sample types based on the metabolites observed. The PLS-DA models of each comparison group were established, and the model evaluation parameters (R 2 , Q 2 ) were obtained through 7-fold cross-validation. R 2 and Q 2 values closer to 1 indicated a model was more stable and reliable. The results of the PLS-DA score map ( Figure 5B) and displacement test for CKY and CKN samples ( Figure 5C) were similar to the PCA analysis results. The model coefficients of CKY and CKN samples were Q 2 Y = 0.94 and R 2 Y = 1.00, which indicated that the model had high predictive ability and goodness-of-fit. In the displacement test, the model Q 2 point was far lower than the rightmost original Q 2 point from left to right. The rightmost R 2 and Q 2 values were greater than 0.9, which indicates that the model had better predictive ability and was effective and usable.
Biomolecules 2020, 10, x FOR PEER REVIEW 10 of 19 of-fit. In the displacement test, the model Q 2 point was far lower than the rightmost original Q 2 point from left to right. The rightmost R 2 and Q 2 values were greater than 0.9, which indicates that the model had better predictive ability and was effective and usable.

Differential Metabolite Analysis
The variable importance in the projection (VIP) value of the first principal component of the PLS-DA model represents the contribution rate of metabolite differences in different groups, which was assessed [21]. Additionally, the fold change (FC) expresses the ratio of the mean of the repeated quantitative values of all organisms in the comparison group for each metabolite. FC and VIP were combined with p-values of t-tests between metabolite expression levels to assess metabolites at thresholds of VIP > 1.0, FC > 1.2 (or FC < 0.833), and p < 0.05 [22][23][24]. In total, 443 metabolites were screened, which revealed 60 differential metabolites through this analysis, including 43 differentially up-regulated metabolites and 17 differentially down-regulated metabolites (Table 5). Among the 60 different metabolites screened, there were 29 phenolic substances up-regulated as well as 8 nucleic acids and their derivatives, 3 amino acids and their derivatives, and 20 other substances (Table 6).

Differential Metabolite Analysis
The variable importance in the projection (VIP) value of the first principal component of the PLS-DA model represents the contribution rate of metabolite differences in different groups, which was assessed [21]. Additionally, the fold change (FC) expresses the ratio of the mean of the repeated quantitative values of all organisms in the comparison group for each metabolite. FC and VIP were combined with p-values of t-tests between metabolite expression levels to assess metabolites at thresholds of VIP > 1.0, FC > 1.2 (or FC < 0.833), and p < 0.05 [22][23][24]. In total, 443 metabolites were screened, which revealed 60 differential metabolites through this analysis, including 43 differentially up-regulated metabolites and 17 differentially down-regulated metabolites (Table 5). Among the 60 different metabolites screened, there were 29 phenolic substances up-regulated as well as 8 nucleic acids and their derivatives, 3 amino acids and their derivatives, and 20 other substances (Table 6).

Analysis of Relative Contents of Key Differential Metabolites
A total of 10 quercetin-related glycosides were found among 29 different phenolic metabolites (Table 6) Figure 6). Among the 10 quercetin glycosides, the quercetin glycoside content of the CKY sample was significantly higher than that of the CKN sample. Thus, we deduced that the quercetin-related compounds from the phenylpropane metabolic pathway were likely to be the key metabolites involved in the formation of 'sunshine muscat' grape russet since quercetin itself is unstable and mostly existed in the form of quercetin glycosides.

Analysis of Relative Contents of Key Differential Metabolites
A total of 10 quercetin-related glycosides were found among 29 different phenolic metabolites (Table 6) Figure 6). Among the 10 quercetin glycosides, the quercetin glycoside content of the CKY sample was significantly higher than that of the CKN sample. Thus, we deduced that the quercetin-related compounds from the phenylpropane metabolic pathway were likely to be the key metabolites involved in the formation of 'sunshine muscat' grape russet since quercetin itself is unstable and mostly existed in the form of quercetin glycosides.

Analysis of Different Metabolites and Metabolic Pathways
Hierarchical cluster analysis was performed on the differential metabolites between the CKY and CKN samples. After clustering analysis of the identified metabolites, clearly up-regulated or down-regulated metabolites were visible, and the two groups could be distinguished according to metabolite expression ( Figure 7A). Phenylpropane biosynthesis, including flavonoid and flavonol Figure 6. Relative content of key differential metabolites in russet CKY versus non-russet CKN.

Analysis of Different Metabolites and Metabolic Pathways
Hierarchical cluster analysis was performed on the differential metabolites between the CKY and CKN samples. After clustering analysis of the identified metabolites, clearly up-regulated or down-regulated metabolites were visible, and the two groups could be distinguished according to metabolite expression ( Figure 7A). Phenylpropane biosynthesis, including flavonoid and flavonol biosynthesis, and differential metabolite accumulation of products of pyrimidine metabolism were related to the 'sunshine muscat' grape russet phenotype ( Figure 7B).
Biomolecules 2020, 10, x FOR PEER REVIEW 13 of 19 biosynthesis, and differential metabolite accumulation of products of pyrimidine metabolism were related to the 'sunshine muscat' grape russet phenotype ( Figure 7B).

Association Analysis between Transcriptomic and Metabolomic Data
Association analysis was performed based on the Pearson correlation coefficient between genes that were significantly different in the transcriptomic analysis, and metabolites that were significantly different in the metabolomic analysis were assessed for the degree of correlation between the differential genes and differential metabolites. There was a significant correlation between the DEGs and differential metabolites in 'sunshine muscat' grape peels with and without russet ( Figure 8A). All the DEGs and metabolites obtained from the russet and non-russet peels of 'sunshine muscat' grapes were mapped to the KEGG pathway database to obtain their common pathways and identify the major biochemical pathways and signal transduction pathways in which differentially expressed metabolites and genes were involved. The transcriptomic and metabolomic differences were mainly enriched in phenylpropane metabolic pathways, which is followed by starch and sugar metabolism pathways, and then the glutathione glycine metabolism pathway ( Figure 8B). Sample clustering is shown in the branching diagram at the top of (A), while the metabolite clustering is shown in the branching diagram at the left side of (A). Shorter branches indicate higher similarity. Red represents up-regulation, while blue represents down-regulation (A). The abscissa in (B) is x/y (where x is the number of differential metabolites in the corresponding metabolic pathway and y is the total number of metabolites identified in the pathway). The colors of the points represent the p-values of a hyper geometric test. The smaller the value, the greater the reliability and statistical significance. The size of the point represents the number of differential metabolites in the corresponding pathway. The larger the point, the more differential metabolites occur in the pathway (B).

Association Analysis between Transcriptomic and Metabolomic Data
Association analysis was performed based on the Pearson correlation coefficient between genes that were significantly different in the transcriptomic analysis, and metabolites that were significantly different in the metabolomic analysis were assessed for the degree of correlation between the differential genes and differential metabolites. There was a significant correlation between the DEGs and differential metabolites in 'sunshine muscat' grape peels with and without russet ( Figure 8A). All the DEGs and metabolites obtained from the russet and non-russet peels of 'sunshine muscat' grapes were mapped to the KEGG pathway database to obtain their common pathways and identify the major biochemical pathways and signal transduction pathways in which differentially expressed metabolites and genes were involved. The transcriptomic and metabolomic differences were mainly enriched in phenylpropane metabolic pathways, which is followed by starch and sugar metabolism pathways, and then the glutathione glycine metabolism pathway ( Figure 8B).
Biomolecules 2020, 10, x FOR PEER REVIEW 14 of 19 Figure 8. Association analysis between differentially expressed metabolites and genes according to KEGG pathways. Correlation analysis (A). Correlation coefficients less than 0 describe negative correlations, shown in blue, while those greater than 0 describe a positive correlation, shown in red.
The branching diagram at the top of (A) represents clusters of differentially expressed genes, while the branching diagram at the left of (A) represents clusters of differentially expressed metabolites. The shorter the clustering branch, the higher the similarity. KEGG pathway analysis (B). The abscissa is the ratio of the differentially expressed metabolites or genes enriched in the pathway to the number of metabolites or genes annotated in the pathway (Ratio), and the ordinate is the KEGG pathway Lignin and quercetin are both products of the phenylpropane metabolic pathway. The expression levels of 4CL, CAD, HCT, CCR, CCo AOMT, and COMT were significantly higher in russet CKY samples compared with non-russet CKN samples. POD expression was observed to be both up-regulated and down-regulated. In addition, PAL expression was significantly lower in CKY samples when compared with CKN samples. These genes are all involved in lignin biosynthesis ( Figure 9). Moreover, PAL and 4CL are not only involved in lignin biosynthesis but also in the formation of quercetin ( Figure 10). Therefore, these genes might be responsible for the 'sunshine muscat' grape russet phenotype. Caffeic acid, chlorogenic acid (3-O-caffeoylquinic acid), 4-O-p-coumaroylquinic acid, and ferulic acid, which are involved in lignin biosynthesis, were up-regulated in russet 'sunshine muscat' grapes (Table 6). Naringenin 7-O-glucoside and quercetin-related glycosides, which are involved in quercetin biosynthesis, were both up-regulated in russet 'sunshine muscat' grapes ( Table 6). Our results suggested that lignin and quercetin are the most important metabolites involved in the formation of 'sunshine muscat' grapes russet. The genes that encode 4CL, CAD, HCT, CCR, CCo AOMT, COMT, POD, and PAL were most related to russet formation, which suggests that the expression of these genes might induce the accumulation of metabolites related to 'sunshine muscat' grape russet. Thus, the results of the transcriptomic and metabolomic analysis were consistent with each other, which indicates that the synthesis of the quercetin and lignin likely underlies the formation of the 'sunshine muscat' grape russet. Lignin and quercetin are both products of the phenylpropane metabolic pathway. The expression levels of 4CL, CAD, HCT, CCR, CCo AOMT, and COMT were significantly higher in russet CKY samples compared with non-russet CKN samples. POD expression was observed to be both upregulated and down-regulated. In addition, PAL expression was significantly lower in CKY samples when compared with CKN samples. These genes are all involved in lignin biosynthesis ( Figure 9). Moreover, PAL and 4CL are not only involved in lignin biosynthesis but also in the formation of quercetin ( Figure 10). Therefore, these genes might be responsible for the 'sunshine muscat' grape russet phenotype. Caffeic acid, chlorogenic acid (3-O-caffeoylquinic acid), 4-O-p-coumaroylquinic acid, and ferulic acid, which are involved in lignin biosynthesis, were up-regulated in russet 'sunshine muscat' grapes (Table 6). Naringenin 7-O-glucoside and quercetin-related glycosides, which are involved in quercetin biosynthesis, were both up-regulated in russet 'sunshine muscat' grapes ( Table 6). Our results suggested that lignin and quercetin are the most important metabolites involved in the formation of 'sunshine muscat' grapes russet. The genes that encode 4CL, CAD, HCT, CCR, CCo AOMT, COMT, POD, and PAL were most related to russet formation, which suggests that the expression of these genes might induce the accumulation of metabolites related to 'sunshine muscat' grape russet. Thus, the results of the transcriptomic and metabolomic analysis were consistent with each other, which indicates that the synthesis of the quercetin and lignin likely underlies the formation of the 'sunshine muscat' grape russet.

Discussion
Thus far, researchers have studied pear and apple russet in detail. Research on 'Dangshansu pear' has shown that phenylpropane metabolic pathways are related to the formation of russet, and lignin biosynthesis can regulate russet formation [8,25]. Enzymes involved in lignin biosynthesis can be divided into three categories [26]: (1) enzymes in the phenylalanine metabolic pathway, e.g., PAL and 4CL, which show expression that is closely related to the total lignin, (2) enzymes related to the regulation of lignin monomer synthesis in lignin-specific synthesis pathways, e.g., COMT and CCo AOMT, (3) enzymes downstream of the lignin synthesis branch pathway and involved in the regulation of enzymes involved in the synthesis and polymerization of lignin monomers, including CCR, CAD, and POD, which are responsible for the ultimate reduction of various hydroxycinnamylcoenzyme A esters into lignin monomers. The monomers are polymerized into lignin, which play a key role in lignin synthesis and metabolism (Figure 9).
Studies on pears and apples have highlighted that lignin synthesis genes, such as COMT and C4H, have increased expression, underlying the lignin accumulation that leads to russet formation [5,27]. Studies on 'Dangshansu Pear' mutant varieties have shown that increased expression of CCoAOMT leads to an increase in the lignin content of the outer peel, which results in russet formation [28]. Similar results were observed in this study. We found that the expression of genes related to lignin synthesis was significantly different between CKY and CKN. Furthermore, we evaluated their expression by qRT-PCR. The expression patterns of these genes obtained from qRT-PCR were highly consistent with those in the RNA-seq data ( Figure 4). Thus, these genes may be considered the key genes conferring the formation of the 'sunshine muscat' grape russet. Hence, peel lignin synthesis appears to underlie the formation of 'sunshine muscat' grape russet ( Figure 9). However, the differences in lignin content and the expression of genes related to lignin synthesis during the rusting process in 'sunshine muscat' grape require further exploration.
Quercetin (3,5,7,3′,4′-pentahydroxyflavones) is an important aglycon in flavonoids, which are almost insoluble in water, and it is mostly present in the flowers, leaves, and fruit, among other tissues [29][30][31]. Quercetin is an important secondary metabolite whose stress resistance determines its distribution in grape fruits. It is mainly concentrated in the outer epidermal cells of fruits with lower content in the pulp and seeds [32]. Quercetin biosynthesis begins with coumaryl Co A and malonyl Co A, which both form chalcone under the catalysis of chalcone synthase (CHS), and chalcone synthesizes naringenin under the catalysis of chalcone isomerase (CHI). In turn, naringenin forms dihydro kaempferol under the action of flavonol-3-hydroxylase and then synthesizes dihydroquercetin under the action of flavonoid 3′hydroxylase. Lastly, dihydroquercetin synthesizes quercetin under the action of flavonol synthase [33] (Figure 10). Studies have shown that FaGT6 (DQ289587) can catalyze quercetin to form 3-O-glucosidethe in strawberries in addition to a small amount of 7-O-, 4′-O-, and 3′-O-glucosides and other glucosides [34]. Cp3GT (ACS15351) specifically catalyzes glucosylation at the 3-O position of quercetin in Citrus paradise [35]. In grape (Vitis vinifera), VvGT5 (AB499074) has quercetin 3-O-glucuronyltransferase activity, while VvGT6 (AB499075) uses quercetin as a substrate and has UDP-glucose and UDP-galactose glycotransferase activities [36].

Discussion
Thus far, researchers have studied pear and apple russet in detail. Research on 'Dangshansu pear' has shown that phenylpropane metabolic pathways are related to the formation of russet, and lignin biosynthesis can regulate russet formation [8,25]. Enzymes involved in lignin biosynthesis can be divided into three categories [26]: (1) enzymes in the phenylalanine metabolic pathway, e.g., PAL and 4CL, which show expression that is closely related to the total lignin, (2) enzymes related to the regulation of lignin monomer synthesis in lignin-specific synthesis pathways, e.g., COMT and CCo AOMT, (3) enzymes downstream of the lignin synthesis branch pathway and involved in the regulation of enzymes involved in the synthesis and polymerization of lignin monomers, including CCR, CAD, and POD, which are responsible for the ultimate reduction of various hydroxycinnamyl-coenzyme A esters into lignin monomers. The monomers are polymerized into lignin, which play a key role in lignin synthesis and metabolism (Figure 9).
Studies on pears and apples have highlighted that lignin synthesis genes, such as COMT and C4H, have increased expression, underlying the lignin accumulation that leads to russet formation [5,27]. Studies on 'Dangshansu Pear' mutant varieties have shown that increased expression of CCoAOMT leads to an increase in the lignin content of the outer peel, which results in russet formation [28]. Similar results were observed in this study. We found that the expression of genes related to lignin synthesis was significantly different between CKY and CKN. Furthermore, we evaluated their expression by qRT-PCR. The expression patterns of these genes obtained from qRT-PCR were highly consistent with those in the RNA-seq data ( Figure 4). Thus, these genes may be considered the key genes conferring the formation of the 'sunshine muscat' grape russet. Hence, peel lignin synthesis appears to underlie the formation of 'sunshine muscat' grape russet ( Figure 9). However, the differences in lignin content and the expression of genes related to lignin synthesis during the rusting process in 'sunshine muscat' grape require further exploration.
Quercetin (3,5,7,3 ,4 -pentahydroxyflavones) is an important aglycon in flavonoids, which are almost insoluble in water, and it is mostly present in the flowers, leaves, and fruit, among other tissues [29][30][31]. Quercetin is an important secondary metabolite whose stress resistance determines its distribution in grape fruits. It is mainly concentrated in the outer epidermal cells of fruits with lower content in the pulp and seeds [32]. Quercetin biosynthesis begins with coumaryl Co A and malonyl Co A, which both form chalcone under the catalysis of chalcone synthase (CHS), and chalcone synthesizes naringenin under the catalysis of chalcone isomerase (CHI). In turn, naringenin forms dihydro kaempferol under the action of flavonol-3-hydroxylase and then synthesizes dihydroquercetin under the action of flavonoid 3 hydroxylase. Lastly, dihydroquercetin synthesizes quercetin under the action of flavonol synthase [33] (Figure 10). Studies have shown that FaGT6 (DQ289587) can catalyze quercetin to form 3-O-glucosidethe in strawberries in addition to a small amount of 7-O-, 4 -O-, and 3 -O-glucosides and other glucosides [34]. Cp3GT (ACS15351) specifically catalyzes glucosylation at the 3-O position of quercetin in Citrus paradise [35]. In grape (Vitis vinifera), VvGT5 (AB499074) has quercetin 3-O-glucuronyltransferase activity, while VvGT6 (AB499075) uses quercetin as a substrate and has UDP-glucose and UDP-galactose glycotransferase activities [36].
Few reports have explored the relationship between quercetin and quercetin glycosides and the formation of russet. Previous studies on apples have found that the contents of quercetin, quercetin-3-rutin, quercetin-galactoside, quercetin-glucoside, quercetin-xylosin, quercetin-arabinoside, quercetin-rhamnoside, and quercetin-arabin-glucoside in the non-russet bud variety 'Feng Shuai' were significantly lower than those in the russet variety 'Jinguan', which indicates that the increase in quercetin and quercetin glycosides promotes the formation of 'Jinguan' apple russet [37]. The results of this study show that 4-coumarin-CoA ligase (100254698) and quercetin 3-O-methyltransferase (100260786) are up-regulated and phenylalanine ammonia-lyase (100241575) is down-regulated in russet 'sunshine muscat' grapes ( Table 4). The differentially expressed metabolites in the peels with and without russet are mainly phenolic substances, while these phenolic substances are mainly quercetin glycosides. Moreover, the content of quercetin glycosides in russet peel is significantly higher than that in non-russet peel ( Figure 6). Therefore, it can be hypothesized that, as quercetin content increases, glycosylation occurs under the action of UDP-glycosyltransferase to form different quercetin glycosylation products during the formation of 'sunshine muscat' grape russet. However, changes in quercetin and quercetin glycosides during the rusting process of 'sunshine muscat' grape require further exploration.

Conclusions
We found that the phenylpropane synthesis pathway is the key metabolic pathway associated with russet formation and that the regulation of genes related to lignin and quercetin synthesis differed, which promotes the formation of russet. In addition, phenols are the key metabolites underlying the formation of 'sunshine muscat' grape russet. Moreover, lignin and quercetin play major roles in phenol production. Thus, the synthesis of lignin and quercetin is responsible for the formation of berry russeting in 'sunshine muscat' grape.