Modeling the dependence of respiration and photosynthesis upon light, acetate, carbon dioxide, nitrate and ammonium in Chlamydomonas reinhardtii using design of experiments and multiple regression.

BACKGROUND: In photosynthetic organisms, the influence of light, carbon and inorganic nitrogen sources on the cellular bioenergetics has extensively been studied independently, but little information is available on the cumulative effects of these factors. Here, sequential statistical analyses based on design of experiments (DOE) coupled to standard least squares multiple regression have been undertaken to model the dependence of respiratory and photosynthetic responses (assessed by oxymetric and chlorophyll fluorescence measurements) upon the concomitant modulation of light intensity as well as acetate, CO₂, nitrate and ammonium concentrations in the culture medium of Chlamydomonas reinhardtii. The main goals of these analyses were to explain response variability (i.e. bioenergetic plasticity) and to characterize quantitatively the influence of the major explanatory factor(s).
RESULTS: For each response, 2 successive rounds of multiple regression coupled to one-way ANOVA F-tests have been undertaken to select the major explanatory factor(s) (1st-round) and mathematically simulate their influence (2nd-round). These analyses reveal that a maximal number of 3 environmental factors over 5 is sufficient to explain most of the response variability, and interestingly highlight quadratic effects and second-order interactions in some cases. In parallel, the predictive ability of the 2nd-round models has also been investigated by k-fold cross-validation and experimental validation tests on new random combinations of factors. These validation procedures tend to indicate that the 2nd-round models can also be used to predict the responses with an inherent deviation quantified by the analytical error of the models.
CONCLUSIONS: Altogether, the results of the 2 rounds of modeling provide an overview of the bioenergetic adaptations of C. reinhardtii to changing environmental conditions and point out promising tracks for future in-depth investigations of the molecular mechanisms underlying the present observations.

Background

In plants and <span class="Species">algae, energy transduction processes involve the respiratory and photosynthetic electron transpn>ort chains, which take place at the level of the mitochondrial inner membrane and the thylakoid, respectively, through chemi-osmotic mechanisms coupling electron transpn>ort and <span class="Chemical">ADP phosphorylation [1]. An ubiquinol-O2 terminal oxidase (alternative oxidase, AOX) which competes with complex III for electrons is also found in the mitochondrial inner membrane and is responsible for a cyanide-insensitive “alternative” respiratory pathway, opposed to the “cytochromial” pathway due to complexes III and IV. AOX activity does not contribute to the building of the electrochemical proton gradient and is therefore qualified as “energy-dissipating”. This enzyme has long been known to be responsible for heat production in the thermogenic tissues of higher plants (the spadix of Araceae) but is also thought to play important roles in non-thermogenic cells in some circumstances by limiting the production of superoxide anion by complexes I and III and accelerating the turnover of reduced cofactors to ensure a continuous operation of the primary metabolism [2].

Over the last century, efficient techniques have been developed to study respiration and photosynthesis in vivo[3]. In green micro<span class="Species">algae, they can easily be characterized in terms of O2 consumption/production rate of cell suspn>ensions using an aqueous phase Clark-type polarographic electrode. Determination of the apparent maximal activities (MA) of the cytochromial and alternative pathways is enabled by the use of specific inhibitors, i.e. <span class="Chemical">cyanide and substituted hydroxamic acids, respectively [4],[5]. Monitoring chlorophyll fluorescence is also a particularly suitable method for studying the functional properties of the photosynthetic apparatus [6]. In this field, pulse-amplitude modulated (PAM) fluorimetry is the tool of choice, since it enables to monitor chlorophyll fluorescence without any interference of the actinic light applied to induce the biological response [7]. This technology gives access to several important parameters characterizing photosynthesis, notably the quantum yield of photosystem II (ΦPSII) and the non-photochemical quenching of chlorophyll fluorescence (NPQ), which is actually made of 3 components: qE (ΔpH-dependent chlorophyll de-excitation mediated by the xanthophyll cycle), qT (transition of light-harvesting complexes from state 1 to 2) and qI (photoinhibition) [8]–[10]. In contrast to higher plants, state transitions have been demonstrated to be very dynamic in green microalgae, so that qT can importantly contribute to the overall NPQ together with qE [11].

The unicellular green alga <span class="Species">Chlamydomonas reinhardtii is considered as a model to study the metabolism and bioenergetics of photosynthetic organisms [12]. As shown in Figure 1, <span class="Species">C. reinhardtii is not only able to grow photoautotrophically by using light energy to fix CO2 into organic molecules, but can also assimilate acetate as an exogenous organic carbon source under the form of acetyl-CoA through ATP-dependent enzymatic reactions [13]. These features enable cells to grow mixotrophically in the light by harnessing inorganic (CO2) and organic (acetate) carbon sources, and even heterotrophically in the dark by oxidizing acetyl-CoA through the glyoxylate and tricarboxylic acid (TCA) cycles to promote the production of reduced cofactors and ATP. The glyoxylate cycle, which bypasses the 2 decarboxylation steps of the TCA cycle, also accounts for the net biomass accumulation because its C4 intermediates can be used in biosynthetic pathways [14]. An important feature of CO2 fixation in C. reinhardtii relies on the carbon concentrating mechanism (CCM), a whole-cell enzymatic machinery enabling to increase CO2 availability in the local environment of Rubisco through the dehydration of accumulated bicarbonate to counterbalance the weak catalytic activity of the enzyme and limit its oxygenase activity under low CO2 conditions [14],[15]. The CCM consists of several isoforms of carbonic anhydrases (CA) catalyzing the interconversion of CO2 and bicarbonate in different sub-cellular compartments but also of diverse inorganic carbon membrane transporters [16],[17].

Figure 1

Assimilatory pathways of light, carbon and inorganic nitrogen in. GC, glyoxylate cycle; METC, mitochondrial electron transport chain; CETC, chloroplastic electron transport chain; e−, electrons; NR, nitrate reductase; NiR, nitrite reductase; 3-PG, 3-phosphoglycerate; G-3-P, glyceraldehyde-3-phosphate; Fdred, reduced ferredoxin; Fdox, oxidized ferredoxin.

Assimilatory pathways of light, <span class="Chemical">carbon and inorganic <span class="Chemical">nitrogen in. GC, glyoxylate cycle; METC, mitochondrial electron transport chain; CETC, chloroplastic electron transport chain; e−, electrons; NR, nitrate reductase; NiR, nitrite reductase; 3-PG, 3-phosphoglycerate; G-3-P, glyceraldehyde-3-phosphate; Fdred, reduced ferredoxin; Fdox, oxidized ferredoxin.

<span class="Species">C. reinhardtii can assimilate <span class="Chemical">nitrate and ammonium as inorganic nitrogen sources, but ammonium is preferred to nitrate when they are present together in the medium, consistently with the lower energy cost of ammonium assimilation [18],[19]. Nitrate assimilation first requires its reduction into ammonium, which can then be incorporated within organic molecules under the form of glutamate through the GS/GOGAT cycle in the chloroplast (Figure 1) [20],[21]. Reduced cofactors and ATP being necessary for inorganic nitrogen assimilation mainly originate from photophosphorylation under photoautotrophic culture conditions, but also from the glyoxylate and TCA cycles coupled to oxidative phosphorylation in mixotrophically-grown cells [22],[23].

As illustrated in Figure 1 and demonstrated by recent bioinformatics- and literature-based models of the <span class="Species">C. reinhardtii metabolic network, the assimilatory pathways of light, <span class="Chemical">carbon and inorganic nitrogen are tightly interconnected through complex exchanges of metabolites, energy and reducing power which are strictly regulated in response to environmental changes to maintain cellular homeostasis [24]–[28]. These features strongly suggest that bioenergetics is likely to be influenced by the cumulative effect of different factors and even by interactions between some of them. However, to our knowledge, most reported studies have only focused on the qualitative or semi-quantitative influence of one or a few environmental factor(s) on respiration and/or photosynthesis while other factors were kept constant.

Design of experiments (DOE) coupled to multiple regression are powerful statistical tools to model the dependence of a physical, chemical or biological process to different intrinsic or extrinsic factors with a limited number of experiments [29]. They are commonly used in diverse applied research fields, in particular for the screening of culture conditions aiming to heighten the production of biomass and/or molecules of interest (metabolites, high-added value compounds, pharmaceutical recombinant proteins, etc.) by diverse organisms [30]–[33]. The methodology has successfully been applied to different micro<span class="Species">algae to optimize culture medium for heterotrophic growth, <span class="Chemical">starch and lipid production, CO2 fixation as well as metabolite extraction for bio-industrial purposes [34]–[38].

In the present work, DOE coupled to standard least squares multiple regression have been used to model the dependence of several respiratory and photosynthetic responses upon the concomitant modulation of light intensity and <span class="Chemical">acetate, <span class="Chemical">CO2, nitrate and ammonium concentrations in the culture medium of C. reinhardtii. Bioenergetic responses of interest have been defined as the dark cellular respiration (CR) and the apparent maximal activities of the cytochromial (MACYT) and alternative (MAALT) respiratory pathways, as well as the quantum yield of photosystem II in saturating light (ΦPSII800), the gross O2 evolution (P800, apparent photosynthetic rate) and the non-photochemical quenching of chlorophyll fluorescence (NPQ800) measured under a light intensity of 800 μmolphotons.m−2.s−1, which is sufficient to saturate photosynthesis but not to induce photoinhibition [39]. The main goals of this study were to determine which environmental factor(s) induce(s) most of the response variability (in other words, which factor(s) mostly account(s) for bioenergetic plasticity) and to characterize quantitatively the influence of these major explanatory factors. Such goals have been achieved through a 2 step approach consisting of a 1st-round of multiple regression aiming to detect the factor(s) of interest, which were then selected for a 2nd-round to generate predictive models enabling to simulate the mathematical profile of their influence. We discuss the results with regards to data reported in literature and we propose biological hypotheses attempting to rationalize the present observations and to provide new promising tracks for future in-depth investigations of the molecular mechanisms underlying bioenergetic plasticity in photosynthetic organisms.

