Model-assisted analysis for tuning anthocyanin composition in grape berries

Abstract

Anthocyanin composition is responsible for the red colour of grape berries and wines, and contributes to their organoleptic quality. However, anthocyanin biosynthesis is under genetic, developmental and environmental regulation, making its targeted fine-tuning challenging. We constructed a mechanistic model to simulate the dynamics of anthocyanin composition throughout grape ripening in Vitis vinifera, employing a consensus anthocyanin biosynthesis pathway. The model was calibrated and validated using six datasets from eight cultivars and 37 growth conditions. Tuning the transformation and degradation parameters allowed us to accurately simulate the accumulation process of each individual anthocyanin under different environmental conditions. The model parameters were robust across environments for each genotype. The coefficients of determination (R²) for the simulated versus observed values for the six datasets ranged from 0.92 to 0.99, while the relative root mean square errors (RRMSEs) were between 16.8 and 42.1 %. The leave-one-out cross-validation for three datasets showed R² values of 0.99, 0.96 and 0.91, and RRMSE values of 28.8, 32.9 and 26.4 %, respectively, suggesting a high prediction quality of the model. Model analysis showed that the anthocyanin profiles of diverse genotypes are relatively stable in response to parameter perturbations. Virtual experiments further suggested that targeted anthocyanin profiles may be reached by manipulating a minimum of three parameters, in a genotype-dependent manner. This model presents a promising methodology for characterizing the temporal progression of anthocyanin composition, while also offering a logical foundation for bioengineering endeavours focused on precisely adjusting the anthocyanin composition of grapes.

INTRODUCTION

Anthocyanins are specialized metabolites essential for the coloration of plant organs, particularly fruits and flowers, whose colour is an important quality criterion for their visual appeal and market value (Zhang et al., 2014). Moreover, anthocyanins possess health-promoting properties, including antioxidant, anti-inflammatory and anti-cancer activities, making them a subject of interest in medicinal and nutraceutical research (Khoo et al., 2017). More than 600 anthocyanins have been identified to date, and they are synthesized from the precursor amino acid phenylalanine through the phenylpropanoid pathway (Zhang et al., 2014). This considerable diversity in anthocyanins is due largely to the biochemical decorations of a common chemical backbone via hydroxylation, methylation, glycosylation and acylation (Jaakola, 2013). These biochemical decorations have important influences on stability, bioavailability and colour hues, forming a nexus of structure and functionality for anthocyanins (Saito et al., 2013; Wen et al., 2020; Houghton et al., 2021). Therefore, developing plants with targeted anthocyanin profiles is of great interest for both scientific and application prospects.

The anthocyanin profile varies greatly among plant species and cultivars, mainly as a result of the activities and specificities of diverse decorating enzymes involved in the biosynthesis pathway, ultimately giving rise to a broad range of colour hues and patterns (Jaakola, 2013). In grape (Vitis vinifera) genotypes, anthocyanins include 3-monoglucosides, which can be further acylated, of five anthocyanidins (Mattivi et al., 2006; He et al., 2010), namely delphinidin (Dp), cyanidin (Cy), petunidin (Pt), peonidin (Pn) and malvidin (Mv), while pelargonidin is rarely present because of the substrate specificity of the enzyme DFR (dihydro flavonol 4-reductase) (Mattivi et al., 2006). Grape berry anthocyanins can be further categorized into di- or tri-hydroxylated anthocyanins based on their hydroxylation via the F3ʹH (flavonoid 3ʹ-hydroxylases) and F3ʹ5 ʹH (flavonoid 3ʹ,5ʹ-hydroxylases) enzymes (Falginella et al., 2010), their methylation by AOMT1 (anthocyanin O-methyltransferase 1) and AOMT2 (anthocyanin O-methyltransferase 2) (Hugueney et al., 2009; Fournier-Level et al., 2011), and their acylation by 3AT (anthocyanin 3-O-glucoside-6″-O-acyltransferase) (Rinaldo et al., 2015). The dominance or absence of a specific category of anthocyanins is genotype-dependent in grapevine (Mattivi et al., 2006). For example, most grape genotypes accumulate predominantly tri-hydroxylated anthocyanins (also known as Dp-derived anthocyanins), while genotypes such as Sangiovese synthesize dominantly di-hydroxylated anthocyanins (also known as Cy-derived anthocyanins) (Mattivi et al., 2006; Pastore et al., 2017). Some cultivars, including Pinot noir, only synthesize unacylated anthocyanins (Dimitrovska et al., 2011; Rinaldo et al., 2015). Moreover, the anthocyanins of most grape cultivars are limited to the berry skin, except in teinturier cultivars, which accumulate anthocyanins in both berry skin and pulp (Kong et al., 2021). These specificities in anthocyanin composition can be considered as a cultivar identification fingerprint (Mattivi et al., 2006; van Leeuwen et al., 2013).

In addition to the genetic determinants, the composition of anthocyanins within a given cultivar undergoes dynamic changes throughout berry development and is influenced by environmental factors such as temperature (Sugiura et al., 2018; de Rosas et al., 2022), water supply (Berdeja et al., 2014), nitrogen supply (Hilbert et al., 2003; Olsen et al., 2009) and light (Keller & Hrazdina, 1998). These complexities hinder the identification of key molecular regulators via association analysis between genotypes and phenotypes for the determination of anthocyanin composition (Costantini et al., 2015), because the traditional one-time-point phenotyping approaches cannot fully reflect the developmental and environmental dynamics of anthocyanin composition. Phenotyping time series across berry development may provide more comprehensive information for deciphering the role of genotype × development × environment interactions for anthocyanin composition. However, there has been limited research into the relative contributions of genotype, development and environment on the dynamics of anthocyanin composition.

