Phenotypic and genetic dissection of component traits for early vigour in rice using plant growth modelling, sugar content analyses and association mapping

Highlight Plant-model-assisted phenotyping and metabolomics enable phenotypic and genetic dissection of early vigour traits in a rice diversity panel.

altitudes, and hydrologies. This broad adaptation as a species is associated with large genetic and phenotypic diversity (McNally et al., 2009). By providing rapid access to aboveand below-ground resources and the ability to compete with weeds, early vigour is a key adaptive behaviour, particularly where levels of mechanization and inputs are low (Asch et al., 1999;Namuco et al., 2009).
Early vigour is commonly estimated by the capacity of seedlings to accumulate shoot biomass rapidly. It is a complex trait, the genetic improvement of which may benefit from its dissection into component traits of lesser genetic complexity. Among such traits, Specific Leaf Area (SLA, cm 2 g -1 ), tillering, leaf size, leaf appearance, and leaf expansion rates were studied for rice and wheat (Rebetzke et al., 2007;Maydup et al., 2012). Ter Steege et al. (2005) reported negative genetic linkages among traits constituting early vigour in wheat, such as leaf size and leaf number. Tisné et al. (2008) reported, for A. thaliana, a genetic negative linkage between leaf area determination at meristem level and whole-plant leaf production.
Eco-physiological models can help to dissect phenotypic traits into genotypic and environmental parameters in a heuristic approach. The model can be as simple as a response curve (Reymond et al., 2003) or more complex (Yin et al., 1999;Nakagawa et al., 2005;Quilot et al., 2005). Associating the fitted, genotypic values of model parameters to molecular loci and polymorphisms opens new roles for plant modelling in genetics and breeding. In silico simulations can potentially predict QTL effects on the phenotype and explore optimal allele combinations for a given environment (Chenu et al., 2009;Bogard et al., 2014).
Recombinant, bi-parental populations were used for the genetic study of early vigour in cereals (Ter Steege et al., 2005;Zhang et al., 2005;Trachsel et al., 2009). Several genes were identified as controlling early vigour-related traits such as plastochron (Ahn et al., 2002), leaf size Cho et al., 2013), and tillering (Lu et al., 2013;Liang et al., 2014). Rosti et al. (2007) studied the effects of a starch deficiency mutant on vegetative growth in rice.
While, so far, plant-model-assisted QTL studies used biparental mapping populations, the approach was recently applied by Bogard et al. (2014) to a diversity panel of about 200 wheat genotypes to study the genetic control of heading date. The Genome Wide Association Study (GWAS) enables to take advantage of a high density of markers. GWAS associates SNP with the phenotype, mostly using diversity and recombination that evolved naturally. Compared with biparental QTL studies, it captures the available allelic diversity on the loci while achieving greater physical resolution due to the smaller linkage disequilibrium (LD) in diversity panels. GWAS is, therefore, an efficient way to dissect the genetic architecture of complex traits and a powerful tool for crop breeding (McNally et al., 2009;Huang et al., 2010;Courtois et al., 2013). Cereals are pioneer crops in this respect (Famoso et al., 2011;Huang et al., 2012;Matsuda et al., 2015) although applications in plant breeding are not yet common. Courtois et al. (2013) recently conducted a GWAS study on a rice tropical japonica panel for root growth and morphology traits. Regarding the diversity of above-ground traits constituting early vigour, the same diversity panel of ca. 200 tropical japonica accessions was used by Rebolledo et al. (2012a). This study confirmed the above mentioned negative correlation between leaf size and number, here represented by leaf appearance rate (syn. development rate, DR). Both traits, however, contributed positively to seedling vigour in terms of rate of shoot biomass accumulation. It was suggested that this tradeoff was due to internal competition for assimilates impacting at meristem level (Luquet et al., 2006;Dingkuhn et al., 2007), a key concept for modelling genotypic diversity in growth and phenotypic plasticity. Rebolledo et al. (2012b) described the relationship between morphological traits constituting early vigour and NSC in source and sink leaves. Highly vigorous genotypes having small but numerous leaves had low starch but high sucrose concentration in leaves. Genotypes deriving vigour from larger but fewer leaves stored more starch in leaves. Starch reserves can serve as a buffer for stress periods but apparently come at the cost of lesser vigour. Luquet et al. (2012) dissected rice early vigour into component traits using Ecomeristem, a sink driven plant growth model simulating organ appearance and growth rates and their regulation by C resource availability. Model parameters governing processes such as leaf and tiller outgrowth and leaf sizing, estimated for 200 japonica rice accessions, had significant correlations with observed morphological traits and NSC content observed on the same genotypes (Rebolledo et al., 2012a, b), suggesting this model can be useful to analyse the genetic and phenotypic diversity of early vigour traits in rice. Thus, eco-physiological models may serve as tools to analyse the genetic and phenotypic diversity of early vigour traits in rice.
The objective of the present study was to compare early vigour-related traits measured directly (morphology, NSC) or extracted by the Ecomeristem model (fitted process based parameters) and evaluate their added value in support to the phenotypic and genetic study of early vigor diversity in a rice japonica panel. The underlying morphological data, as well as NSC data for a subset of 50 accessions, have previously been reported by Rebolledo et al. (2012a, b). The heuristic approach for plant-model-assisted phenotyping with Ecomeristem was described by Luquet et al. (2012), and the genomics data obtained with GBS were described by Courtois et al. (2013). The hypothesis behind the present study, not previously tested, was that, by better taking into account trait-by-trait interactions and trade-offs among traits constituting early vigour, metabolic and plant-model-assisted phenotyping may add value to conventional phenomics in terms of GWAS results.