Results

Design of experiments

A DOE assuming first- and second-order effects was constructed on the basis of the features summarized in Table 1 to define a limited number of combinations of values for the different environmental factors which altogether were sufficient to cover the whole design space. JMP calculated that the DOE had to contain at least 27 combinations for an unbiased subsequent modeling, and we decided to raise this number up to 42 (+50%) to ensure more confidence toward the analyses. The resulting DOE is a 2 level fractional factorial design with additional center points (i.e. combinations for which all values are equal to the center of the working range) with some extra-points typically characterizing the central composite and box-Behnken types of designs [see Additional file 1 for a complete list of all DOE items] [33]. A 3 dimensional representative example of factor dispatching within the design space is illustrated in Figure 2 for acetate, ammonium and nitrate concentrations.

Table 1

Characteristics of the factors in DOE

Factor	Type	Unit	x min / M 1	x max / M 2
Light intensity	Continuous	μmolphotons.m−2.s−1	0	200
Acetate concentration	Continuous	g.L−1	0	1
CO2 concentration	Ordinal (2 M)	%	0.035	1.5
Nitrate concentration	Continuous	mM	0	20
Ammonium concentration	Continuous	mM	0	15

M, modality; x and x, minimal and maximal values of the working range (for continuous variables); M1 and M2, 2 modalities of CO2 concentration.

Figure 2

Factor dispatching within the design space for acetate, ammonium and nitrate concentrations. Each point is an item of the DOE and its occurrence is indicated by a surrounding framed number. Black points correspond to the fractional factorial basis of the design with additional center points. Grey and white points are extra-points characterizing the central composite and box-Behnken types of design, respectively.

M, modality; x and x, minimal and maximal values of the working range (for continuous variables); M1 and M2, 2 modalities of <span class="Chemical">CO2 concentration.

Factor dispatching within the design space for <span class="Chemical">acetate, <span class="Chemical">ammonium and nitrate concentrations. Each point is an item of the DOE and its occurrence is indicated by a surrounding framed number. Black points correspond to the fractional factorial basis of the design with additional center points. Grey and white points are extra-points characterizing the central composite and box-Behnken types of design, respectively.

Pearson's correlation coefficients (r) were calculated to verify that altogether the different DOE items were uniformly covering all directions of the design space. These analyses demonstrated that the factors did not relevantly correlate to each other (|r| ≤ 0.0382, data not shown). However, while also considering second-order effects, a much more important correlation was detected in a small number of cases (Table 2). Remarkably, an absolute correlation of approximately 0.70 was systematically observed between each single factor and its combination with <span class="Chemical">CO2 concentration. This feature does not reflect a sub-optimal experimental design (which would generate a lack of information in some regions of the design space), but is rather an artificial bias attributable to the definition of <span class="Chemical">CO2 concentration as an ordinal factor. Due to this particularity, contrarily to the other second-order interactions (which all imply 2 continuous factors), the combination between CO2 concentration and another factor cannot be expressed as a continuous product term, but rather as a 2 component term with a specific modality for CO2 concentration and a specific numeric value for the other interacting factor. In such a situation, calculating the correlation between a single factor and its interaction with CO2 concentration implies to correlate this factor with itself. The observation of an absolute correlation coefficient inferior to 1 is only attributable to the definition in DOE of a different number of points for the 2 CO2 modalities for each value of the other interacting factor. As shown in Table 2, an important correlation (|r| = 0.5940) could also be noticed between squared acetate and nitrate concentrations, which indicates a lack of test points in a particular region of the design space. Such a collinearity can generate a bias in modeling (due to leverage effects of isolated points) only if these squared concentrations are implied together in the statistical analyses (as discussed below).

Table 2

Significantly correlated effects and associated absolute Pearson's correlation coefficients (| )

Factor 1	Factor 2	|r|
[Acetate]*[CO2]	[Acetate]	0.7135
Light*[CO2]	Light	0.6968
[NH4+]*[CO2]	[NH4+]	0.7092
[NO3−]*[CO2]	[NO3−]	0.7135
[Acetate]*[Acetate]	[NO3−]*[NO3−]	0.5940

The symbol "*" is used to represent second-order effects (quadratic or interactions) of individual factors.

Measurements of the responses

Bioenergetic responses were measured for the 42 combinations of factors defined in DOE and are summarized in Additional file 1. Measurements were not replicated because a global compensation of the individual experimental errors was expected in subsequent regression processes. Figure 3 presents typical oxymetric and <span class="Chemical">chlorophyll fluorescence traces obtained for one of the center points (DOE item 4) and explains the calculations of respiratory (panel a) and photosynthetic (panel b) responses. It has to be noticed that CR and MACYT could not be measured in a few cases (DOE items 14 and 28 for CR, 26 for MACYT) but that the amount of experimental data (40 and 41, respectively) remains nevertheless largely sufficient for an unbiased subsequent modeling. As emphasized by the high relative standard deviations of the data (<span class="Disease">RSD = 29% for ΦPSII800 and ≥50% for the other responses, data not shown), every response is distributed within a wide range of values, which provides a highly favorable experimental background for regression processes.

Figure 3

Oxymetric (a) and chlorophyll fluorescence (b) traces and related calculations for a DOE center point. Displayed measurements are those of item 4 in Additional file 1. The small vertical bold dotted lines define the zones which were used for rate determinations. AU, arbitrary unit. [P], protein concentration of algal suspension. SP, saturating light pulse.

Oxymetric (a) and <span class="Chemical">chlorophyll fluorescence (b) traces and related calculations for a DOE center point. Dispn>layed measurements are those of item 4 in Additional file 1. The small vertical bold dotted lines define the zones which were used for rate determinations. AU, arbitrary unit. [P], protein concentration of algal suspn>ension. SP, saturating light pulse.

Selection of the major explanatory factor(s) (1st-round of modeling)

The 1st-round of modeling was performed using the available set of experimental data to detect the factor(s) accounting for most of the response variability. First- and second-order effects were assumed, i.e. it was hypothesized that the responses could depend on each factor linearly or quadratically (for continuous factors) but also on the interaction between 2 factors. This feature rose up to 19 the number of effects which had actually to be involved in the modeling processes, making necessary to get a previous selection to avoid any possible underestimation of important contributions. For this purpose, stepwise regressions were undertaken to define the subset of effects which would provide altogether the smallest AICc (corrected Akaike information criterion) in subsequent modeling. The selected effects are listed in Table 3 with the corresponding AICc, clearly showing that the latter is always dramatically lower than the AICc including all 19 initial effects.

Table 3

Summary of the 1st-round of modeling

CR				MA CYT				MA ALT
AICcALL = 263.56 AICcMOD = 203.75				AICcALL = 285.34 AICcMOD = 228.73				AICcALL = 264.77 AICcMOD = 219.39
R2 = 0.83 R2 adjusted = 0.81				R2 = 0.70 R2 adjusted = 0.64				R2 = 0.74 R2 adjusted = 0.68
RMSEF = 2.7 (37% of the average scale)				RMSEF = 3.4 (45% of the average scale)				RMSEF = 2.7 (40% of the average scale)
Whole-model ANOVA: p < 0.0001*				Whole-model ANOVA: p < 0.0001*				Whole-model ANOVA: p < 0.0001*
Lack-of-fit: p = 0.4295				Lack-of-fit: p = 0.9282				Lack-of-fit: p = 0.6110
	Factor	β	p-value		Factor	β	p-value		Factor	β	p-value
1	[Acetate]	0.660	<0.0001*	1	[Acetate]	0.608	<0.0001*	1	[NH 4 + ]	−0.465	<0.0001*
2	[NH 4 + ]	−0.471	<0.0001*	2	Light	0.405	0.0002*	2	Light	0.391	0.0001*
3	Light	0.367	<0.0001*	3	[NH 4 + ]	−0.355	0.0007*	3	[Acetate]	0.348	0.0004*
4	[Acetate]*[NH 4 + ]	−0.292	0.0002*	4	[NH4+]*[NO3−]	0.155	0.1125	4	[Acetate]*[Acetate]	−0.292	0.0026*
5	[Acetate]*Light	0.160	0.0297*	5	[Acetate]*Light	0.136	0.1635	5	[Acetate]*[NH 4 + ]	−0.252	0.0078*
				6	[NO3−]	−0.095	0.3243	6	[NH4+]*[NO3−]	0.248	0.0090*
								7	[NH4+]*Light	−0.171	0.0637
								8	[NO3−]	0.100	0.2671
ΦPSII 800				NPQ 800				P 800
AICcALL = −53.30 AICcMOD = −110.93				AICcALL = −25.05 AICcMOD = −78.73				AICcALL = 422.01 AICcMOD = 373.30
R2 = 0.80 R2 adjusted = 0.76				R2 = 0.52 R2 adjusted = 0.44				R2 = 0.75 R2 adjusted = 0.68
RMSEF = 0.056 (28% of the average scale)				RMSEF = 0.081 (40% of the average scale)				RMSEF = 16.4 (37% of the average scale)
Whole-model ANOVA: p < 0.0001*				Whole-model ANOVA: p < 0.0001*				Whole-model ANOVA: p < 0.0001*
Lack-of-fit: p = 0.1462				Lack-of-fit: p = 0.2226				Lack-of-fit: p = 0.0874
	Factor	β	p-value		Factor	β	p-value		Factor	β	p-value