Mathematical models for metabolic pathway analysis (Morgan & Rhodes, 2002; Baghalian et al., 2014) play an important role in experimental biology for advancing the boundaries of plant science (Marshallcolon et al., 2017). These models may help to dissect complex traits, such as sugar concentrations in fruits, into simple traits by integrating time-series of phenotyping data, consequently facilitating the association analysis between genotypes and phenotypes (Génard et al., 2010; Prudent et al., 2011). Moreover, these models can also provide a rational analysis of biochemical pathways, identify potential limiting metabolic steps, screen candidate intervention points for bioengineering, as well as generate novel hypotheses for testing by wet-lab experiments (Rios-Estepa et al., 2008; Farré et al., 2014; Wang et al., 2019, 2022). For example, several mathematical models have been developed for lignin biosynthesis (Lee & Voit, 2010; Faraji & Voit, 2017; Wang et al., 2018; Matthews et al., 2020, 2021). Faraji et al. (2018) investigated several computational models of lignin biosynthesis in various plant species, including black cottonwood (Populus trichocarpa), alfalfa (Medicago truncatula), switchgrass (Panicum virgatum) and the grass Brachypodium distachyon. Their findings indicated that the intermediates of the lignin heteropolymer biosynthetic pathway are similar, while the enzymatic reactions of the pathway exhibit significant variation, providing valuable clues for targeted bioengineering of lignin composition (Matthews et al., 2021). Recently, a qualitative metabolic model has been constructed for the metabolic pathway of anthocyanins (Wheeler & Smith, 2019). Wheeler et al. (2020) developed a kinetic model for the anthocyanin pathway, which explained the production of blue, purple and red anthocyanin pigments with multiple branches and substrate competition. This model was used to investigate the evolution of the anthocyanin pathway through fixed mutations and provided a theoretical framework to predict the consequences of new mutations in pigment phenotypes and pleiotropy (Wheeler & Smith, 2019). Notwithstanding its predictive value in assessing the effects of mutations on pigment production, this anthocyanin model provides mostly qualitative results and lacks quantitative outputs for a specific genotype over development under various environments.

The present study aimed to develop a dynamic model of anthocyanin composition based on the biosynthetic pathway (Boss et al., 1996), which should be robust for a wide spectrum of grape genotypes and environmental conditions. The Dynamic Anthocyanin Composition Model (DACM) introduced here was calibrated and validated through the utilization of observed individual anthocyanin concentrations in berries of diverse grape cultivars under varying conditions. Calibration was achieved using five publicly available datasets and one dataset that has not been published. A global sensitivity analysis was then conducted to identify the key parameters controlling the concentration of individual anthocyanins. Model-based virtual experiments were finally utilized to explore strategies for fine tuning the anthocyanin composition with specific targets of biochemical decorations, providing possible intervention points for reorienting the metabolic fluxes within the pathway.

MATERIALS AND METHODS

Model description

The DACM is a computational tool that enables the simulation of the accumulation profiles of individual anthocyanins during the development of a particular grape genotype from veraison (the onset of ripening) to maturity under varying environmental conditions (Hernández-Montes et al., 2021), using one set of genotype-specific parameters. The structure of the model (Fig. 1B) is based on an explicit anthocyanin biosynthesis pathway (Boss et al., 1996), which encompasses up to 15 different anthocyanins (3-glucosides) with various biochemical decorations of hydroxylation, methylation and acylation (Fig. 1A). Multiple enzymatic reactions are concatenated to simplify the model while preserving the crucial topology of the anthocyanin biosynthesis pathway (Fig. 1B).

The time-course changes in the quantity of a given anthocyanin were described as a result of three processes: biosynthesis, conversion to other anthocyanin forms and degradation to non-anthocyanin metabolites. The biosynthesis of one individual anthocyanin or the conversion of one individual anthocyanin (V_Cy, V_Dp) to another was modelled as the product of a relative rate constant (r_i, i from 1 to 13) multiplied by the substrate quantity, following the mass balance reaction principle. Similarly, the degradation of one anthocyanin to non-anthocyanin metabolites was modelled by multiplying the quantity of the specific anthocyanin with a relative degradation constant (kd). The degradation of anthocyanins is not well understood in grape berry, nor in other plant species, and the specific degradation rate for each individual anthocyanin is not available in the literature. For the sake of simplicity, we assumed that kd is the same for all degradation reactions, as previously proposed (Guardiola et al., 1995). To verify the robustness of this assumption, we also tested the model performance with different kd values together with model comparisons using the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) (Burnham and Anderson, 2002).

Taken together, the dynamic changes of each individual anthocyanin during berry development can then be described by the following 15 ordinary differential equations (ODEs):

where r_i (i from 1 to 13) are the relative transformation rate constants (g anthocyanin g whole berry^–1 d^–1); kd is the relative degradation constant (g g^–1 d^–1), V_Cy is the influx of cyanidin-based anthocyanins (g d^–1), V_Dp is the influx of delphinidin-based anthocyanins (g d–1). Cy_glc (g per berry) is the quantity of cyanindin-3-glucoside per berry; Dp_glc for delphinidin-3-glucoside; Pn_glc for peonidin-3-glucoside per berry; Pt_glc for petunidin-3-glucoside per berry; Mv_glc for malvidin-3-glucoside per berry; Pn_ac for peonidin-3-acetylglucoside per berry; Pn_cou for peonidin-3-coumaroyl glucoside per berry; Pt_ac for petunidin-3-acetylglucoside per berry; Pt_cou for petunidin-3-coumaroyl glucoside per berry; Mv_ac for malvidin-3-acetylglucoside per berry; Mv_cou for malvidin-3-coumaroyl glucoside per berry; Cy_ac for cyanidine-3-acetylglucoside per berry; Cy_cou for cyanidine-3-coumaroyl glucoside per berry; Dp_ac for delphinidin-3-acetylglucoside per berry; and Dp_cou for delphinidin-3-coumaroyl glucoside per berry. The quantity of each individual anthocyanin was expressed in equivalents of Mv_glc.

The input to the model was the total anthocyanin influx into the anthocyanin-specific pathway on each day (V_in, g d⁻¹), which was the sum of influxes into the two anthocyanin biosynthesis branches (V_Cy and V_Dp).

where δ is an allocation coefficient (from 0 to 1) between the two branches of the metabolic pathway.

V_in can be calculated with the following equations:

where $dTAobs/dt$ is the newly accumulated total anthocyanin per berry at each time step. It was estimated from the first-order derivative of the observed total anthocyanin (TA_obs) curves (Supplementary Data Fig. S1) in the skin or pulp. TA (g per berry) is the instantaneous total anthocyanins per berry in the simulation system. To take into account the developmental changes in anthocyanin degradation (Mori et al., 2007; Movahed et al., 2016), kd was assumed to change as a linear function of berry development (DOY: day of the year) with parameters a_kd and b_kd. The variable a_kd (g g d^–2) represents how fast the anthocyanin degradation constant changes as a function of berry development, while b_kd (g g⁻¹ d⁻¹) represents the basal degradation constant across the whole berry development (Table S1).