Plant material
The panel used in this study included 123 accessions, all belonging to the tropical japonica group except three belonging to the temperate japonica group (see Supplementary Table S1 at JXB online). All were lines or varieties adapted to rainfed upland ecosystems, except three that were adapted to irrigated flooded systems and four to rainfed-lowland systems. Roughly, half of the panel consisted of improved genetic materials and the other half of traditional varieties or landraces. Seeds of the accessions were obtained either from the Centre de Ressources Biologiques Tropicales de Montpellier or from the International Rice Research Institute (IRRI) gene bank. For each accession, the seeds were produced by single seed descent over two generations in a CIRAD Montpellier greenhouse to ensure that the samples were homogeneous.

Phenotypic data
Phenotypic data used in the present study included different types of traits (morphological, metabolic, and model parameters) that were reported in three previous publications by the same authors on 202 accessions (Luquet et al., 2012;Rebolledo et al., 2012a, b). The subset of 123 accessions selected for this study was composed in order to dispose of complete data sets for the three types of traits for each accession, and to achieve satisfactory fitting of plant-model parameters to the observed data (details in section on model-based analyses).
A greenhouse experiment with two replications was carried out from 9 February 2009 to 12 May 2009 (late winter and early spring) at CIRAD (Montpellier, France) on the panel described above. Mean daily photosynthetically active radiation (PAR) was 4.6 MJ m -2 d -1 . Air humidity and temperature were regulated by the adiabatic method and were set to 25/22 °C (day/night) and 50/90% air humidity. Seedlings were grown in a germination chamber at 29 °C. When seedlings had reached 3 cm in height, five seedlings per pot were transplanted in 1.0 l drained pots (for details, see Rebolledo et al., 2012a). The date of transplanting was variable depending on the accession and its time of germination. Pots were placed on flooded tables with a water depth of 5 cm. Plants were thinned to one plant per pot at the 4-leaf stage. Within a replication, each accession was represented by one potted plant in well-watered conditions. Pots contained about 450 g (dry weight) of a mixed soil consisting of 20% peat and 80% loamy sandy clay soil.

Growth and development measurements
For each plant, the date of appearance of leaf 6 on the main stem was noted. Tiller number (NBT), last ligulated leaf length (LLL, cm), appeared leaf number on the main stem, appeared leaf number on the whole plant (NBL), and shoot dry weight (SDW, g plant -1 ) were measured at the end of the experiment, 30-40 d after transplanting depending on the accession and its development rate. Final sampling was not performed at the same date for all accessions because of differences in phenology. Accessions had 8-9 appeared leaves on the main stem at final sampling. In order to enable comparison among accessions, observations (tiller and leaf number) were normalized by the photo-thermal time experienced by each plant from germination to final sampling, whereas LLL was normalized by its rank (for more detail, see Rebolledo et al., 2012a). Development Rate [DR (°C d -1 ), syn. phyllochron -1 ], defined as number of leaves appeared on the main stem per unit thermal time) and Relative Growth Rate (RGR, g g -1 °Cd -1 ) from germination to final sampling were computed as described by Rebolledo et al. (2012a).