1	Light	0.804	<0.0001*	1	[NO 3 − ]*[NO 3 − ]	−0.418	0.0011*	1	Light	0.579	<0.0001*
2	Light*[NH4+]	−0.243	0.0029*	2	[Acetate]*Light	−0.360	0.0041*	2	[NO 3 − ]*[NO 3 − ]	−0.755	0.0002*
3	[Acetate]*[Acetate]	−0.225	0.0062*	3	[Acetate]	−0.354	0.0046*	3	[Acetate]*[Acetate]	0.554	0.0039*
4	[Acetate]	0.161	0.0424*	4	[NH4+]	−0.202	0.0928	4	[NO 3 − ]	0.273	0.0043*
5	[CO2]	−0.153	0.0534	5	Light	−0.126	0.2919	5	[NH4+]*[NO3−]	−0.251	0.0083*
6	[NH4+]	0.044	0.5675	6	[NO3−]	−0.114	0.3388	6	[Acetate]*Light	−0.216	0.0214*
								7	[CO2]	0.152	0.0991
								8	[NH4+]	−0.129	0.1555
								9	[Acetate]	−0.048	0.5948

AICcALL and AICcMOD correspond to hypothetical models which would comprise every 19 initial effects and the 1st-round models which only contain the stepwise-selected effects, respectively. RMSEF standardized in terms of percentage of the average scale (i.e. the difference between the mean and minimal experimental response values) is displayed to facilitate error comparison among models. Numbers ranging from 1 to 9 classify the different effects by increasing order of individual p-value, and β-weights are also provided. Statistically significant p-values (p ≤ 0.05) are surrounded by *. Effects which are highlighted in bold were considered as major explanatory ones and selected for the 2nd-round of modeling.

AICcALL and AICcMOD correspond to hypothetical models which would comprise every 19 initial effects and the 1st-round models which only contain the stepwise-selected effects, respectively. RMSEF standardized in terms of percentage of the average scale (i.e. the difference between the mean and minimal experimental response values) is displayed to facilitate error comparison among models. Numbers ranging from 1 <span class="Species">to 9 classify the different effects by increasing order of individual p-value, and β-weights are also provided. Statistically significant p-values (p ≤ 0.05) are surrounded by *. Effects which are highlighted in bold were considered as major explanatory ones and selected for the 2nd-round of modeling.

Modeling was performed through standard least squares multiple regression with the stepwise-selected effects to establish a predictive mathematical equation [provided in Additional file 2] associating a theoretical response with each of the 42 experimental values [exhaustive list in Additional file 1]. In equations, <span class="Chemical">CO2 concentration (which is a factor for ΦPSII800 and <span class="Chemical">P800) is present under the form of an additional extension term due to its definition as an ordinal factor: this extension is equal to zero for 0.035% CO2 and different from zero for 1.5% CO2 [Additional file 2]. As summarized in Table 3, R2, R2 adjusted (which introduces a trade-off for the number of effects to enable proper comparisons among models) and fitting root-mean-square error (RMSEF) were calculated to evaluate the goodness of fit, and whole-model and lack-of-fit ANOVA tests were performed to assess the statistical significance of the models. RMSEF was also standardized in terms of percentage of the average scale of the response (i.e. the difference between the mean and the minimal experimental values) in order to facilitate comparisons of regression error among models. Remarkably, the statistical insignificance of lack-of-fit tests (p > 0.05 in every model) suggests that no important effect remained unselected through the stepwise regression processes. The values of R2, R2 adjusted (≥0.70 and 0.64, respectively, except in case of NPQ800 for which R2 = 0.52 and R2 adjusted = 0.44) and RMSEF (≤45% of the average scale) as well as the p-values of the whole-model ANOVA tests (p < 0.0001 in every model) tend to indicate that altogether the effects selected through stepwise regression account for an important part of the response variability. The relative contribution of the different effects to the models was further investigated by calculating ANOVA tests and β-weights (i.e. regression coefficients which would result from modeling with previous standardization of all variables to a mean of 0 and a variance of 1) for individual effects (Table 3). The absolute value of the latter parameter participates (together with the p-value of the ANOVA test) to characterize the extent to which the effect contributes to the model, and its sign assesses whether the effect exerts a positive or negative influence on the response. As shown in Table 3 (in which the effects are classified by increasing order of p-value), high absolute β-weights are always associated with small p-values. In most cases, β-weights and individual ANOVA tests lead to establish the same order of importance for the different effects, with the exception of P800 for which light intensity exhibits a smaller p-value than squared nitrate concentration (0.0001 versus 0.0002) but not a higher absolute β-weight (0.579 versus −0.755); however, the 2 parameters similarly point out that these effects are major explanatory ones.

In case of <span class="Chemical">P800, 2 second-order effects which had been shown to correlate to each other with |r| = 0.5940, i.e. squared <span class="Chemical">acetate and nitrate concentrations (Table 2), were selected through the stepwise regression process. These effects exhibit high absolute β-weights and small p-values (β = 0.554/p = 0.0039 and β = −0.755/p = 0.0002, respectively) comparatively to most of the other effects of the model (Table 3), suggesting that the regression could have been affected by a lack of experimental data in some regions of the design space. However, the observation of an insignificant lack-of-fit ANOVA test (p = 0.0874) and of similar R2 adjusted and standardized RMSEF values (R2 adjusted = 0.68; RMSEF = 37% of the average scale) than those of the other models (except NPQ800 which exhibits a lower R2 adjusted; see Table 3) tend to indicate that no bias was introduced due to this correlation.

As shown in Table 3, one or several effect(s) which do not exert any statistically significant influence (p ≤ 0.05 cutoff) can be found within the models, and considerable differences can even be observed among the relative contribution of the statistically significant effects as indicated by the important heterogeneity of their p-values (<0.0001 to 0.0424) and absolute β-weights (0.160 to 0.804). These features suggest that some effect(s) within the models do not substantially influence the responses despite their selection through the stepwise regression processes. Consequently, to avoid overfitting and select the effects which altogether are sufficient to explain most of the response variability, a trial-and-error method consisting of several steps of multiple regression was used: 1st-round effects were removed successively by descending order of p-value (from 9 to 1, see Table 3) and modeling was tested with the remaining effects at each step of the process. The remaining effects were considered as major explanatory ones when the removal of the less important of them (i.e. the one exhibiting the highest p-value) led to the observation of a significant lack-of-fit ANOVA test (p ≤ 0.05) and/or a R2 coefficient lower than 0.60 (except in case of <span class="Chemical">NPQ800 for which R2 = 0.52 in the 1st-round model) (data not shown). As expected, the selected effects (bold characters in Table 3) had been shown to exhibit small p-values (p < 0.01) and high absolute β-weights (≥0.252) in the 1st-round models comparatively to the unselected effects, which confirms their particularly important contribution.

As emphasized in Table 3 (bold characters), maximum 3 different environmental factors over 5 (with squared effects and second-order interactions in some cases) appear to be sufficient to explain an important part of the response variability. <span class="Chemical">Acetate concentration and light intensity are major explanatory factors for every response except for ΦPSII800, for which light intensity is the only one. In addition, <span class="Chemical">ammonium and nitrate concentrations seem to influence to a large extent respiratory (CR, MACYT, MAALT) and other photosynthetic (NPQ800, P800) responses, respectively, unlike CO2 concentration which does not appear to be an important effector of bioenergetic plasticity.

Mathematical simulation of the influence of the major explanatory factor(s) (2nd-round of modeling)

The 2nd-round of multiple regression was performed with the effects selected through the 1st-round (bold characters in Table 3) in order to build a simplified model for each response. A few unselected first-order effects were also included in the modeling processes for <span class="Chemical">NPQ800 (light intensity and <span class="Chemical">nitrate concentration) and P800 (acetate concentration) because of the selection of second-order effects containing these factors (Table 3). This process led to establish a simplified mathematical expression (Equations 1 to 6) associating a new predicted response with each of the 42 experimental values:

An exhaustive list of the 2nd-round predicted responses is provided in Additional file 1, which also summarizes the 1st-round predicted responses obtained prior to the restriction of the models to the major explanatory factor(s). AICc, R2, R2 adjusted, RMSEF, β-weights and ANOVA tests for whole-model, individual effects and lack-of-fit were calculated and are summarized in Figure 4 and Additional file 3. Comparatively to the 1st-round models, the observations that R2 and R2 adjusted decreased by maximum 0.15 and 0.14 units, respectively, and that AICc and RMSEF did not increase by more than 8.6 units and 8% of the response average scale, respectively, tend to confirm that the eliminated effects did not substantially account for the response variability. As expected, β-weights and p-values of the individual ANOVA tests are similar to those of the 1st-round [see Additional file 3], so that the relative contribution of the essential effects described in Table 3 seems to be globally conserved.

Figure 4

Experimental values () as a function of predicted responses () for the 2nd-round of modeling. The diagonal full straight line and the curved dotted lines indicate a perfect match of the model (y = ŷ) and the 95%-confidence intervals of the whole-model ANOVA test, respectively. The horizontal dotted straight line represents the mean of the response as experimentally observed. The value indicated in parentheses is the RMSEF standardized in terms of percentage of the average scale of the response (i.e. the difference between the mean and the minimal experimental values). R2 adj., R2 adjusted.