The integration of ODEs was utilized to determine the quantity of accumulated individual anthocyanins in each berry. The resulting amount of anthocyanin was then divided by the fresh weight of the grape berry or specific tissues (Supplementary Data Fig. S2), such as skin or pulp, to calculate the concentration of each individual anthocyanin.

Experimental data

Six datasets were utilized in this study obtained from six separate experiments encompassing eight distinct grapevine cultivars (Supplementary Data Table S2), two rootstocks, and 37 growth conditions that included variations in growing seasons, leaf-to-fruit ratio, water stress, nitrogen supply and light levels (Table 1). Datasets 1–5 were derived from previously published studies, and their respective experimental designs are briefly summarized below.

The first experiment investigated the effects of nitrogen (N) availability (0.34, 1.7 or 3.4 g N per plant as NH₄NO₃ applied at bloom) and light intensity (3 weeks at 100, 20 or 2 % of full sunlight, modulated outdoors with shade cloth starting from veraison) on growth and fruit ripening of pot-grown cv. Cabernet Sauvignon (V. vinifera L.) vines (Keller & Hrazdina, 1998). The second experiment contained three nitrogen treatments applied from fruit set to leaf fall (1.4, 3.6 and 7.2 mm N, denominated N1, N2, N3I, respectively) with cv. Merlot (V. vinifera) vines in a glasshouse (Hilbert et al., 2003). The third dataset illustrated the effects of different rootstocks and water supply levels on berry growth and anthocyanins in cv. Pinot noir (V. vinifera), which was grafted onto either rootstock 110R (drought tolerant, medium to high vigour) or 125AA (drought sensitive, high vigour) during three growing seasons in the field under control or water shortage conditions (Berdeja et al., 2014). The fourth experiment investigated the effects of two leaf-to-fruit ratio levels on the quality of berries of cvs. Cabernet Sauvignon and Sangiovese (V. vinifera), respectively (Bobeica et al., 2015). The fifth experiment showed the anthocyanins in skin and pulp of grape berries collected from cvs. Gamay, Gamay de Bouze and Gamay Freaux (V. vinifera) grown in a glasshouse (Kong et al., 2021). The sixth experiment was conducted in the current study. Briefly, two cultivars, Cabernet Sauvignon and Tempranillo (V. vinifera), were grown in a common garden vineyard named ‘Vitadapt’ (Suter et al., 2021) under standard viticultural practices. Berries were collected at 10-d intervals from veraison to maturity in two growing seasons with three biological replicates (30 berries of each replicate) to measure berry fresh weight with a high-precision balance and anthocyanin composition with high-performance liquid chromatography (HPLC) as described in Kong et al. (2021).

These datasets were chosen based on four criteria, including: (1) there must be measurements of individual anthocyanins in grape berries with specific analytical equipment (HPLC); (2) the dataset should contain different environmental treatments, as this model aims to simulate the accumulation of anthocyanins under varying environmental conditions; (3) the anthocyanin composition must be measured throughout berry ripening, covering at least four different developmental stages; and (4) the berry and/or skin fresh weight should be measured, enabling calculation of the total anthocyanin content per berry at each sampling date. Accordingly, the aforementioned six experiments all determined the developmental dynamics of grape growth (e.g. fresh weight, Supplementary Data Fig. S2) and quantified the anthocyanin composition from veraison to maturity with 4–10 sampling dates (Fig. S1). Moreover, the analysis of individual anthocyanins for all six datasets was performed by HPLC, and the quantification was based on calibration curves with malvidin-3-glucoside as an external standard, with all other individual anthocyanins being expressed as malvidin-3-glucoside equivalents.

Model resolution, parameterization and validation

The dynamic model was simulated with a 1-d time step and implemented using R software (R Core Team, 2013). The ODEs of the model were numerically integrated using the Euler method.

For parameterization, a genetic algorithm (GA) was used to estimate all the model parameters for each cultivar in the dataset (Supplementary Data Table S1). The GA function in R (R Core Team, 2013) was used to minimize an objective criterion, which was defined as:

where n is the number of conditions for each cultivar, m is the sampling number of each condition throughout berry development, and y_o and y_s are the observed and simulated values of anthocyanins, respectively.

The parameterization process with GA was repeated five times to obtain five sets of parameters to assess parameter stability. The values of the set of parameters that gave the smallest criterion value was used as the best estimated parameter values for each cultivar in each dataset (Supplementary Data Tables S1 and S3).

During the process of parameterizing the model for a specific cultivar under diverse environmental conditions it was discovered that the allocation coefficient (δ) varied in relation to the levels of nitrogen and light intensity. Subsequently, an equation was developed to illustrate the responses of δ to different nitrogen and light intensity conditions for the cultivar in question:

where a_δ and b_δ are the coefficients of nitrogen and light levels, respectively and c_δ is a constant. N_rel was the relative nitrogen availability under optimum growing condition (considered as 1) for the investigated cultivar, while light_rel was the relative light intensity under full sunlight (considered as 1).

To assess the performance and predictive ability of the model, various statistical measures were used for each cultivar, including the root mean squared error (RMSE), relative root mean squared error (RRMSE) and coefficient of determination for the linear correlation (R²). The model’s prediction quality was evaluated through the leave-one-out cross-validation method (Wallach et al., 2006). Specifically, this process was carried out for Cabernet Sauvignon (dataset 1), Pinot noir (dataset 2) and Merlot (dataset 3), which each contained more than three growth conditions. For each cultivar, a set of optimal parameters was estimated using all observed data except for one condition, and these parameters were then used to predict the anthocyanin composition under the omitted condition. This validation process was repeated multiple times, with nine runs for Cabernet Sauvignon, 12 runs for Pinot noir and three runs for Merlot. The RMSE and RRMSE were calculated for each validation condition, and averaging these values across all growing conditions provided an overall estimate of prediction quality for each cultivar (Wallach et al., 2014).

Global sensitivity analysis