Model-based analysis of phenotypes
Ecomeristem is a deterministic crop model simulating plant morphogenesis at organ level at daily time steps using weather and soil hydrological properties as input. The model version used here only simulates the vegetative growth phase of rice. It was described in detail by Luquet et al. (2006Luquet et al. ( , 2012. The plant is simulated as an average individual of a population forming a canopy. In the absence of drought, plant organogenetic and morphogenetic processes are driven by sinks (demand for current organ growth, depending on their appearance rate and potential size) and an incremental, aggregate C assimilate source (gross photosynthesis as a function of light interception and a genotypic conversion efficiency, minus maintenance respiration). The resulting daily plant C supply/ demand ratio defines Ic, a state affecting growth and development processes depending on genotypic, process-specific parameters. Genotypic parameters describe growth, development, and potential resource acquisition rates, subject to Ic and environment (temperature and PAR in the present study, as water supply was optimal). These parameters are presented in Supplementary Table S1 at JXB online. The model was calibrated for each accession using the RGenoud package of the R software (Luquet et al., 2012). For this purpose model simulations of final NBT, NBL, LLL, SDW, and leaf number on the main stem at two dates were fitted to corresponding observations under the environmental conditions encountered by each accession. In contrast to the calibration presented in Luquet et al. (2012) (model calibration on average phenotypic data of the two replicates for a given accession), parameters were estimated in the present study by fitting the model conjointly on the two replications for each accession (multi-fitting). Five of the 12 genotypic parameters described in Supplementary Table S1 at JXB online were used to decribe each accession: MGR (Meristem Growth Rate, cm); PLASTO (°Cd), SLAP (cm 2 g -1 ), ICT (unitless), EPSIB (g MJ -1 m 2 ). Only accessions for which the estimated parameters resulted in less than 15% deviation, on average, between simulated and observed variables were retained. This limited the GWAS study on model parameters to 123 accessions.
As MGR and PLASTO showed a positive correlation in the present study (r=0.32, P <0.01, with the linear regression PLASTO=44.5+1.16 MGR), similar to previous reports (Luquet et al., 2012, Rebolledo et al., 2012a on the corresponding morphological data (negative correlation between DR and LLL), the deviation from the linear regression between these two parameters was estimated to quantify the ability of a given accession to combine large leaves (LLL, MGR) and high development rate (DR, or small PLASTO) compared with the observed, average relationship regression line between PLASTO and MGR across all accessions on the panel. This deviation was called DEV_PLASTO_MGR (equation 1). The more negative DEV_PLASTO_MGR, the more the accession is able to combine large leaves and rapid development. (1) With 44.5 and 1.16, the parameters of the linear regression explaining PLASTO from MGR, both model parameters were estimated for each of the 123 accessions. Also the higher the value of DEV_ PLASTO_MGR, the lower the ability of an accession to produce large leaves at a high rate, when compared with the average linear regression.

Sugar concentration analyses
This work analysed NSC content in source leaves for the 123 accessions, thus for a larger panel than the one (50 accesions) described in Rebolledo et al (2012a).The two last ligulated leaves (youngest source leaves) were sampled in the morning (before 10.00 h) to analyse NSC content: hexoses (glucose GLU and fructose FRU), sucrose (SUC), and starch (STA). Sugar contents were analysed by enzymatic reaction with a spectrophotometric method at 340 nm using 96-well micro plates and a Janus pipetting robot as described in Gibon et al. (2009). The results are expressed as GLU equivalents per unit dry matter (mg g -1 ) for starch and as soluble sugar per unit dry matter (mg g -1 ) for hexoses and sucrose.

Statistical analyses
Statistical analyses (ANOVA, Tukey test, Multivariate Analysis) of phenotypic data were performed with statistical software R (http:// www.R-project.org) and R package FactoMineR v 1.27 as previously described in Rebolledo et al. (2012a, b). Histograms were performed with Sigma Plot software (Version 11, 2008 Systat Software, Inc). Pearson correlations and bar plots were performed using XLSTAT.

Genotyping and GWAS analysis
Genotyping was conducted at Diversity Arrays Technology Pty Ltd. (DArT P/L, Australia), using a method combining Diversity Arrays Technology (DArT Markers) and a next-generation sequencing technique based on GBS called DArTseq™ (SNPs markers). The original dataset described in Courtois et al. (2013) contains 16 664 markers (9 727 DArTs and 6 717 SNPs) genotyped in 167 japonica accessions, with a minor allele frequency (MAF) above 2.5%. To perform GWAS in the 123 tropical japonica subset, SNPs with MAF below 5% within these 123 samples were filtered leaving 12 221 SNPs (7 345 DArTs and 4 876 SNPs). No significant difference was detected between the structure observed in the 167 tropical japonica panel (Courtois et al., 2013) and the one used in this study (data not shown). The number of SNPs available for analysis, resulting in a mean distance of 22.5 kb between markers, was sufficient largely to cover the genome because linkage desequilibrium was spanning 150 kb, on average, in this population (Courtois et al., 2013). Therefore a distance of ±20 kb around each marker, corresponding to the average distance between markers, was used to identify genes of interest. SNPs linked with a strong LD (r 2 >0.7) were considered as a unique region. To find consistent and significant associations between the genetic markers and phenotypic traits, a mixed model involving population structure as fixed effect and a Kinship matrix (K) as random effect) in Tassel V4.3.6 was used. Population structure was represented by the matrix of accession scores in a Principal Component Analysis (PCA) performed on markers in Tassel using six components to make it consistent with the structure analysis reported for the same tropical japonica panel by Courtois et al. (2013). K was generated using the pairwise Identity By State distance approach implemented in Tassel. Plots representing GWAS results (Manhattan and Quantile-Quantile plots) were performed using the package QQman in R. Only associations with P-values below the thresholds of 1 × 10 -4 (named in this study as significant association) and 5 × 10 -4 (named in this study as suggestive association) were considered. The corresponding q-values for all traits were computed using the software Q-value v1.0 (Storey and Tibshirani, 2003). Q values measure the False Discovery Rate, which is the proportion of false positives among the tests found to be significant. Figure 1 shows the distribution of a wide range of values found in the diversity panel for key traits. The observed variables DR, LLL, and NBT ( Fig. 1A, B, C, respectively) were normally distributed. Among the NSC concentrations in  Table 1) measured in the japonica rice diversity panel of 123 accessions.

Phenotypic variation
source leaves (Fig. 1D, E, F), SUC showed a comparatively narrow band of relative variation (factor 2), whereas FRU concentration varied from 2 to 20 mg g -1 (factor 5) and starch from 0 to 50 mg g -1 on a dry weight basis. The distributions of genotypic model parameters (Fig. 1G, H, I) were normal with a generally smaller relative range of variation: MGR from 7 to 14 cm, PLASTO from 50 to 85 °Cd, and ICT from 0.8 to 1.8. The distribution of other traits can be found in Supplementary Fig. S1 at JXB online. To identify the relation between traits, a Multiple Factorial Analysis (MFA, Fig. 2) was performed using the variables listed in Table 1.
Supplementary Table S2 at JXB online shows the corresponding matrix of correlation. Two dimensions explained 43% of the phenotypic diversity (Fig. 2). Dimension 1 (explaining 26%) combined effects of RGR and DR (positive) and DEV_PLASTO_MGR and PLASTO (negative). RGR, the relative growth rate which was also associated with SDW, clustered with variables related to organ number (DR, NBL, and NBT). The correlations between RGR and SDW versus variables related to organ number were significant (P <0.01) (see Supplementary Table S2 at JXB online). Regarding the NSC traits, SUC was the only variable having an effect in Dimension 1 (positive). This was associated with a significant positive correlation (P <0.01) between SUC and both DR and SDW (see Supplementary Table S2 at JXB online). Dimension 2 (explaining 18%) was mainly composed of sugar-related traits (STA, NSC, GLU, FRU; all positive) and organ size (MGR as model parameter, LLL for measured trait; both positive). This is confirmed in Supplementary  Table S2 at JXB online by positive, significant correlations  between FRU and LLL, and between NSC and MGR. SDW was more explained by organ number (Dimension 1) than by organ size-related traits (Dimension 2). In general, few significant correlations were found between sugar-related traits and other types of traits.

Association mapping
Figure 3 presents Quantile-Quantile (QQ) plots and Manhattan plots for the traits with the most significant detection levels (FRU, PLASTO, DEV_PLASTO_MGR, SDW, NBL, DR). Supplementary Fig. S2 at JXB online shows QQ plots and Manhattan plots for the other traits. A total of 23 QTL were identified that were significantly associated with either morphological traits, sugar-related traits or model parameters at P <1 × 10 -4 . For approximately two-thirds of the significant associations, the P-values corresponded to q-values below 0.10 ( Table 2). Five QTL were significantly associated with morphological variables, eight with model parameters and seven with metabolic traits, the latter mainly on chromosome (chr) 11 and related to FRU and GLU in source leaves.

QTL, trait and gene associations
A total of 120 genes for the significant QTL (P <1 × 10 -4 ) were found within ±20 kb regions. Three QTL for the model parameter (DEV_PLASTO_MGR), the morphological variable (DR), and the metabolic variable (FRU_GLU) were associated with cloned genes having known phenotypic functions on panicle development, floral meristem size or meiosis (Table 2). A total of 314 genes for all the suggested QTL (P <5 × 10 -4 ) reported in Supplementary Table S3 at JXB online were found within the ±20 kb regions. Nine suggestive QTL (P <5 × 10 -4 ), related either to early vigour (SDW or RGR) or its component morphological variables (LLL), model parameters (SLAp, MGR, DEV_PLASTO_MGR) or metabolic traits (NSC), were associated with cloned genes with a known function on vegetative organs (seedling vigour, cell production, and cell size) or reproductive organs (anther dehiscence) (Table 3). Two QTL were significantly associated with SDW on chr 1 and chr 5 and a total of eight QTL were significantly associated with morphological component traits (Table 2). NBL was significantly associated with a QTL on both chr 6 and another on chr 1 ( Table 2). For NBL, a strong LD (r 2 >0.7) on chr 1 was observed between two markers (735 758 and 754 061) ( Table 2). 109 accessions carrying the T allele in position 735 758 and the A allele in position 754 061 had a lower leaf number than 14 accessions having the A allele in position 735 758 and the G allele in position 754 061 (see Supplementary  Fig. S3A, B at JXB online). However, no gene with known function was reported within this QTL region. DR was significantly associated at QTL12 on chr 7 at position 19 079 282 (Table 2) and colocated with the Bub1-Related Kinase 1 (BRK1) gene involved in meiosis (19 012 bp distance) ( Table 3) (Wang et al., 2012). For this QTL, the accessions having the A allele showed lower DR values (see Supplementary  Fig. S3F at JXB online) than accessions having the T allele. Only one QTL was significantly associated with LLL on chr 9 Fig. 3. QQ plots and Manhattan plots for traits associated to markers with the highest level of significance: FRU, PLASTO, DEV_PLASTO_MGR, SDW, NBL, DR. In the Manhattan plots the red line indicates the threshold for significant SNP association (P <1 × 10 -4 ) and the blue line indicates the threshold for suggestive SNP association (P <5 × 10 -4 ). (Table 2). LLL, DR or NBT were never associated in the same QTL region. Some of the suggestive QTL (P <5 × 10-4) (see Supplementary Table S3 at JXB online) were also associated with genes with known function (Table 3): a QTL associated with RGR on chr 2 colocated with the KRP1 gene, increasing cell size and reducing cell production (Barroco et al., 2006), and with the SGL1 gene controlling grain elongation via a SG1-Like Protein1 (Nakagawa et al., 2012).
A total of seven QTL were significantly associated with GLU and FRU only (Table 2) and were never associated together with morphological variables or model parameters. GLU and FRU were significantly associated together with QTL11 (Table 2). At this locus, the 18 accessions carrying allele A showed higher hexose concentrations (11.17 mg g -1 for FRU and 17.3 mg g -1 for GLU) compared with accessions with allele G (see Supplementary Fig. S3G at JXB online). Significant associations for GLU and FRU were also observed on chr 6 at position 9 355 224 ( Table 2). The FLO1 gene controlling a protein kinase, CLAVATA1 (Suzaki et al., 2004), involved in floral meristem size and floral organ number, was 14 050 bp distant from this marker. FRU was also significantly associated with QTL19, 14, and 22 on chr 11 (Table 2). No gene with a known function was previously reported in this region. Only one suggested association was observed for NSC on chr 5 (see Supplementary  Table S3 at JXB online); it was within the 20 kb region from the Rice S1Z1 (SUMO E3) ligase 1 protein gene that controls anther dehiscence (Thangasamy et al., 2011) (Table 3).
Among the model parameters, eight significant QTL were detected, mainly for (DEV_PLASTO_MGR) ( Table 2). DEV_PLASTO_MGR was significantly associated with four markers on chr 1, 3, 6, and 10; PLASTO was significantly associated with markers on chr 3, 5, and 10; and MGR was significantly associated with two markers on chr 6 and 9 (Table 2). Regarding PLASTO, the strongest association was observed on chr 5 (Fig. 3) at position 62 353. Accessions with the T allele showed higher PLASTO, thus lower leaf appearance rates, than accessions with the C allele at this locus (see Supplementary Fig. S3C at JXB online). PLASTO colocalized with NSC on chr 7 at position 22 559 297 but this association had P <5 × 10 -4 (see Supplementary Table S3 at JXB online). However, the phenotypic effects of these traits were located on different MFA dimensions (see Supplementary  Fig. S3D, E at JXB online). Allele A had a negative effect on NSC (105 mg g -1 versus 120 mg g -1 ) and a positive effect on PLASTO (62.79 °Cd versus 59.78 °Cd). No cloned gene was positioned near this marker.
Model parameters significantly colocalized with morphological traits for only two QTL: on chr 9 (QTL15: LLL and MGR, both related to leaf size and positively correlated; see Supplementary Table S3 at JXB online) and on chr 6 at position 9 355 224 (QTL8: RGR and DEV_PLASTO_MGR) ( Table 2), with a similar level of significance. The latter codetection for QTL8 suggests that high RGR is associated with the plant's ability to combine both high leaf production rate and large leaf size (equation 1). This QTL8 was located at 18 865 bp distance from the widely known HD1 gene controlling heading date (Taoka et al., 2011) and panicle size (Weng et al., 2014), via a CCT/B-box zinc finger protein (Table 3). In this study, 57 accessions with the T allele at QTL8 showed lower DEV_PLASTO_MGR values (0.1) and higher RGR values (1.014) compared with accessions with the C allele (see Supplementary Fig. S3D, E at JXB online). The wide distribution observed for DEV_PLASTO_MGR within each group carrying each allele explained the high standard deviation observed (data not shown).

Discussion
In spite of the considerable knowledge on traits conferring early vigour, there is little evidence on the added value of such secondary traits when explicitly included in rice breeding programmes. This may be explained by the lack of knowledge on: (i) the compensation among traits, (ii) their genotypic diversity and genetic architecture, and (iii) their behaviour in different environments. The present study aimed to combine morphological, metabolic, and model parameters for the phenotypic and genetic dissection of early vigour. It was expected from 'hidden' phenotypic markers (NSC variables and process-based model parameters heuristically estimated) to provide further insight into underlying trait-trait correlations and trade-offs. The approach was tested on a japonica rice diversity panel by performing GWAS for these traits.

Trait associations and contributions to early vigour
The three types of traits addressed in the present study were confirmed to contribute to early vigour, characterized by biomass accumulation (SDW). SDW was significantly correlated with morphological traits (LLL, NBT, NBL, DR), model parameters (MGR, PLASTO, ICT, DEV_PLASTO_MGR), and metabolic traits (SUC). This confirmed the results reported by Luquet et al. (2012) for model parameters and Rebolledo et al. (2012a, b) for metabolic traits and morphological traits in a subset of the japonica panel (50 accessions). Using MFA, the relative contribution of the various traits to SDW could be evaluated in this study. Development rate (DR, positive) or its reciprocal model parameter PLASTO (negative) contributed to early vigour. The strongest contribution Drought tolerance (Wang et al., 2011) Cloned genes in significantly associated markers (P <1 × 10 -4 ) DEV_PLASTO_ MGR Meiosis (Wang et al., 2012) to vigour was observed, however, for the derived model parameter DEV_PLASTO_MGR, expressing the combined ability of a genotype to produce both large and numerous leaves on a given culm. Plant leaf number (NBL) and tiller number (NBT) contributed similarly to vigour, whereas leaf length (LLL) and its functional model equivalent MGR alone showed a small contribution. Among the metabolic traits, leaf sucrose content showed the highest (positive) contribution to vigour whereas starch and hexose contents contributed strongly (positively) to leaf size (observed LLL and MGR model parameter) but did not have a significant effect on vigour.
It therefore appears that early vigour was favoured first, by a high development rate associated with a high sucrose concentration in source leaves, and was further increased by large leaves which were associated with high starch and hexose concentrations in the source leaves. This was translated in the model by a combination of low PLASTO and high MGR values. In the model, effects on the phenotype of PLASTO (affecting leaf appearance rate), MGR (affecting leaf size), and ICT (affecting tillering) are subject to compensatory responses due to assimilate limitation. It is interesting that PLASTO and MGR, which act at the culm scale, together had a dominant impact on early vigour, whereas ICT, which acts on tiller number at the plant scale, had a comparatively small effect. Since these parameter effects were evaluated on the basis of parameter fitting to each of the accessions, they reflect an observed pattern across the panel.

GWAS for traits contributing to early vigour
In general, traits targeted in the present study were associated with numerous markers with significant but not very high effects, confirming the evidence that early vigour is controlled by many QTL with small effect. This was also the case in other studies on other complex rice seedling traits such as root architecture (Courtois et al., 2013;Sandhu et al., 2015). The population used here for GWAS was small (123) compared with the diversity panels commonly used for such studies (>200 accessions: Huang et al., 2010;Yang et al, 2014), which probably limited the attainable levels of significance for the associations. However, this study dealt with physiological traits more difficult to phenotype than plant height or grain weight (Huang et al., 2010). Conclusive GWAS using small populations have previously been reported: betweeen 87 and 97 accessions for the components of a rice panel  or 95 for Arabidopsis (Aranzana et al., 2005).
A total of 23 QTL were identified having significant associations with P <1 × 10 -4 (Table 2). Colocalizations for different traits were rare (maximum 2). QTL for all types of traits were identified, but three model parameters (MGR, PLASTO and the derived, aggregated parameter DEV_PLASTO_MGR) and leaf fructose concentration (FRU) gave the strongest associations. This suggests that component traits derived heuristically by the Ecomeristem model or NSC measurement strengthened the analysis of early vigour. In fact, only one QTL on chr 1 (Table 2) was significantly associated with SDW, the variable defining early vigour and being interpreted here as the complex trait functionally composed of the component traits. This QTL did not colocalize with any of the other traits and could not be associated with any published, cloned gene (Table 3). It is thus an interesting marker that might be associated with a component trait of early vigour that was not measured in this study. Rebolledo et al. (2012a) suggested that DR is a key component trait explaining the genotypic variation of early vigour, through rapid deployment of organs and thus, sinks. This study indicated that genetic control of DR is related to two QTL having P <1 × 10 -4 . This result potentially opens perspectives for marker-assisted selection for DR in breeding. Two recent field studies, focusing mainly on morphological traits, evaluated biparental japonica×indica crosses (Zhou et al., 2007;Sandhu et al., 2015) and reported QTL for early vigour or seedling vigour. Several QTL identified in the present study by either morphological or model parameters were situated near genes controlling seed size (Table 3) which may then contribute to seedling vigour .

Relevance of model-based and metabolic traits compared with morphological traits
The negative correlation between leaf size and number (LLL versus DR) can be interpreted as a trade-off between two traits that contribute to early vigour, because both increase the sink strenght for a finite assimilate resource (Luquet et al., 2012;Rebolledo et al., 2012a). GWAS results did not show any shared genomic association of these traits, nor did the related model parameters MGR and PLASTO. Due to the functional similarity between MGR and LLL, they shared QTL15 on chr 9. However, the aggregated model parameter DEV_PLASTO_MGR combining leaf size and number effets was associated with three QTL that colocalized neither with MGR nor with PLASTO. It was associated, together with RGR, with QTL6. This QTL colocalized with the HD1 gene (Table 3), related to heading date, plant height, and seed weight (Lu et al., 2007). Other (non-cloned) genes were also reported in this region (Table 3). Accordingly, this region is promising to search for genes related to early vigour combining large and numerous leaves. It would also be of interest to determine if the effects of QTL6 on earliness, plant height, seed weight, and early vigour (or RGR) rely on the same gene or a cluster of complementary genes, as sometimes observed (Zhang et al., 2012). If this is the case, the trade-off between leaf size and number would not only be physiological but might also represent a physical linkage between different but complementary genes, as suported by Tisné et al. (2008) and Ter Steege et al. (2005), for A. thaliana and wheat, respectively, suggesting a strong negative genetic linkage among these traits, constraining their co-selection. Identifying QTL controlling specifically or conjointly leaf size and leaf deployment rate, as in the present study, is a crucial step to overcome such genetic linkages for improving a complex trait such as early vigour.
Metabolic traits were previously reported in rice to participate in the expression of complex traits such as early vigour (Rebolledo et al., 2012a) and grain yield (Ishimaru et al., 2007). One difficulty for the genetic dissection of this type of trait is that tissue concentrations fluctuate rapidly in response to environment (Rebolledo et al., 2012b) and plant development (He et al., 2005;Luquet et al., 2008). In this study, QTL were associated with sugar concentrations, mainly on chr 11 for hexoses (GLU, FRU) in source leaves. This region was narrow (position 7 603 153 to 7 609 931 bp) and different from that associated with total NSC concentration. On chr 6, FRU and GLU were associated with the cloned gene Floral Organ Number 1 controlling meristem size (Table 3). Hexose concentration was mainly correlated with LLL (Fig. 2), suggesting that pools of different sugars might have different physiological roles and also depend on organ type and age Rebolledo et al., 2012a). One QTL related to total NSC in the present study was associated with the PLASTO model parameter, and allelic effects were in the same direction for both traits. This is in agreement with Rebolledo et al. (2012a) who reported that high PLASTO is correlated with high leaf starch concentration. Physiologically, this was interpreted as high leaf starch concentration resulting from low sink activity for assimilates in the plant and, consequently, low vigour. However, neither starch nor sucrose concentrations were significantly associated with any genetic marker in the present study. Further studies using a larger panel and more SNP markers might elucidate QTL controlling leaf sucrose and starch concentration with greater power of trait association (Korte and Farlow, 2013).

Added value from model-assisted phenotyping
The use of modelling to support genetic studies has already proved its potential in the case of bi-parental populations, based on less-complex models and applied to less-complex traits (Yin et al., 1999;Reymond et al., 2003). One recent study combined GWAS with model-assisted phenotyping to dissect the physiological and genetic architecture of a complex trait (heading date in wheat: Bogard et al., 2014). The reasons for using heuristic applications of models in phenotyping are several (Dingkuhn et al., 2005). For example, the genetic variation of a physiological process or response trait may be difficult to measure ('hidden trait') because it is masked by environmental noise and trait×trait interactions, or it constitutes a response to environment itself (Reymond et al., 2003), requiring a model to extract the genotypic information. Once model parameters related to QTL and allele effects, the physiological model should be used subsequently to predict QTL effects on the expression of a complex trait of recombinants and in a different environment, by means of in silico extrapolation (Chenu et al., 2009). Such predictions can only be accurate if the genetic linkages and pleiotropic effects among the simulated traits are accounted for in the modelling process, and if physiological interactions with other traits are simulated (Cooper et al., 2014). With this respect, the present study provides further demonstration of the added value of modeling, and provides further validation of the physiological and genetic meaning of Ecomeristem model parameters, as previously suggested by Luquet et al. (2012). More research is needed to attain the skill to combine genetic studies and crop modelling to build in silico genetics and physiology-based ideotypes accurately. Based on direct measurements and on the heuristic estimation of model parameters, revealed new component traits of early vigour, involving phenotypic trade-offs and controlled by different QTL. Once validated through fine mapping, expression analysis, and/or development of NIL populations carrying combinations of alleles at QTL controlling leaf size and leaf development, the results might help breeders to develop plants having increased early vigour.

Conclusions
This study suggested that the complex genetic architecture of rice early vigour is controlled by several genomic regions with small but significant effects on the phenotype. It confirmed a physiological trade-off between leaf size and development rate and provided insight into their individual and common genetic control. Early vigour might be improved by combining both large leaves and high leaf deployment rate or development rate. Metabolic traits, mainly hexose content in source leaves, and Ecomeristem model parameters, in particular those controlling leaf size and leaf appearance rate, helped to understand the genetic and physiological architecture of early vigour. As a next step, Ecomeristem can be used to predict QTL effects on rice seedling performance depending on the environment. The results should be extended and validated for a wider range of genetic diversity and using a higher molecular marker density. Marker effects should be also validated in the field.

Supplementary data
Supplementary data can be found at JXB online.
Supplementary Table S1. A detailed description of the 123 accessions (name, group, system, and type).
Supplementary Table S2. The correlation matrix for all the variables studied.
Supplementary Fig. S2. shows the Manhattan and Q-Q Plot for other variables studied. Supplementary Fig. S3. The allele effects for significant QTl.