In Figure 4, experimental values (y) are plotted as a function of predicted responses (ŷ) calculated from Equations 1 to 6. The y = ŷ straight line (corresponding to a perfect match of the model) and the 95%-confidence intervals of the whole-model ANOVA tests are also represented for visual evaluation of <span class="Disease">goodness of fit and statistical significance, respectively (the more the confidence intervals are close to the y = ŷ straight line, the more the statistical significance of the model is important). R2, R2 adjusted and RMSEF of the different responses (quantifying the <span class="Disease">goodness of fit) are shown to be ≥0.60, ≥0.54 and ≤47% of the response average scale, respectively, except in case of NPQ800 for which R2 = 0.48 and R2 adjusted = 0.41. ANOVA test p-values (characterizing the statistical significance of the models) appear to be inferior or equal to 0.0002 [see Additional file 3]. In Figure 5, contour plots representing the evolution of the experimental (panel a) and predicted (panel b) responses as a function of the 2 factors exhibiting the smallest individual effect p-values [see Additional file 3] are compared and reveal closely related global profiles. For ΦPSII800, contour plots were replaced by line plots since this response was modeled with only 1 factor (light intensity) in the 2nd-round. These different observations tend to confirm that altogether the effects included in the 2nd-round of modeling are sufficient to explain an important part of the response variability.

Figure 5

Comparative profiles of experimental (a) and predicted (b) responses for the 2nd-round of modeling. Graphs were drawn using JMP on the basis of the 2 explanative factors with the highest individual ANOVA p-value. Data points are positioned on the graphs and their occurrence is indicated by surrounding numbers. For ΦPSII800, line plots are presented instead of contour plots because light intensity is the only factor of the model.

Figure 6 presents a graphical simulation of the influence of each factor as predicted by the 2nd-round models with the corresponding β-weight(s) and individual effect ANOVA test p-value(s) [see also Additional file 3]. The mathematical profile of each factor was established by applying a specific value (i.e. the arithmetic mean of all DOE items) to the other factors of the model within Equations 1 to 6 to generate a new single-unknown predictive equation. In case of quadratic profile, the factor value for which the response is maximal or minimal (for concave and convex shapes, respectively) is displayed and the β-weights and p-values are provided for the first- and second-order parameter estimates. The profile of light intensity is not presented for <span class="Chemical">NPQ800 since this factor does not relevantly account for response variability per se (β = −0.120/p = 0.3280); its presence in the model is exclusively due to a second-order interaction with <span class="Chemical">acetate concentration as detailed below. As expected, in case of linear dependence, increasing and decreasing profiles are associated with positive and negative β-weights, respectively; for quadratic relationships, the type of dependence (convex or concave) is pointed out by the sign of the β-weight of the second-order effect (characterizing the second degree coefficient of the polynomial simulation equation), which is positive for convex profiles and negative for concave ones.

Figure 6

Simulation of the influence of each major explanatory factor as predicted by the 2nd-round models. Panels a to c are for respiratory responses (CR, MACYT and MAALT, respectively). Panels d to f are for photosynthetic responses (ΦPSII800, NPQ800 and P800, respectively). The black vertical dotted lines indicate the arithmetic DOE mean, and the small black vertical full straight lines the factor value for which the response is maximal or minimal in case of quadratic profile. The grey curved dotted lines are the 95%-confidence intervals of the simulation. β-weights and individual ANOVA p-values [see also Additional file 3] are displayed and highlighted in bold characters while statistically significant (p ≤ 0.05). For quadratic profiles, these parameters are provided for the second- and first-order parameter estimates (symbolized by x and x, respectively).

As illustrated in Figure 6, CR linearly depends on <span class="Chemical">acetate concentration (β = 0.652), light intensity (β = 0.374) and <span class="Chemical">ammonium concentration (β = −0.466) with p < 0.0001 (Figure 6a). MACYT exhibits a similar dependence upon these factors with the following parameters: β = 0.602/p < 0.0001 for acetate concentration, β = 0.419/p = 0.0001 for light intensity and β = −0.337/p = 0.0015 for ammonium concentration (Figure 6b). MAALT is linearly modulated by light intensity (β = 0.395/p = 0.0004) and ammonium concentration (β = −0.450/p < 0.0001), and is also quadratically influenced by acetate concentration (β = −0.296/p = 0.0058 and β = 0.344/p = 0.0016 for the second- and first-order parameter estimates, respectively; maximum for 0.623 g.L−1) (Figure 6c). ΦPSII800 and P800 linearly increase with light intensity (β = 0.810 and 0.572, respectively) with p < 0.0001. In addition P800 also quadratically depends on nitrate (β = −0.875/p = 0.0002 and β = 0.272/p = 0.0145 for the second- and first-order parameter estimates, respectively; maximum for 10.5 mM) and acetate (β = 0.657/p = 0.0034 and β = −0.049/p = 0.6448 for the second- and first-order parameter estimates, respectively; minimum for 0.497 g.L−1) (Figure 6d and f). NPQ800 is quadratically modulated by nitrate concentration (β = −0.414/p = 0.0016 and β = −0.120/p = 0.3260 for the second- and first-order parameter estimates, respectively; maximum for 9.1 mM), and also linearly decreases with acetate concentration (β = −0.360/p = 0.0049) (Figure 6e).

In most cases, mathematical profiles (Figure 6) are not influenced by the other factors of the model since modifying their value only induces a translation of the graphs along the ordinate axis without alteration of their general shape. However, as revealed by the occurrence of product terms within Equations 1 to 6 (second-order interactions), a mutual influence between the individual effects of acetate and ammonium concentrations can be observed in case of CR (β = −0.297/p = 0.0003) and MAALT (β = −0.263/p = 0.0127), as well as a strong influence of light intensity on the individual effect of acetate concentration in case of NPQ800 (β = −0.366/p = 0.0043). The negative sign of β-weights indicates that heightening the value of one of the interacting factors leads to a decrease of the slope (i.e. the first degree coefficient of the simulation equation) characterizing the profile of the other factor. In Figure 7, the influence of each interacting factor is independently simulated for 2 different values of the other one (i.e. the minimal and maximal values of the DOE range; see Table 1) by the same procedure as described for Figure 6 (the third, non interacting factor is kept equal to the arithmetic DOE mean). No graph is presented for light intensity in case of NPQ800 since this factor poorly contributes to the model (β = −0.120/p = 0.3280) whatever the acetate concentration is. Remarkably, for CR and MAALT, the effects of acetate and ammonium concentrations appear to be weakened and strengthened, respectively, by increasing the value of the other factor (Figure 7a and b). Concerning MAALT, rising ammonium concentration also generates a displacement of the optimal acetate concentration towards smaller values (maximum for 0.727 and 0.514 g.L−1 in case of 0 and 15 mM ammonium, respectively). For NPQ800, a relevant effect of acetate concentration is only observed upon high irradiance (Figure 7c).

Figure 7

Second-order interactions characterizing the 2nd-round models. Panels a, b and c are for CR, MAALT and NPQ800, respectively. The β-weight and individual ANOVA p-value are provided for each interaction. In panel b, the small black vertical full straight lines with surrounding numbers indicate the optimal acetate concentrations.

Second-order interactions characterizing the 2nd-round models. Panels a, b and c are for CR, <span class="Chemical">MAALT and NPQ800, respectively. The β-weight and individual ANOVA p-value are provided for each interaction. In panel b, the small black vertical full straight lines with surrounding numbers indicate the optimal acetate concentrations.

Validation of the 2nd-round models: cross-validation and experimental validation on new random combinations

In parallel to the main goals of the present study (i.e. to explain response variability and to characterize the influence of the major explanatory factors), we wondered whether the 2nd-round models could also be used as tools to predict the responses associated with any combination of factors within the range of DOE (Table 1).

In order to investigate the predictive ability of the 2nd-round models, experimental data were divided into k = 5 subsets (defined in Additional file 1) and k-fold cross-validation was undertaken. Each of the 5 test sets was compared to the corresponding training model (built from all data except those of the test set) in terms of cross-validation mean absolute error (MAECV) and root-mean-square error (RMSECV). Table 4 presents for each response the average MAECV and RMSECV summarizing the deviation of the 5 test data sets to their respective training model, as well as the average R2, R2 adjusted, MAEF (<span class="Disease">fitting mean absolute error) and RMSEF summarizing the analytical error of the 5 training models (for an exhaustive description of cross-validation results, see Additional file 4). The average scales of the different responses are also provided as a reference to assess the extent of the deviations. As shown in Table 4, the R2, R2 adjusted, MAEF and RMSEF of the training models are identical or very similar to those of the 2nd-round models (which are also presented in Table 4 to facilitate result overview). Remarkably, the MAECV and RMSECV characterizing the deviation of the test data sets to their respective training model do not exceed the analytical error of more than 37% for MAECV and 21% for RMSECV, and these proportions are reduced to 29% and 12%, respectively, if <span class="Chemical">P800 is not taken into account.

Table 4

Summary of cross-validation and experimental validation of the 2nd-round models

	2nd-round models				K-fold cross-validation (average values fork  = 5)						Experimental validation sets		Average scale
					Training models				Test sets
	R 2	R 2 adj.	MAE F	RMSE F	R 2	R 2 adj.	MAE F	RMSE F	MAE CV	RMSE CV	MAE EV	RMSE EV