To determine the key metabolic steps influencing anthocyanin composition, global sensitivity analysis using the Morris method (Morris, 1991) was performed. The investigation involved the use of 8–18 parameters for different anthocyanin compositions to study their effects on model outputs. The assumption was made that all the investigated parameters were uniformly distributed within the range of 0.9–1.1 times of the default values. The sensitivity index (SI) and its standard deviation (ST) obtained from the Morris method were used to assess parameter sensitivity. The stability of the sensitivity analysis was ensured by verifying the convergence of sensitivity ranking based on the SI of each parameter by gradually increasing the number of trajectories and, finally, 300 trajectories were applied to provide a stable sensitivity ranking. The sensitivity analysis was performed using the ‘Morris’ function in the ‘sensitivity’ package of the R language (Iooss & Lemaître, 2015).

RESULTS

Overview of the influence of genotype and environment on total anthocyanins and anthocyanin composition

The concentrations of total anthocyanins showed clear differences between cultivars and growing seasons and were increased by high light intensity (Supplementary Data Fig. S1A) and leaf-to-fruit ratio (Fig. S1D), but were decreased by high N levels (Fig. S1A, B) and water supply (Fig. S1C). Here, we focused on the responses of anthocyanin composition related to their biochemical decorations of hydroxylation (Fig. 2), methylation (Fig. S3) and acylation (Fig. S4).

Based on their hydroxylation, individual anthocyanins can be grouped into di- and tri-hydroxylated anthocyanins (Fig. 2). The ratio of di- to tri-hydroxylated anthocyanins (RDT) represents one of the most important properties for anthocyanin composition, indicating whether the grape colour is reddish or dark blue (Castellarin et al., 2006). This ratio was largely determined by genotypes: the Sangiovese and the pulp of Gamay de Bouze and Gamay Freaux had RDT > 1 (Fig. 2D and E), while the remaining genotypes, including Cabernet Sauvignon, Gamay, Merlot, Pinot noir and Tempranillo, had RDT < 1 (Fig. 2A–C and F). These results showed that Sangiovese skin and the pulp of Gamay de Bouze and Gamay Freaux accumulated dominantly di-hydroxylated anthocyanins, while the remaining genotypes accumulated dominantly tri-hydroxylated anthocyanins. Over berry development, the RDT gradually increased in genotypes with RDT > 1 (Fig. 2D and E), while it gradually decreased in genotypes with RDT < 1 (Fig. 2A–C and F). However, the RDT varied less in response to the explored combinations of growth conditions (Fig. 2) than the concentrations of total anthocyanins (Supplementary Data Fig. S1). The RDT in Cabernet Sauvignon was not significantly different between light and nitrogen conditions, except the extreme condition under the highest light intensity (L100) and lowest nitrogen supply (N1) (Fig. 2A). The RDT of the extreme condition (L100-N1) was >1 while in other conditions it was <1 (Fig. 2A). Moreover, the effects of nitrogen and light intensity on RDT were opposite, as the RDT responded positively to increasing light intensity and negatively to increasing N levels (Fig. 2A). A similar negative effect of nitrogen on the RDT was also observed in Merlot (Fig. 2B). The RDT of Pinot noir was hardly affected by rootstock, water stress or growing seasons, particularly around maturity (Fig. 2C). The RDT was affected by the leaf-to-fruit ratio in a genotype-dependent manner, with the RDT of Cabernet Sauvignon being systematically decreased by lower leaf-to-fruit ratio while the RDT of Sangiovese was hardly affected by the leaf-to-fruit ratio (Fig. 2D). The developmental profiles of RDT were different in 2013 and 2014 for Tempranillo and Cabernet Sauvignon; however, they overlapped when resynchronized with veraison dates (Fig. 2F).

Based on the methylation of individual anthocyanins, they could be grouped into methylated and unmethylated anthocyanins (Supplementary Data Fig. S3). Seven out of the eight investigated cultivars possessed predominantly methylated anthocyanins, with a ratio of methylated to unmethylated anthocyanins ranging between 2 and 200 (Fig. S3), while only Sangiovese had more unmethylated anthocyanins at maturity with a methylated to unmethylated ratio <1 (Fig. S3D). In response to various growing conditions, this ratio was increased by conditions that reduced the total anthocyanins (Fig. S3), such as low light (Fig. S3A), high nitrogen (Fig. S3A and S3B) and low leaf-to-fruit ratio for Cabernet Sauvignon (Fig. S3D).

Based on their acylation, individual anthocyanins could be grouped into acylated and unacylated anthocyanins (Supplementary Data Fig. S4). Pinot noir and Sangiovese do not accumulate acylated anthocyanins (Mattivi et al., 2006; Rinaldo et al., 2015) and thus showed a ratio of acylated to unacylated anthocyanins of 0 (Fig. S4C and S4D). For the other cultivars, this ratio ranged from 0.05 to 1.0 (Fig. S4B, S4D, S4E and S4F), with higher values under conditions that reduced the total anthocyanins, such as high nitrogen (Fig. S4B) and low leaf-to-fruit ratio (Fig. S4D).

Calibration of the DACM

To simulate the developmental dynamics of anthocyanin composition in different cultivars under various growth conditions, the DACM was developed based on mass balance of biochemical reaction rules with the total anthocyanin influx as input and reaction rates as parameters. The DACM was calibrated with actual measurements of each individual anthocyanin along berry development from the six datasets (Table 1). With a unique set of parameter values for each cultivar in each dataset (Supplementary Data Table S3), the model simulations were highly aligned with the observed results for most genotype × growth conditions (Figs 3 and 4; Figs S5–S12). The performance of the DACM for each dataset is briefly described below.

For the dataset collected during the current study (dataset 6), the DACM precisely simulated the developmental dynamics of the observed concentrations of 11 individual anthocyanins in Cabernet Sauvignon in 2013 and 2014 with RRMSEs of 25.2 and 25.6 %, respectively, and R² of 0.96 in both years (Fig. 3). The model also accurately reproduced the developmental dynamics of the observed concentrations of nine individual anthocyanins in Tempranillo in both 2013 and 2014, with RRMSE of 23.7 and 17.0 %, and R² of 0.96 and 0.99, respectively (Fig. 3). These results were obtained with the same set of parameter values for a given cultivar in both years (Supplementary Data Table S3), indicating the stability of parameters across growing seasons.