CR	0.81	0.78	2.2	2.9	0.81	0.79	2.2	2.8	2.6	3.1	1.9	2.2	7.4
MA CYT	0.64	0.61	2.7	3.5	0.66	0.62	2.6	3.5	3.2	3.9	1.7	2.1	7.5
MA ALT	0.64	0.59	2.2	3.0	0.65	0.59	2.1	3.0	2.7	3.2	1.8	1.9	6.8
ΦPSII 800	0.66	0.65	0.053	0.068	0.65	0.64	0.053	0.068	0.054	0.067	0.068	0.078	0.201
NPQ 800	0.48	0.41	0.062	0.084	0.50	0.41	0.061	0.083	0.076	0.093	0.039	0.044	0.205
P 800	0.60	0.54	13.6	19.6	0.60	0.54	13.4	19.6	18.4	23.7	17.3	23.4	43.9

Cross-validation was performed by the k-fold method with k = 5; the R2, R2 adjusted, MAE and RMSE presented for the training models (“F”) and the test data sets (“CV”) are the average values compiling the 5 iterations. For each response, the average scale (i.e. the difference between the mean and minimal experimental values of the DOE study) is also provided as a reference to assess the importance of the deviations. R2 adj., R2 adjusted.

In addition to cross-validation, the predictive ability of the 2nd-round models was also evaluated by measuring the bioenergetic responses for 6 new randomly-generated combinations of factors (different from any DOE item) with 3 experimental replicates each (independent cultures and functional measurements) (Table 5). Every algal culture was undertaken with 1.5% CO2 for technical convenience since CO2 concentration was not included in the 2nd-round models. Because of the quite large standard deviations related to NPQ800 measurements (RSD = 19.7% on average), 2 more combinations were tested for photosynthetic responses (total number = 8) in order to increase confidence toward the general tendency of data. The predicted responses associated with the random combinations were calculated from Equations 1 to 6, as summarized in Additional file 4. The deviation of the experimental validation data sets (containing 6 or 8 items) to the 2nd-round models was assessed in terms of experimental validation mean absolute error (MAEEV) and root-mean-square error (RMSEEV) (Table 4); this deviation is also illustrated in Figure 8, in which the experimental values are plotted as a function of the predicted ones. As shown in Table 4, the MAEEV and RMSEEV of the experimental validation data are inferior to the analytical error of the 2nd-round models for CR, MACYT, MAALT and NPQ800; for ΦPSII800 and P800, they exceed the analytical error of only 28 and 27% for MAEEV and 15 and 19% for RMSEEV, respectively. As clearly evidenced while comparing Figure 4 and Figure 8a, it has to be noticed that the responses measured for the new random combinations do not cover the full range observed in the DOE study (except for P800). This feature could be attributable to the random choice of the factor values for the 6 or 8 validation points without consideration of the responses predicted by the 2nd-round models.

Table 5

Randomly-generated combinations of factors and associated mean experimental responses used for experimental validation

	[Ac.]	Light	[NH 4 + ]	[NO 3 − ]	[CO 2 ]	CR	MA CYT	MA ALT	ΦPSII 800	NPQ 800	P 800
1	0.58	61	10.0	7.0	1.5	8.1 ± 0.5	7.3 ± 0.2	10.1 ± 0.2	0.273 ± 0.019	0.218 ± 0.020	50.3 ± 1.7
2	0.32	58	15.0	3.0	1.5	6.7 ± 0.3	5.8 ± 0.3	6.6 ± 0.4	0.250 ± 0.025	0.183 ± 0.028	58.9 ± 2.4
3	0.21	90	15.0	9.0	1.5	6.4 ± 0.5	6.1 ± 0.2	8.0 ± 1.0	0.369 ± 0.005	0.235 ± 0.052	64.0 ± 0.4
4	0.75	65	12.0	5.0	1.5	10.6 ± 0.8	8.8 ± 0.6	8.3 ± 0.8	0.307 ± 0.064	0.166 ± 0.075	48.8 ± 1.6
5	0.35	144	0.0	10.0	1.5	9.0 ± 1.3	8.3 ± 0.8	12.9 ± 0.9	0.428 ± 0.013	0.244 ± 0.042	134.9 ± 9.1
6	0.25	75	7.5	10.0	1.5	5.2 ± 0.7	5.2 ± 0.7	7.2 ± 0.6	0.270 ± 0.022	0.273 ± 0.025	71.3 ± 3.9
7	0.82	167	0.0	0.0	1.5				0.467 ± 0.054	0.097 ± 0.021	15.2 ± 3.2
8	1.00	21	15.0	0.0	1.5				0.260 ± 0.015	0.133 ± 0.024	23.2 ± 1.9

Data result from 3 independent sets of cultures and measurements. [Ac.], acetate concentration.

Figure 8

Mean experimental values as a function of predicted responses for random combinations in Table5. Mean experimental values (ȳ) result from 3 independent sets of cultures and measurements. Predicted responses (ŷ) were calculated from Equations 1 to 6 (2nd-round models). The diagonal full straight line represents a perfect match of the model (y = ŷ). In panel a, axis scaling is identical to that of Figure 4 to facilitate visual comparison. The framed zone in panel a is enlarged in panel b, in which scaling is adapted to data points and standard deviations are drawn except in case of they are comprised within the limits of the markers. For P800, standard deviations are presented in panel a since it was not necessary to adapt scaling.

Data result from 3 independent sets of cultures and measurements. [Ac.], <span class="Chemical">acetate concentration.

Mean experimental values as a function of predicted responses for random combinations in Table5. Mean experimental values (ȳ) result from 3 independent sets of cultures and measurements. Predicted responses (ŷ) were calculated from Equations 1 to 6 (2nd-round models). The diagonal full straight line represents a perfect match of the model (y = ŷ). In panel a, axis scaling is identical to that of Figure 4 to facilitate visual comparison. The framed zone in panel a is enlarged in panel b, in which scaling is adapted to data points and standard deviations are drawn except in case of they are comprised within the limits of the markers. For <span class="Chemical">P800, standard deviations are presented in panel a since it was not necessary to adapt scaling.

Altogether, these different validation results indicate that the 2nd-round models can be used as tools to predict the responses associated with any combination of factors, with an inherent average deviation being quantified by the MAEF and RMSEF characterizing the analytical error of the models. As deduced from Table 4, this deviation is comprised between 26 and 36% of the response average scale for MAEF and 34 and 47% for RMSEF.

Discussion

In the present work, DOE coupled to standard least squares multiple regression have been used to model the dependence of different respiratory (CR, MACYT, <span class="Chemical">MAALT) and photosynthetic (<span class="Chemical">P800, ΦPSII800, NPQ800) responses upon the concomitant modulation of light, carbon and inorganic nitrogen sources in the culture medium of C. reinhardtii. This methodology was applied to characterize the extent to which the different environmental factors contribute to bioenergetic plasticity (through a 1st-round of modeling) as well as the mathematical profile of their influence for those accounting for most of the response variability (through a 2nd-round of modeling). Altogether, these analyses provide an overview of the bioenergetic adaptations resulting from global changes in culture conditions. This type of sequential statistical approach, which is commonly undertaken for the optimization of industrial production yields and the design and analysis of “-omics” experiments, had never been used to characterize the bioenergetic plasticity of photosynthetic cells. The individual influence exerted by one or a few environmental factor(s) (maintaining the others constant) on the cellular bioenergetics and metabolism had extensively been studied independently, but little information was available concerning their cumulative effect and their relative contribution to bioenergetic plasticity in a context in which they vary concomitantly in the medium.

The present analyses demonstrate that maximum 3 environmental factors over 5 are sufficient to explain most of the response variability (Table 3) and remarkably evidence squared effects and second-order interactions in some cases (Figure 6; Figure 7). As shown in Figure 4, comparatively to the other responses, lower R2 and R2 adjusted characterize the 2nd-round models obtained for NPQ800 (this difference was also noticed in the 1st-round models; see Table 3). Such discrepancies could (at least partly) be attributable to the higher experimental error inherent to NPQ800 measurements (see in Table 4 and Figure 8b the large standard deviations among the 3 experimental replicates of the random combinations tested for model validation: on average, RSD = 19.7% for NPQ800 and ≤9.0% for the other responses).

In order to check whether the 2nd-round models could also be used to predict the responses associated with any combination of factors within the range of DOE (Table 1), k-fold cross-validation and experimental validation tests on new random combinations have been undertaken. The similarity between the deviation of the validation points (quantified using MAECV/MAEEV and RMSECV/RMSEEV) and the analytical error of the training or 2nd-round models (for cross-validation and experimental validation, respectively) tend to confirm the predictive ability of the 2nd-round models (Table 4). It must nevertheless be emphasized that deviations of 26 to 36% of the response average scale (in terms of MAEF) or 34 to 47% (in terms of RMSEF) corresponding to the analytical error of the 2nd-round models are inherent to the predictions.

In literature, O2 evolution commonly appears to be normalized in terms of <span class="Chemical">chlorophyll concentration. However, in the present analyses, it was rather chosen to use protein concentration because of the high dependence of <span class="Species">C. reinhardtii pigment content upon culture conditions, particularly light and acetate [40]–[42]. Our unconventional normalization strategy could therefore generate apparent discrepancies between previously reported studies and the present results in some cases.

In the following sections, the authors will attempt to propose literature-based hypotheses addressing the possible biological implications of their observations. They insist on emphasizing that these hypotheses must not be considered as firm assertions, but rather aim to provide tracks for future in-depth molecular investigations.

Light stimulates CO2 fixation through the Calvin cycle and provides mitochondrial respiration with oxidizable substrates

The gross photosynthetic O2 evolution and the quantum yield of photosystem II of <span class="Species">C. reinhardtii cells adapted to moderate light intensities (0–200 μmolphotons.m−2.s−1) have been measured under 800 μmolphotons.m−2.s−1. Under such a so-called “saturating” intensity, the electron transpn>ort rate (ETR) is not limited by light availability but rather by the capacity of downstream metabolic pathways that consume photo-generated reductant and <span class="Chemical">ATP (such as the Calvin cycle). In these conditions, the gross O2 evolution (which is partly mediated by the rate of water photolysis) can primarily be considered as representative of the capacity of these pathways, even if numerous studies indicate that complex photosynthesis-associated O2-consuming processes (particularly PTOX chlororespiration and Mehler reaction) can also importantly contribute in some circumstances to modulate this response in C. reinhardtii[43].

The present analyses indicate that light intensity exerts a positive linear influence on ΦPSII800 and <span class="Chemical">P800 (Figure 6d and f). Accordingly, the maximal gross O2 evolution was previously reported to be doubled in <span class="Species">C. reinhardtii cells grown under illumination of 400 μmolphotons.m−2.s−1 comparatively to a lower illumination of 50 μmolphotons.m−2.s−1[40]. These observations could (at least partly) be attributable to the well-known stimulation of the expression and activity of Calvin cycle enzymes by light [44], in good agreement with the higher CO2 fixation rates observed upon increasing illumination in C. reinhardtii[45]. This improvement of CO2 fixation by light was reported to be correlated to higher cellular metabolite content, respiratory O2 consumption and TCA cycle-mediated CO2 production, in line with the linear stimulation of CR, MACYT and MAALT by light which could also be detected here (Figure 6a to c). For MACYT and MAALT, the term “apparent” is used because measurements were carried out on entire cells but not on isolated mitochondria. The availability of respiratory substrates could therefore not be directly controlled, so that the measured maximal activities could have been underestimated comparatively to the actual capacities if the intracellular reductant concentration was insufficient to saturate the mitochondrial electron transport chain in the presence of KCN or salicylhydroxamic acid (SHAM).

In apparent contradiction with these considerations, the present analyses did not retain <span class="Chemical">CO2 concentration as a major explanatory factor of bioenergetic plasticity (Table 3). Such an absence of influence had already been highlighted for the maximal gross O2 evolution in a previous study, in which the sum between the net O2 evolution monitored under 600 μmolphotons.m−2.s−1 and the dark respiration measured before illumination was shown to be similar in low and high <span class="Chemical">CO2-grown C. reinhardtii cells [46]. These observations could be explained by the existence of a low CO2-inducible CCM in C. reinhardtii, by which a high CO2 availability for Rubisco is maintained in low CO2 condition. Several transcriptomic analyses demonstrated that adaptation to different CO2 concentrations mainly occurs through the regulation of the genetic expression of CCM components but not Calvin cycle enzymes in C. reinhardtii[47]–[49]. Moreover, transferring C. reinhardtii cells from high to low CO2 external concentration was shown to result in a transient decrease of the amount of the small and large Rubisco subunits before returning (within the time period required to induce CCM) to the levels characterizing high CO2-grown cells [50]. Altogether, these different observations and the present ones tend to indicate that adaptation to low CO2 environment in C. reinhardtii principally occurs through CCM induction but not Calvin cycle regulation.

Acetate down-regulates the capacity of the Calvin cycle and promotes its own uptake and storage to counteract the osmotic stress associated with high extracellular acetate concentrations

The present analyses indicate that <span class="Chemical">P800 depends on <span class="Chemical">acetate concentration following a quadratic convex profile with a minimal value for 0.497 g.L−1 (Figure 6f). In a previous study, the net O2 evolution measured under 600 μmolphotons.m−2.s−1 was shown to decrease with acetate concentrations ranging from 0 to 1.75 g.L−1, but rates had been normalized in terms of chlorophyll concentration and cultures conducted under the same saturating light intensity than that of measurements (600 μmolphotons.m−2.s−1) [51].

In <span class="Species">C. reinhardtii and <span class="Species">Chlorogonium elongatum (a closely related unicellular green alga), acetate is known to repress the expression of the genes encoding the small and large Rubisco subunits (rbcS and rbcL, respectively), thereby lowering the capacity for CO2 fixation through the Calvin cycle [51]–[53]. In this context, carbon originating from acetate can substitute for up to half the photoautotrophically-generated biomass content [51],[54], and light-driven photosynthetic reactions importantly contribute to provide reductant and ATP for biosynthetic acetate assimilation (as shown in Chlamydomonas mundane) [55]. In heterotrophically-grown C. reinhardtii cells, acetate storage as starch is also known to be promoted through the improvement of the expression and activity of enzymes of the glyoxylate cycle (as isocitrate lyase, ICL) and gluconeogenesis [56],[57]. In parallel to its influence on carbon metabolism, acetate inhibits C. reinhardtii heterotrophic growth beyond 0.4 g.L−1 in the medium (“substrate inhibition”) [58]. From this concentration, the osmotic potential reaches a critical value beyond which active transport processes are impaired and energy requirements for cellular maintenance are considerably heightened. Interestingly, for P800, the present analyses point out a “concentration of inflexion” (0.5 g.L−1 approximately) which is very close to the critical substrate inhibition concentration of 0.4 g.L−1 (Figure 6f). This observation tends to indicate that P800 could be influenced by 2 independent acetate-responsive metabolic processes consuming photo-generated reductant and ATP: the Calvin cycle (repressed while increasing acetate concentration due to Rubisco down-regulation) and the biosynthetic assimilation of acetate (stimulated while increasing acetate concentration, especially beyond 0.5 g.L−1, to promote acetate uptake and storage in order to attenuate the osmotic stress).

Acetate stimulates mitochondrial respiration by heightening the intracellular reductant content and the capacity of the cytochrome pathway

The present analyses indicate that <span class="Chemical">acetate concentration exerts a positive linear influence on CR and MACYT (Figure 6a and b). Accordingly, when grown in an <span class="Chemical">acetate-containing medium, C. reinhardtii cells were previously shown to exhibit a twice-enhanced respiratory rate (partly due to the improvement of the intracellular reductant content in mixotrophic condition) [59] as well as increased transcript levels for diverse components of oxidative phosphorylation, suggesting a higher capacity of the cytochrome pathway [60]. In parallel, MAALT depends on acetate concentration following a quadratic concave profile with an optimum for 0.623 g.L−1. As illustrated in Figure 6c, this response can nevertheless be considered as linearly stimulated up to 0.5 g.L−1 acetate without further increase beyond this concentration. This observation tends to indicate that substrate-saturation of the alternative pathway in the presence of KCN could be reached beyond 0.5 g.L−1 acetate, which would imply that AOX capacity is not responsive to acetate concentration. Accordingly, enhancement of the capacity of the cytochrome pathway was already suggested to contribute to the acetate-induced improvement of dark respiration without concomitant modification of AOX capacity [61].

Acetate inhibits NPQ through repression of the LHCSR3-dependent qE component

The present analyses demonstrate that <span class="Chemical">acetate concentration exerts a negative linear influence on NPQ800 (Figure 6e), but only in case of high light intensity (Figure 7c). Interestingly, the extent of qE has recently emerged as being dramatically lowered by the presence of acetate in the growth medium, as notably evidenced by Finazzi and co-workers who demonstrated that qT is the major contributor to the global NPQ in mixotrophically-grown C. reinhardtii cells [62]. Even if the molecular mechanisms underlying the functional relationship between NPQ and acetate are not yet understood, the present results tend to confirm these findings and indicate that the magnitude of the inhibitory effect of acetate on qE could depend on its external concentration. Recently, qE has been proposed to be mediated by LHCSR3, a light-harvesting complex orthologue which is only expressed upon high irradiance [63],[64]. NPQ plasticity induced in response to changing environmental conditions (such as different acetate concentrations) could therefore be disabled in the dark and low light intensities due to the down-regulation of LHCSR3. These interpretations must be considered with caution due to the impossibility to distinguish the contributions of qE and qT to the global NPQ here.

Mitochondrial respiration contributes to provide nitrate assimilation with reductant through the acetate-dependent activity of AOX

<span class="Chemical">Ammonium concentration is shown here to exert a negative linear influence on CR, MACYT and <span class="Chemical">MAALT (Figure 6a to c). In Selenastrum minutum (another green alga), mitochondrial respiration was proposed to play a role in nitrate assimilation by acting as a trigger factor for the TCA cycle. This would in turn promote the production and export of reductant in the cytoplasm and the chloroplast to support nitrate reduction [65]. In C. reinhardtii, the enzymatic activity and genetic expression of proteins involved in nitrate assimilation are known to be repressed by ammonium [20],[21]. These regulatory events are responsible for a strict control of inorganic nitrogen uptake and assimilation by ammonium availability and enable to preferentially exploit this reduced N form if nitrate is also present in the medium [18]. Such a primary control of nitrate assimilation by ammonium could rationalize the present observations with regards to the postulated role of mitochondrial respiration in this metabolic process.

For <span class="Chemical">MAALT, <span class="Chemical">ammonium concentration is the factor which explains the highest proportion of response variability (β = −0.450/p < 0.0001; Figure 6c). Interestingly, the gene encoding AOX (Aox1) is known to be located within a gene cluster which also encodes components of the nitrate assimilatory pathway and is tightly regulated by the nitrogen source [66]; consequently, AOX expression and capacity were shown to be induced by nitrate and repressed by ammonium in a concentration-dependent manner [61]. With regards to the peculiar genetic localization and regulation of Aox1, the postulated role of mitochondrial respiration in nitrate assimilation was proposed to be essentially mediated by AOX, as also indicated by a recent comparative proteomic study published by our group [67]. The present results are in good agreement with these findings.