Moreover, we tested the capability of the model to simulate different anthocyanin categories according to their molecular decorations, including di- vs tri-hydroxylated, and methylated vs unmethylated anthocyanins (Fig. 4). The model precisely simulated the developmental dynamics of the observed concentrations of the four anthocyanin categories in Cabernet Sauvignon and Tempranillo in both 2013 and 2014, with RRMSEs ranging from 10.2 to 15.4 % and R² ranging from 0.98 to 1. These results show that the model accuracy is higher for anthocyanin categories with distinct decorations (Fig. 4) than for individual anthocyanins (Fig. 3).

Similarly, the model simulation agreed well with experimental observations for Cabernet Sauvignon under nine growth conditions covering three canopy light intensities and three levels of soil nitrogen supply (dataset 1, Supplementary Data Figs S5 and S6) with RRMSE ranging from 12.8 to 44.9 % and R² ranging from 0.82 to 0.97. In detail, the model simulation performed better in the three N conditions under high (light 3) and moderate light (light 2) with RRMSE ranging from 16.9 to 23.5 %, than in the low light conditions (light 1) with RRMSE ranging from 28.8 to 44.9 % (Fig. S5, Keller). The larger RRMSE under low light intensity (light 1) was attributed mainly to the underestimation of the malvidin-3-glucoside (Mv_glc) concentration in the rapid accumulation stage, when their counterparts under higher light (light 2 and light 3) had already plateaued. On the other hand, the model performed similarly under the three N levels under a given light intensity in dataset 1 (Fig. S5, Keller). The second dataset (Figs S5 and S6, Hilbert) consisted of three soil nitrogen supply levels in Merlot with 13 individual anthocyanins, and the model precisely reproduced the developmental dynamics of individual anthocyanins or anthocyanin categories with RRMSE ranging from 18.4 to 28.8 % and R² ranging from 0.95 to 0.98.

The third dataset (Supplementary Data Figs S7 and S8) consisted of two rootstocks and two water supply conditions for Pinot noir in three growing seasons (2009–2011), which contained five individual anthocyanins with RRMSE for the simulation results ranging from 8.0 to 38.7 %, and R² ranging from 0.82 to 1.00.

The fourth dataset (Supplementary Data Figs S9 and S10) consisted of two leaf-to-fruit ratios for Cabernet Sauvignon and Sangiovese, which contained 11 and five individual anthocyanins respectively with RRMSE for the simulation results ranging from 13.4 to 61.2 %, and R² ranging from 0.86 to 0.97. Under low leaf-to-fruit ratio (three leaves per cluster), the large RRMSE for Cabernet Sauvignon was attributed mainly to the model overestimating delphinidin-3-glucoside (Dp_glc) and peonidin-3-glucoside (Pn_glc) during berry development (Fig. S8, 3L CS), while the large RRMSE for Sangiovese was attributed mainly to the model underestimating cyanidin-3-glucoside (Cy_glc) in the late berry development stages (Fig. S8, 3L S).

The fifth dataset (Supplementary Data Figs S11 and S12) consisted of three cultivars (Gamay, Gamay de Bouze and Gamay Freaux) in two berry tissues (skin and pulp) which contained seven individual anthocyanins with RRMSE for the simulation results ranging from 15.4 to 50.1 %, and R² ranging from 0.94 to 0.98. The large RRMSE for pulp was attributed mainly to the model underestimating peonidin-3-glucoside (Pn_glc) and overestimating malvidin-3-glucoside (Mv_glc) during berry development (Fig. S11, pulp). Because of the close genetic relationship among the three cultivars (Kong et al., 2021), the anthocyanins in the skin were simulated with the same set of parameter values for all cultivars (Table S3), while those in the pulp were simulated with a set of parameter values different from the skin (Table S3). These results indicated that the skin and pulp of teinturier cultivars (i.e. cultivars with anthocyanins in both the skin and pulp) need to be considered separately.

Validation of the dynamic anthocyanin composition model