Interestingly, a mutual influence could be detected between the individual effects of <span class="Chemical">acetate and <span class="Chemical">ammonium concentrations for CR and MAALT (second-order interactions). As illustrated in Figure 7a and b (left panels), the effect of acetate concentration is attenuated by ammonium. For MAALT, there is also a displacement of the optimal acetate concentration toward smaller values with increasing ammonium concentration (0.727 and 0.514 g.L−1 for 0 and 15 mM ammonium, respectively). These results are consistent with a negative influence of ammonium concentration on AOX capacity. They also tend to confirm that the involvement of mitochondrial respiration in nitrate assimilation is essentially mediated by AOX since no relevant second-order interaction was retained for MACYT. Reciprocally, as illustrated in Figure 7a and b (right panels), ammonium concentration exerts a relevant influence on CR and MAALT only upon high acetate concentration (these responses exhibit a basal ammonium-independent value in the absence of acetate). This observation tends to indicate that acetate assimilation could provide the TCA cycle with oxidizable substrates to support the involvement of AOX in nitrate assimilation.

Nitrate assimilation is retro-inhibited to prevent the deleterious effects of nitrite and ammonium intracellular accumulation

The present analyses demonstrate that <span class="Chemical">nitrate concentration exerts a quadratic concave influence on <span class="Chemical">P800 with an optimum for 10.48 mM (Figure 6f). In C. reinhardtii, photosynthesis is known to contribute to provide nitrate reduction with electrons (together with mitochondrial respiration as stated beyond) [22],[23],[67], so that the rate of nitrate assimilation can influence P800 in the same way as for the Calvin cycle. The effect of nitrate on P800 can therefore be thought to result (such as for acetate) from 2 distinct metabolic processes of which the relative importance varies with nitrate concentration: substrate stimulation of reductase activity (predominant from 0 to 10 mM) and retro-inhibition of nitrate reduction by nitrate-derived intracellular ammonium (predominant beyond 10 mM). Such a retro-inhibition could attenuate the production of nitrite and ammonium (despite the higher nitrate availability) and prevent the deleterious effects which would result from their intracellular accumulation (nitric oxide overproduction and buffering disturbance, respectively) [20],[68],[69].

Similarly to <span class="Chemical">P800, <span class="Chemical">NPQ800 also depends on nitrate concentration following a quadratic concave profile with an optimum for 9.09 mM (Figure 6e). Assuming that the NADPH-to-ATP stoechiometric ratio of nitrate assimilation is superior to the yield of photosynthesis, the reoxidation of photo-generated reductant may not be paralleled with ATP turnover. This feature could result in heightening ΔpH across the thylakoid membrane, which would in turn stimulate high energy state chlorophyll de-excitation (qE) in an extent depending on the rate of nitrate assimilation.

Conclusions

In the present work, DOE coupled to standard least squares multiple regression have been applied to model the dependence of respiration and photosynthesis upon light, <span class="Chemical">carbon and inorganic <span class="Chemical">nitrogen sources in C. reinhardtii through a 2 step approach consisting of 2 successive rounds of modeling. This methodology enabled to demonstrate that maximum 3 environmental factors over 5 account for most of the variability of the different responses (i.e. can induce a relevant bioenergetic plasticity) and also permitted to obtain a mathematical simulation of the influence of the major explanatory factors. Altogether, these results provide an overview of the adaptations of C. reinhardtii bioenergetic pathways to changing culture conditions and point out new promising tracks for future more specific investigations. In order to further characterize the molecular adaptations underlying the present functional observations, we have undertaken a DOE-based comparative proteomic analysis (using two dimensional-differential in-gel electrophoresis) coupled to the determination of pigment and lipid composition by chromatography, of which the results are currently being modeled.

Methods

Algal cultures

The cw15 mt+ wall-less strain 83 of <span class="Species">C. reinhardtii[70] was pre-cultured at 25°C in 1 L Erlenmeyer flasks under <span class="Disease">orbital agitation (120 rpm), moderate light intensity (75 μmolphotons.m−2.s−1) and ambient air in 400 mL of a classical pre-culture medium (NaNO3 20 mM, K2HPO4 5.4 mM, KH2PO4 4.6 mM, MgSO4 1.4 mM, CaCl2 450 μM, oligo-elements (H3BO3 180 μM, ZnSO4 75 μM, MnCl2 25 μM, FeSO4 18 μM, CoCl2 6.8 μM, CuSO4 6.3 μM, (NH4)6Mo7O24 890 nM), Tris–HCl 20 mM pH 7.2) (photoautotrophic growth). After 6 days, pre-cultured algae were span down by centrifuging at 1,000 g for 3 min and washed once in the targeted culture medium by successive resuspension and centrifugation to eliminate pre-culture medium. Algae were finally resuspended in 75 mL of targeted culture medium, transferred to lab-scale tubular photobioreactors (Multi-Cultivators MC 1000, Photon System Instruments) allowing an accurate control of light intensity (from 0 to 200 μmolphotons.m−2.s−1) and cultured for 48 h at 25°C prior to performing functional measurements. Algal cultures were bubbled either with ambient air or with a mixture composed of 98.5% ambient air and 1.5% CO2. Culture media invariably contained K2HPO4 5.4 mM, KH2PO4 4.6 mM, MgSO4 1.4 mM, CaCl2 450 μM, oligo-elements (see above), Tris–HCl 20 mM pH 7.2. Acetic acid (0 to 1 g.L−1), NaNO3 (0 to 20 mM) and/or NH4Cl (0 to 15 mM) could also be added to reach a definite concentration depending on the specificity of each culture medium. Given the high dependence of algal growth rate upon environmental conditions, algal cultures were inoculated at a variable initial cellular density allowing to reach a dry biomass concentration of 250 μg.mL−1 (approximately corresponding to an absorbance of 1.1 at 750 nm) after 48 h. To measure the dry biomass concentration, samples of algal cultures were washed twice with milliQ water and dried for 24 h at 75°C.

Measurements of respiratory responses: CR, MACYT, MAALT

48 h-old algal cultures were concentrated up to a 500 μg.mL−1 <span class="Disease">dry biomass to optimize resolution of functional measurements. For this purpose, <span class="Species">algae were span down by centrifuging at 1,000 g for 3 min, half of the culture medium was removed and algal pellets were resuspended in the remaining volume. Respiratory parameters were determined by oxygen concentration measurements using Clark electrode and related oxygraph device and software (Hansatech). Measurements were carried out in the dark at 25°C. Oxygen concentration was monitored for 4 min before adding SHAM 2 mM (as a specific inhibitor of the alternative pathway) or KCN 2 mM (as a specific inhibitor of the cytochrome pathway). After 4 more min, KCN 2 mM or SHAM 2 mM was further added, respectively, and oxygen concentration was monitored for 4 min before ending record (total time: 12 min). The experiment was carried out twice by reversing the addition order of inhibitors (KCN + SHAM and SHAM + KCN). Cellular respiration (CR) was considered as the oxygen consumption rate in the absence of inhibitors (mean of the 2 independent measurements). Apparent maximal activities of cytochromial (MACYT) and alternative (MAALT) pathways were assessed as the oxygen consumption rates in the presence of SHAM or KCN alone, respectively. In each case, the weak oxygen consumption rate remaining after adding the 2 inhibitors (residual respiration due to inhibitor-insensitive cellular oxidases) was substracted. Oxygen consumption rates were reported to protein concentration of algal suspension and expressed in nmolO2.min−1.mgproteins−1.

Measurements of photosynthetic responses: ΦPSII800, NPQ800, P800

48 h-old algal cultures were diluted down to 80 μg.mL−1 <span class="Disease">dry biomass (approximately corresponding to 8 μg/mL <span class="Chemical">chlorophylls) to avoid any light screen effect which would otherwise affect functional measurements. For this purpose, a culture sample was centrifuged at 10,000 g for 3 minutes to spin down algae and the supernatant, only containing culture medium, was added to an untreated culture sample following a 2:1 volume ratio to generate a 3-fold dilution of algal suspension. Photosynthetic parameters were determined using a PAM fluorimeter coupled to an oxymetric device (Clark electrode and oxygraph) (Hansatech). Measurements were carried out at 25°C in the presence of NaHCO3 10 mM to avoid any experimental bias due to CO2 limitation. Chlorophyll fluorescence and oxygen concentration were monitored for 6 min without actinic light and for 6 more min under a saturating actinic light of 800 μmolphotons.m−2.s−1 (total time: 12 min). 2 pulses of saturating light (5,000 μmolphotons.m−2.s−1 for 2.5 s) were applied after 5 and 11 minutes to close transitorily all photosystem II QA centers, resulting in a brief fluorescence raise before returning instantaneously to the baseline level. Chlorophyll fluorescence traces were analyzed to determine the value of F (fluorescence baseline under 800 μmolphotons.m−2.s−1) as well as FM and FM’ (maximum fluorescence values reached during the saturating pulse in the dark and under 800 μmolphotons.m−2.s−1, respectively). These values were used to calculate ΦPSII800 = (FM’-F)/FM’) [8] and NPQ800 = (FM-FM’)/FM’) [9]. Oxygen concentration traces were employed to calculate P800 by adding the dark oxygen consumption rate to the net oxygen release rate monitored under 800 μmolphotons.m−2.s−1. Oxygen consumption rate was reported to protein concentration of algal suspension and expressed in nmolO2.min−1.mgproteins−1.