To assess the prediction quality of the model, three datasets with more than three environmental conditions for each dataset (datasets 1, 2 and 3) were tested with the leave-one-out cross-validation (Wallach et al., 2006). The model prediction quality was high under the nine combinations of light and N levels in Cabernet Sauvignon, with RRMSEP = 26.38 % and mean R² = 0.91 (Table 2). Similarly, the model prediction quality was high under three N levels in Merlot, with mean RRMSEP = 32.91 % and mean R² = 0.96 (Table 2). The model prediction quality was also high under the 12 combinations of water stress and rootstocks in different vintages in Pinot noir, with mean RRMSEP = 28.80 % and mean R²` = 0.99 (Table 2). Overall, the leave-one-out cross-validation results (RMSEP, RRMSEP) were comparable to model calibration results (RMSE, RRMSE), indicating the DACM possessed very high prediction quality.

Sensitivity analysis of the model

To investigate the influence of different parameters on model outputs, a global sensitivity analysis was conducted using the Morris method (Morris, 1991). The model parameters were categorized into three groups based on their association with conversion rate (r_i), degradation rate (kd) or allocation coefficient (δ). Four cultivars with distinct anthocyanin composition, namely Pinot noir, Sangiovese, the skin of Gamay Freaux and Cabernet Sauvignon from dataset 4, were chosen for the sensitivity analysis (Supplementary Data Table S2). The results indicated that the parameter sensitivity for any given anthocyanin remained relatively stable throughout berry development under the growing conditions investigated (Figs S13–S16). Therefore, only the most sensitive parameters for the concentrations of each individual anthocyanin at maturity were examined in detail (Fig. 5; Figs S17–S20).

Both Pinot noir and Sangiovese berries accumulate exclusively unacylated anthocyanins (Supplementary Data Fig. S4C and S4D; Mattivi et al., 2006; Rinaldo et al., 2015), and therefore their berries biosynthesize only five individual anthocyanins that were modelled here. Moreover, Pinot noir predominately accumulates tri-hydroxylated anthocyanins (RDT < 1, Fig. 1C) and Sangiovese predominately accumulates di-hydroxylated anthocyanins (RDT > 1, Fig. 1D; Table S2; Mattivi et al., 2006). Therefore, we compared parameter sensitivities in these two cultivars for the five individual anthocyanins (Fig. 5).

In Pinot noir (Fig. 5A; Supplementary Data Fig. S17), degradation (a_kd, b_kd) had the greatest effect on the two unmethylated anthocyanins (Cy_glc and Dp_glc), with a smaller effect from influx to the pathway (c_δ). The concentration of mono-methylated Pn_glc was mainly affected by biosynthesis (r₁), influx to the pathway (c_δ) and degradation (a_kd, b_kd). The concentration of mono-methylated Pt_glc was primarily affected by competition with Mv_glc biosynthesis (r₃) from the common substrate Dp_glc, followed by degradation (a_kd, b_kd), influx to the pathway (c_δ) and biosynthesis (r₂). The concentration of di-methylated and usually dominant Mv_glc was primarily affected by degradation (a_kd, b_kd), followed by influx to the pathway (c_δ), competition (r₂) and biosynthesis (r₃). The ratio of di- to tri-hydroxylated anthocyanins was only affected by influx to the pathway (c_δ), while the ratios of methylated to unmethylated anthocyanins were mainly affected by degradation (a_kd, b_kd), followed by r₃ and r₁.

In Sangiovese (Fig. 5B; Supplementary Data Fig. S18), all five individual anthocyanins were mainly affected by degradation (a_kd, b_kd) and influx to the pathway (c_δ), while their biosynthesis, conversion or competition played minor roles. The ratio of anthocyanins with two different decorations showed the same model parameter sensitivity as in Pinot noir.

In the berry skin of Gamay Freaux (Fig. 5C; Supplementary Data Fig. S19), degradation (a_kd, b_kd) had the greatest effect on the two unmethylated anthocyanins (Cy_glc and Dp_glc), with a smaller effect from biosynthesis (c_δ). The concentration of methylated and three unacylated anthocyanins (Pn_glc, Pt_glc, Mv_glc) was mainly affected by influx to the pathway (c_δ), followed by biosynthesis and competition. The concentration of acylated anthocyanins was mainly affected by degradation (a_kd, b_kd), followed by competition and biosynthesis. The ratio of di- to tri-hydroxylated anthocyanins was only affected by influx to the pathway (c_δ), while the ratios of methylated to unmethylated anthocyanins and acylated to unacylated anthocyanins were mainly affected by degradation (a_kd, b_kd).

In Cabernet Sauvignon (Fig. 5D; Supplementary Data Fig. S20), all individual anthocyanins were mostly influenced by their biosynthesis, followed by degradation (a_kd, b_kd) and competition. The ratio of the three differently decorated anthocyanins showed the same model parameter sensitivity as in Gamay Freaux.

Virtual experiment for targeted tuning of anthocyanin composition

A virtual simulation (Fig. 6) was conducted to explore potential strategies aiming at tuning the anthocyanin composition for targeted objectives, which were set to increase the proportions of tri-hydroxylated and methylated anthocyanins that are more stable and provide specific colour hues (He et al., 2010; Houghton et al., 2021; Liu et al., 2022). The virtual experiment was implemented in Sangiovese and Cabernet Sauvignon (Supplementary Data Table S2) by modifying the top seven most sensitive parameters identified via the global sensitivity analysis (Fig. 5). Modulating these seven parameters around their default values by −10, −5, 0, +5 and + 10 % while keeping other initial, input and parameter values in the default condition resulted in 78 125 possible combinations. As expected from the global sensitivity analysis, the proportions of tri-hydroxylated anthocyanins were exclusively affected by the parameter c_δ (Fig. 5A, B), which reflects the influx at the branching point of the anthocyanin metabolism pathway (Fig. 1B). In contrast, the proportions of methylated anthocyanins were influenced by all seven parameters with different intensities (Fig. 6A, B). Interestingly, the same variation range (±10 %) in c_δ led to different magnitudes of changes in the tri-hydroxylated anthocyanins in Sangiovese and Cabernet Sauvignon, with variations of ±22 % in Sangiovese but only ±1.3 % in Cabernet Sauvignon (Fig. 6A, B). Similar differences were also observed in the proportion of methylated anthocyanins, where Sangiovese showed a greater variation range (–25 to ~55 %) than those in Cabernet Sauvignon (only about ±6 %) (Fig. 6A, B).

The global optimal combinations for the highest proportions of tri-hydroxylated and methylated anthocyanins were obtained when all seven parameters were simultaneously set to their largest boundaries (namely ±10 %, max7 in Fig. 6E, F). These parameter combinations resulted in a 22 % increase in the proportion of tri-hydroxylated anthocyanins and 58 % increase in the proportion of methylated anthocyanins in Sangiovese, while the increases were only 1.3 and 5.5 %, respectively, in Cabernet Sauvignon (Fig. 6A, B). However, modifying seven parameters that may represent seven enzymatic steps will be almost unfeasible via current bioengineering technologies, which are more suitable to handle 1–3 genes or enzymes (Noda et al., 2017; Zhu et al., 2021). Therefore, we further explored what might be the minimum number of adjusted parameters to reach about 90 % of the improvements in the global optimal combinations (Fig. 6C, D). It was found that the targets were approached when changing 1–3 parameters, after which further increasing the number of adjusted parameters brought minor improvements. In fact, changing only three parameters (max3) yielded a performance of 92.5 % of the global optimal combination (max7) in Sangiovese and 79.4 % in Cabernet Sauvignon (Fig. 6E, F). Moreover, these local optimal combinations (max3) were reached by adjusting different sets of three parameters in Sangiovese and Cabernet Sauvignon, namely c_δ (−10 %), a_kd (−10 %) and b_kd (−10 %) for Sangiovese, and c_δ (−10 %), a_kd (−10 %) and r₃ (+10 %) for Cabernet Sauvignon, respectively (Fig. 6E, F).

We then verified which/how individual anthocyanins were modified for increasing the proportions of tri-hydroxylated and methylated anthocyanins in the simulations. Under different parameter combinations, the changes in proportions of individual anthocyanins in Sangiovese were mainly attributed to a significant decrease in the unmethylated anthocyanin Cy_glc, while its downstream product Pn_glc, which is a methylated anthocyanin, increased significantly (Fig. 6G). In Cabernet Sauvignon, which has a broader spectrum of individual anthocyanins, the changes were relatively small compared to Sangiovese. Specifically, the contents of the unmethylated anthocyanins Cy_glc and Dp_glc, and the methylated anthocyanin Pt_glc, all decreased slightly, while the methylated anthocyanins Pn_glc and Mv_glc and their downstream products all increased slightly in Cabernet Sauvignon (Fig. 6H).

DISCUSSION

The present study parameterized and tested a mechanistic model that accurately simulated the dynamic accumulation of individual anthocyanins in ripening grape berries. To calibrate and validate the model, six datasets with eight V. vinifera cultivars and 37 environmental conditions were utilized. The model parameters remained consistent and reliable across varying environmental conditions within each dataset for a given cultivar. This genotype-dependent but environment-independent property of model parameters enables the current model to serve as a novel phenotyping tool to dissect the complex traits of dynamic anthocyanin profiles into simple traits (Bertin et al., 2010; Génard et al., 2010; Dai et al., 2017). Combining modelled traits with QTL/GWAS analysis could help us to unveil the metabolic steps responsible for fine-tuning anthocyanin composition in grape berries and, eventually, other crops. Such a tool would facilitate genotype-to-phenotype analysis and prediction (Chenu et al., 2018).

The biosynthesis and subsequent decorations of anthocyanins, as for most specialized metabolites, are arranged in metabolic pathways with various topological structures, including linear, cyclical, branched or 3-D grids (Farré et al., 2014). For example, the anthocyanin decorations in Arabidopsis seem to be arranged in a highly branched 3-D grid (Saito et al., 2013), while the anthocyanin decorations follow strict orders in a linear way with several metabolic branches in Petunia spp. (Provenzano et al., 2014), as well as in grapes (V. vinifera) (Fig. 1; Ford et al., 1998; Hugueney et al., 2009). As a result, attempts to modify anthocyanin composition, as well as other specialized metabolites, often suffer from high uncertainties due to pathway complexity (Zhang et al., 2014; Noda et al., 2017; Zhu et al., 2021; Wang et al., 2022). Rational analysis of a metabolic pathway with mathematic models may aid in overcoming such difficulties by evaluating flux distributions and identifying candidate intervention points for bioengineering in order to enrich desirable compounds (Farré et al., 2014; Faraji & Voit, 2017; Wang et al., 2019). To this end, we first developed a dynamic anthocyanin composition model (DACM) and then used it in a virtual experiment to explore possible strategies to enrich tri-hydroxylated and methylated anthocyanins in different genetic backgrounds (Fig. 6). The variation in the relative proportions of di- to tri-hydroxylated anthocyanins was primarily caused by the allocation coefficient (δ) between the two anthocyanin biosynthesis branches, in agreement with the suggestion that anthocyanin hydroxylation is mainly determined by the relative expression of VvF3ʹH and VvF3ʹ5 ʹH in grapevine (Castellarin et al., 2007). On the other hand, fine-tuning methylated anthocyanins might be more challenging, because the proportion of methylated anthocyanins exhibited more complex reactions in a genetic background-dependent manner. This seems mainly due to the higher number of parameters regulating anthocyanin methylation than hydroxylation and due to the fact that the cumulative effects of parameters produce greater changes for the ratio of methylated-to-unmethylated anthocyanins. Interestingly, genotypes with complex anthocyanin profiles, such as Cabernet Sauvignon with at least 14 individual anthocyanins, exhibited greater composition stability in response to parameter perturbations than those with simple anthocyanin profiles, such as Sangiovese, with five individual anthocyanins. These results highlight the trade-offs that may occur in complex anthocyanin biosynthesis networks, because of the competition for shared substrates and/or enzymes (Wheeler & Smith, 2019). With these trade-offs, the shifts towards producing one type of pigment can result in a reduction in the production of other pigments due to limited substrates or modified substrate specificity of enzymes (Wheeler and Smith, 2019). Within a more complex anthocyanin biosynthesis network, the number of metabolic steps/branches with competing substrates and/or enzymes will increase, and consequently increase the probabilities to mitigate the effects of exogenous perturbations. Keeping this in mind, manipulating the anthocyanin composition to a specific target through genetic engineering may be easier in cultivars with simpler anthocyanin profiles (Zhu et al., 2021). These cultivars can offer reduced regulatory complexity, making them easier to achieve desired alterations without unintended effects, and provide more predictability for the outcomes (Shimada et al., 2001; Lin-Wang et al., 2014).

In theory, the enrichment of a desirable metabolite in a pathway may be reached by increasing its biosynthesis precursor and catalytic enzyme activities while decreasing its degradation and/or competitive pathway (Farré et al., 2014; Manela et al., 2015; Wang et al., 2021). However, it is not always straightforward to predict which steps or combinations of steps are the most pertinent strategy to reconstruct or reorient a metabolic pathway for producing targeted products. Our virtual experiments showed that the target, which was set to simultaneously increase the proportions of tri-hydroxylated and methylated anthocyanins, was optimally achieved by collectively modifying up to seven parameters. Considering the feasibility in most circumstances for plant-based bioengineering, we further explored the minimal set of parameters required to reach about 90 % of the optimal achievement. We found that three parameters will largely fulfil the target, but the exact three parameters differed between the two tested genotypes, namely Cabernet Sauvignon and Sangiovese. For Cabernet Sauvignon, the best parameter combination constituted the influx at branching points of the pathway and the degradation of anthocyanins. For Sangiovese, the two first parameters were the same as in Cabernet Sauvignon, but the third parameter (r₃) was related to reducing a competitive branch. Interestingly, the parameter r₃ may represent the step of anthocyanin methylation, primarily involving two O-methyltransferases (AOMTs) in grape berries. The VvAOMT1 is known to preferentially catalyse the 3ʹ,5ʹ methylation (Hugueney et al., 2009), while VvAOMT2 prefers 3ʹ methylation (Fournier-Level et al., 2011). In a study mimicking an increase in r₃ through the over-expression of VvAOMT1 in Petunia spp., a higher percentage of methylated anthocyanins was observed (Provenzano et al., 2014), in agreement with model predictions. These results offer a significant direction for the development of more efficient strategies for fine tuning anthocyanin composition, and highlight again the importance of the genetic background of the host plants.

Both our model analysis and virtual experiments showed that the anthocyanin degradation-related parameters (a_kd, b_kd) played an important role in the model for different anthocyanin compositions (Fig. 6). Similarly, previous studies have suggested that reducing the in vivo process of anthocyanin degradation could increase crop pigmentation and prevent colour degradation (Oren-Shamir, 2009; Zipor et al., 2015; Liu et al., 2022). In fact, anthocyanin degradation can be induced by spontaneous reactions, enzymatic activity or both (Oren-Shamir, 2009). Peroxidases were discovered to be involved in anthocyanin degradation in grape (Calderon et al., 1992) and Brunfelsia calycina flowers (Zipor et al., 2015). Moreover, the VviPrx31 peroxidase may act as a candidate gene involved in anthocyanin degradation in ripening grape berries under high temperature (Movahed et al., 2016). However, the nature of anthocyanin degradation in grape and other plants is far from being fully understood. To further explore this point, we tested the model performance using two sets of kd values for the di- and tri-hydroxylated anthocyanins. The results showed that the models with two sets of kd values generally provided slightly better reproductions of the observations (mean RRMSE = 29.85 %, Supplementary Data Table S4) than models with one set of kd values (mean RRMSE = 31.7 %, Table S4). However, when simultaneously considering model precision and model complexity, the two-kd model performed worse than the one-kd model (Table S4), as indicated by the AIC and BIC (Burnham and Anderson, 2002). Following Occam’s razor principle (Forster, 2000), we opted for a simpler model with a single set of kd. Though this point remains to be tested experimentally, this analysis suggests that the degradation constant might be similar for distinct anthocyanins and the one-kd assumption seems reliable. Moreover, the only study that conducted substrate specificity analysis for one anthocyanin degradation enzyme, BcPrx01 from Brunfelsia calycina, also showed that the rates of peroxidase-catalysed degradation were similar for all individual anthocyanins (Zipor et al., 2015). These findings underline the complexity of anthocyanin degradation and signal the need for further investigation into this process in grapes and other plants.

Despite the robustness of our dynamic anthocyanin composition model, the model may not fully account for alterations in anthocyanin composition under extreme treatment conditions, such as dark or severe nutrient stress. These severe stresses may result in significant variability in the pattern of anthocyanin composition across different cultivars. This suggests that the model parameter values may vary with the extreme environmental perturbations and that this type of response should be studied to improve the model prediction. Moreover, given the laborious measurement of anthocyanin concentrations, which serve as the input for the current model, further improvements to model usability are currently being explored. In particular, the development of a total anthocyanin prediction model, which incorporates cultivar, sugars, light, temperature (Sugiura et al., 2018) and other factors as inputs, may enhance the predictive power of the model. Overall, the model holds promise as a useful tool in phenotyping time-series of anthocyanin measurements, as well as a rationale for bioengineering applications aiming to fine tune anthocyanin composition.

SUPPLEMENTARY DATA

Supplementary data are available online at https://dbpia.nl.go.kr/aob and consist of the following.

Fig. S1: Accumulation of total anthocyanins in various grape genotypes and environments. Fig. S2: Berry fresh weight in various grape genotypes and environments. Fig. S3: Ratio of methylated to unmethylated anthocyanin dynamics in various grape genotypes and environments. Fig. S4: Ratio of acylated to unacylated anthocyanin dynamics in various grape genotypes and 56 environments. Fig. S5: Comparison between observed and simulated concentrations of various individual anthocyanins in response to the combinations of light and N treatments (dataset 1, Keller) or solely N supplies (dataset 2, Hilbert). Fig. S6: Comparison between observed and simulated concentrations of four anthocyanin types. Fig. S7: Comparison between observed and simulated concentrations of various individual anthocyanins in response to the combinations of rootstock and water supply treatments in different vintages (dataset 3). Fig. S8: Comparison between observed and simulated concentrations of four anthocyanin types. Fig. S9: Comparison between observed and simulated concentrations of various individual anthocyanins in responses to the combinations of leaf-to-fruit ratio levels and cultivars (dataset 4, Bobeica). Fig. S10: Comparison between observed and simulated concentrations of four anthocyanin types. Fig. S11: Comparison between observed and simulated concentrations of various individual anthocyanins in response to the combinations of tissues and cultivars (dataset 5, Kong). Fig. S12: Comparison between observed and simulated concentrations of four anthocyanin types. Fig. S13: Parameter sensitivity for anthocyanin compositions during berry development in Pinot noir. Fig. S14: Parameter sensitivity for anthocyanin compositions during berry development in Sangiovese. Fig. S15: Parameter sensitivity for anthocyanin compositions during berry development in Gamay Freaux skin. Fig. S16: Parameter sensitivity for anthocyanin compositions during berry development in Cabernet Sauvignon. Fig. S17: The most sensitive parameters for each individual anthocyanin at maturity for Pinot noir. Fig. S18: The most sensitive parameters for each individual anthocyanin at maturity for Sangiovese. Fig. S19: The most sensitive parameters for each individual anthocyanin at maturity for Gamay Freaux skin. Fig. S20: The most sensitive parameters for each individual anthocyanin at maturity for Cabernet Sauvignon. Table S1: List of the parameters used in the model. Table S2: Characteristics of anthocyanin compositions in different cultivars. Table S3: List of model parameter values of each condition in all the datasets. Table S4: Evaluation of model goodness-of-fit for model comparison.

ACKNOWLEDGMENTS

This research was supported partly by the National Key R&D Program of China (2021YFE0109500), National Natural Science Foundation of China (U20A2041), CAS Youth Interdisciplinary Team (JCTD-2022-06), Agricultural Breeding Project of Ningxia Hui Autonomous Region (NXNYYZ202101) and CAS Youth Interdisciplinary Team (JCTD-2022-06). Research was conducted as part of the LIA INNOGRAPE International Associated Laboratory.

AUTHOR CONTRIBUTIONS

The study was designed by YW and ZD; BS, MG, GH, SD, EG, SP, MK, CR and ZD contributed to the dataset collection; YW, JC and ZD constructed the model and wrote the simulation code; YW undertook the model testing and refinement. YW, BS and ZD contributed to analysing the simulation results; YW wrote the draft and all authors contributed to revising the paper; all authors approved the final manuscript.

LITERATURE CITED