Determination of protein concentration

In case of respiratory and photosynthetic measurements, 1 mL and 6 mL of algal suspension, respectively, were centrifuged at 10,000 g for 5 minutes. Algal pellets were invariably resuspended in 1 mL extraction buffer (<span class="Chemical">NaCl 150 mM, <span class="Chemical">EDTA 1 mM, Triton X-100 1%, Tris–HCl 50 mM pH 7.5) and 25 mg polyvinylpolypyrrolidone (PVPP, insoluble in aqueous solution) were further added to complex polyphenols which would otherwise interfere with protein assay. Samples were sonicated at 3 Amp for 30 s on ice (Sonifier Cell Disruptor B-12, Branson) and thoroughly vortexed for 5 min at 4°C. PVPP was span down by centrifuging at 10,000 g for 5 min and the supernatant was used for determination of protein concentration (mg.mL−1) by a Reagent Compatible/Detergent Compatible protein assay kit (Bio-Rad) derived from the Lowry-Ciocalteu method (for details see the manufacturer’s instructions) [71].

Statistical analyses related to DOE and multiple regression

All statistical analyses were carried out using the JMP 10 software from SAS Institute. The DOE was constructed using the “custom designer” module with a single continuous response as described in Table 1. First- and second-order effects were considered, i.e. it was assumed that the response could depend on each factor linearly or quadratically (for continuous factors) but also on the interaction between 2 factors. The presence of 7 “center points”, i.e. combinations of factors for which all the values are equal to the center of the working range, was imposed in DOE.

Stepwise regressions were undertaken with minimum AICc as stopping rule and the following features as selection process: forward (progressive selection of the effects) and combine (automatic selection of a first-order effect in case of a second-order one containing the same environmental factor is selected). Comparatively to AIC (the corresponding “uncorrected” Akaike information criterion), AICc introduces a supplemental penalty while increasing the number of factors and experimental data to avoid overfitting. It is believed to exhibit theoretical and experimental advantages in the field of regression-based modeling in comparison with other likelihood-derived statistical parameters such as the Bayesian information criterion (BIC) [72].

Modeling was performed through standard least squares multiple regression to minimize the error mean square between experimental (y) and predicted (ŷ) responses by adjusting the coefficients of the following type of equation:

where b is the intercept (constant term), x the factors, b, b and b the respective linear, interaction and squared parameter estimates and e the residual.

One-way ANOVA statistical F-tests, coefficients of multiple determination (R2 and R2 adjusted) and fitting root-mean-square error (RMSEF) were calculated on the basis of total, model and error sum of squares (SS) resulting from the distances illustrated in Figure 9. ANOVA tests were undertaken on the basis of mean squares (MS) calculated from model and <span class="Disease">error SS with k and n-k-1 degrees of freedom (DF), respectively, where k is the number of factors and n the number of experimental values (MS = SS/DF). R2 was assessed as the quotient between model and total SS to estimate the proportion of the variability of the response which can be attributed to the model. R2 adjusted was calculated as <span class="Disease">1-error MS/total MS (with n-1 DF for determination of total MS) to introduce a trade-off for the number of factors, as required to perform proper comparisons of goodness of fit among different models. RMSEF was obtained from error MS and standardized in terms of percentage of the average scale (i.e. the difference between the mean and the minimal values of experimental data).

Figure 9

Distances used for SS determination in the context of R, Radjusted, RMSEand ANOVA test calculations. Distances 1, 2 and 3 were employed for the assessment of total, model and error SS, respectively.

Distances used for SS determination in the context of R, Radjusted, RMSEand ANOVA test calculations. Distances 1, 2 and 3 were employed for the assessment of total, model and <span class="Disease">error SS, respectively.

Whole-model ANOVA tests were applied to characterize the statistical significance of the models. The F-ratio was calculated as the quotient between model and <span class="Disease">error MS and the null hypothesis (i.e. MS equality) was rejected for p ≤ 0.05, indicating the variability of the response to be predominantly attributable to one or several first- or second-order effect(s).

Individual effect β-weights and ANOVA tests were calculated to compare the relative contribution of the different effects to the model. β-weights are the regression coefficients which would be obtained if the different variables had been standardized to a mean of 0 and a variance of 1 prior to modeling. The more the absolute value of this parameter is high, the more the weight of the effect within the model is important. Individual effect ANOVA tests characterize the extent to which the <span class="Disease">regression error would be increased if modeling was performed without the effect. Effect SS was assessed as the difference between the <span class="Disease">error SS of the model deprived of the effect and the actual error SS of the model, 1 DF being used for subsequent MS determination. The F-ratio was calculated as the quotient between effect and error MS and the null hypothesis was rejected for p ≤ 0.05, indicating the effect to be statistically significant (the smaller the p-value, the stronger the contribution to the model).

Lack-of-fit tests were carried out to estimate the likelihood for the model to lack one or several important effect(s). Pure <span class="Disease">error SS was assessed as the <span class="Disease">error SS of a saturated model, and lack-of-fit SS as the difference between error and pure error SS; for detailed explanation, see [73] and the “Regression Reports” webpage of the JMP online support [74]. The F-ratio was calculated as the quotient between lack-of-fit and pure error MS and the null hypothesis was rejected for p ≤ 0.05, indicating at least one important effect to be missing in the model.

Model validation

Cross-validation was performed by the k-fold method with k = 5 using the Statistica 10 software from Statsoft. The deviation of each test data set to its training model was quantified in terms of MAECV (cross-validation mean absolute error) and RMSECV (cross-validation root-mean-square error). R2, R2 adjusted, MAEF (<span class="Disease">fitting mean absolute error) and RMSEF inherent to the training models were also calculated to assess the <span class="Disease">goodness of fit. In addition, responses were measured for several randomly-generated combinations of factors (different from any combination included in DOE) and the deviation of these experimental validation data sets to the 2nd-round models was assessed in terms of MAEEV (experimental validation mean absolute error) and RMSEEV (experimental validation root-mean-square error).

Artwork and graph designs

Artwork was undertaken using the PowerPoint and Picture Manager softwares of the 2007 Microsoft Office suite. Data obtained from functional measurements were processed using the Excel software from the same company to calculate respiratory and photosynthetic responses and display them graphically. Graphs reporting results of multiple regressions and mathematical simulations (including line plots and contour plots) were drawn using the JMP 10 software from SAS Institute.

Availability of supporting data

The data sets supporting the results of this art<span class="Disease">icle are included within the art<span class="Disease">icle and its additional files.

Abbreviations

ΦPSII: Quantum yield of photosystem II

AICc: Corrected Akaike information criterion

<span class="Gene">AOX: Mitochondrial alternative oxidase

CA: <span class="Chemical">Carbonic anhydrase

CCM: <span class="Chemical">Carbon-concentrating mechanism

CR: Cellular respiration

DOE: Design of experiments

ETR: Electron transport rate

MA: Apparent maximal activity

MAECV: Cross-validation mean absolute error

MAEEV: Experimental validation mean absolute error

MAEF: <span class="Disease">Fitting mean absolute error

MS: Mean square

<span class="Chemical">NPQ: Non-photochemical quenching of <span class="Chemical">chlorophyll fluorescence

P: Gross O2 evolution

PAM: Pulse-<span class="Chemical">amplitude modulated

RMSECV: Cross-validation root-mean-square error

RMSEEV: Experimental validation root-mean-square error

RMSEF: Fitting root-mean-square error

<span class="Disease">RSD: Relative standard deviation

<span class="Chemical">SHAM: <span class="Chemical">Salicylhydroxamic acid

SS: Sum of squares

<span class="Chemical">TCA: <span class="Chemical">Tricarboxylic acid

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

GM is the author of the original idea of the work. GM and SG conceived the DOE. SG performed algal cultures and functional measurements according to FF’s expertise, as well as data analyses and modeling with GM’s contribution for statistics. SG wrote the manuscript with helpful advice from FF and GM in regards to biological and technical aspects, respectively. All authors read and approved the final manuscript.

Additional file 1

Table summarizing DOE as well as experimental and predicted responses. Predicted responses of the 1st- and 2nd-round of modeling were calculated from Additional file 2 and Equations 1 to 6, respectively. Predicted responses of the 1st-round are those obtained before the restriction of the models to the major explanatory factor(s) (this restriction was subsequently applied through the 2nd-round). Data subsets defined for k-fold cross-validation are also provided.

Click here for file

Additional file 2

Mathematical equations resulting from the 1st-round of modeling. These equations were obtained prior to the restriction of the models to the major explanatory factor(s) (as done in the 2nd-round of modeling). Due to its definition as an ordinal factor, <span class="Chemical">CO2 concentration is present as an extension term which is equal to zero for 0.035% <span class="Chemical">CO2 and different from zero for 1.5% CO2.

Click here for file

Additional file 3

β-weights, whole-model and individual effect ANOVA tests of the 2nd-round of modeling. Numbers ranging from 1 to 5 classify the different effects by increasing order of individual p-value. P-values which are surrounded by * are considered as statistically significant (p ≤ 0.05).

Click here for file

Additional file 4

Supplementary information concerning the validation of the 2nd-round models. Sheet 1 describes the parameters (R2, R2 adjusted, MAEF/MAECV, RMSEF/RMSECV and whole-model p-value) characterizing each of the 5 training models and test sets of k-fold cross-validation (set numbering corresponds to that in Additional file 1). Sheet 2 lists the randomly-generated combinations of factors used for experimental validation tests and the associated predicted responses (obtained from Equations 1 to 6).

Click here for file