Estimation of Genetic Parameters of Biomass Production and Composition Traits in Miscanthus sinensis Using a Staggered-Start Design

Traits for biomass production and composition make Miscanthus a promising bioenergy crop for different bioconversion routes. They need to be considered in miscanthus breeding programs as they are subjected to genetic and genetic × environment factors. The objective was to estimate the genetic parameters of an M. sinensis population grown during 4 years in two French locations. In each location, the experiment was established according to a staggered-start design in order to decompose the year effect into age and climate effects. Linear mixed models were used to estimate genetic variance, genotype × age, genotype × climate interaction variances, and residual variances. Individual plant broad-sense heritability means ranged from 0.42 to 0.62 for biomass production traits and were more heritable than biomass composition traits with means ranging from 0.26 to 0.47. Heritability increased through age for most of the biomass production and composition traits. Low genetic variance along with large genotype × age and genotype × climate interaction variances tended to decrease the heritability of biomass production traits for young plant ages. Most of the production traits showed large interaction variances for age and climate in both locations, while biomass composition traits highlighted large interaction variances due to climate in Orléans. The genetic and phenotypic correlations between biomass production and composition traits were positive, while hemicelluloses were negatively correlated with all traits. Selection is difficult on young plants as the heritability is too low. The joint improvement of biomass production and composition traits would help provide a better response of miscanthus to selection.


Introduction
High-yielding crops with low-nutrient inputs and minimal negative environmental impacts are required for sustainable agriculture [1,2]. As such, species from the Miscanthus genus are perennial C4 grasses native to East Asia and identified as a promising lignocellulosic biomass crop in Europe and North America [3][4][5][6]. Therefore, this crop is a feedstock for a broad range of applications: heat and electricity production, gardening and animal bedding, as well as recent applications such as biofuel and biomaterial production [7][8][9].
Miscanthus × giganteus (M. × giganteus) [10], an interspecific hybrid between M. sacchariflorus and M. sinensis, is currently being used for most biomass production in Europe and North America with a single clone [11,12]. This clone can be grown for 20-25 years in the field, for the oldest plantations [4]. However, growing this single clone is risky, because pests and diseases could spread in non-native areas and economically threaten the overall biomass production. Moreover, the narrow genetic base of this interspecific hybrid [13] and the sterility of the cultivated clone hampers its breeding. Therefore, the development of new miscanthus cultivars is targeted with a wide adaptation to a range of different environments [14,15].
The diploid Miscanthus sinensis and the tetraploid Miscanthus sacchariflorus represent substantial genetic and undomesticated diversity, making them capable of adapting and growing in various geographic areas [16,17]. M. sinensis has a high yield potential and is adapted to drained soils, while M. sacchariflorus develops long rhizomes with an ability to spread, especially in wetland areas [16,18], which is not desirable due to the risk of it becoming invasive. For these reasons, M. sinensis is more suitable for improving biomass production in arable lands. The genetic improvement of M. sinensis will thus serve to enhance its own potential biomass production as well as M. × giganteus yield, while preserving the sterility of the new varieties in order to prevent potential invasiveness [14,19]. The breeding of M. sinensis can be indeed targeted as a high-yielding crop or as a good parent.
M. sinensis exhibits considerable intraspecific variability, with major genetic groups originating from regions across China, Korea, and Japan [20][21][22]. The breeding of new miscanthus cultivars would make it possible to extend the genetic base in miscanthus breeding programs of western countries where it is mainly obtained from ornamental accessions originating from Japan [19,21]. This would be possible by crossing distinct ornamental cultivars or exploiting wild available germplasm [14,19,23].
To target miscanthus end-uses such as heat, biofuel, biomaterial, or animal bedding, high-yielding cultivars are needed with a quality well suited for each application targeted [14]. Regarding biomass production-related traits, some key traits such as plant height, stem diameter, lateness at panicle emergence or flowering, and growth rate are positively related to biomass production [24]. The consideration of plant stem number is also of interest as the aboveground volume-calculated as the product between plant height, plant stem number, and stem section-is shown to be a good predictor of the aboveground biomass [25,26]. Plant biomass is mostly composed of the cell wall (CW). Plant cell walls contain fibrils of cellulose deeply inserted in a network of structural polysaccharides (hemicelluloses and pectins), lignins, and proteins [27]. Biomass production and biomass composition traits have to be evaluated across contrasted environments. The phenotype or expression of each trait differs according to the available genetic variability, itself determined by the genetic material studied. It also differs according to the environmental variability, determined for the years and locations considered. As miscanthus is a perennial crop, this effect of year is composed of a climatic condition effect confounded with an age effect.
Many studies have evaluated M. × giganteus, M. sinensis, or M. sacchariflorus genotypes (accessions, hybrids, etc.) for yield and morphological traits. For these traits, genotype × year interactions have been highlighted in a period covering at least 3 years after implantation. Genotype × location interactions have been evidenced as well: breeding efficiency has been shown to potentially increase if miscanthus genotypes are evaluated in different locations, as some genotypes reached a better biomass yield in one location compared to others [23,[28][29][30][31][32]. Other studies have observed genotype × year and genotype × location interactions while evaluating biomass composition traits, from at least the first to the third year after establishment [31,[33][34][35][36].
This suggests that genetic variability and environmental variability are key features which must be taken into account in the evaluation of traits for miscanthus breeding. For that reason, their characterization based on the study of a mapping population is required in order to determine the proportion of phenotypic variance that is the result of genetic factors, known as broad-sense heritability, for each trait. Broad-sense heritability estimates make it possible to predict the response to selection [37]. The higher the heritability, the more the genetic improvement will be possible through breeding. In addition to the investigation of the phenotypic correlations [38,39], it is interesting to investigate the genetic correlations between biomass composition traits, and their relationship with yield and yield-component traits. The correlations among traits provide information about the opportunity to improve a given trait along with another trait, which is potentially possible when the correlation is positive. When negative, the improvement of a trait takes place to the detriment of the other trait considered. Strong correlations can favor the improvement of both traits together, when positive, or be unfavorable when negative. Genetic correlations are most interesting to interpret, because the effect of the environment is not considered [40].
Miscanthus genetic parameters of biomass production and composition traits were studied prior to molecular-trait associations or genotypic performance evaluations [41][42][43][44][45][46][47][48]. Broad-sense heritability was evaluated for biomass production traits [43][44][45][46][47], biomass composition traits [48], or both [41,42]. Successive years and one location were considered in these studies, except Clark et al. [45] who considered several locations. Biomass production traits were found to be highly heritable (i.e., heritability value higher than 0.40) and tended to increase from the second year after establishment [41,42,[45][46][47]. Biomass composition traits were also highly heritable, but did not systematically increase over the years after establishment [41,42,48]. Only Clark et al. [45] have evaluated M. sinensis accessions in six locations and have highlighted combined genetic group × location and genotype within genetic group × location effects that represent a substantial proportion of the total variation for biomass yield (32%) in the third year after planting.
The objectives of this study were to estimate the genetic parameters, genetic and phenotypic correlations in terms of biomass production and composition traits together, during four consecutive years in a diploid biparental M. sinensis population in two contrasted locations. In a perennial crop such as miscanthus, environmental variability is itself determined by specific conditions such as location and year, and the year encompasses the age and climate effects. If we already know that the genetic parameters and correlations vary according to genetic and environmental variability, the extent to which this variation is due to specific conditions such as age, year, and location is not yet known in miscanthus as the previous studies did not decompose the year effect into age and climatic condition effects. In particular, is the selection of young plants during the yield-building phase equivalent to old plants which reached their full growth potential? Are the age effects on the genetic parameters of the same magnitude of climatic condition effects? Does the trend vary between locations? For that purpose, the population was tested according to a staggered-start design [49] that established the population twice in 2014 and 2015 in both locations. With such a design, the "year" effect can be decomposed into an "age" effect and a "climate" effect [50,51]. To date, it is the first study in miscanthus which estimates the genetic parameters and correlations when considering both biomass production and composition traits and distinguishes within the year, age, and climate effects. This will substantially broaden knowledge of the genetic evaluation of biomass traits, according to the effects of age, climate, or location.

Plant Material
The mapping population was generated by crossing two diploid ornamental M. sinensis cultivars, "Malepartus" (Mal) and "Silberspinne" (Sil), to obtain an F1 progeny. These two parents were purchased in a French nursery for a previous study [24] and their initial provenance was identified to be from central and southern Japan [21,22]. They were then evaluated in the field within a set of 21 miscanthus clones [24,31]. They were selected for their highly contrasted plant stem numbers, and then reciprocally crossed in isolated greenhouses during the fall of 2007, at the INRAE BioEco-Agro research center of Estrées-Mons in Northern France. All the seeds resulting from the cross were then stored in the refrigerator in suitable packaging at 4 °C. From 2012 to 2015, each seed was germinated in vitro and the corresponding plant was propagated in vitro according to a protocol of shoot organogenesis and regeneration [52]. This provided clonal replicates of each single genotype originating from a seed. Then, all seedlings were planted in a greenhouse that offered suitable growing conditions and grew in the greenhouse for a year before being transplanted to a field. An average of four replicates was finally obtained per genotype, with variations for some genotypes due to variable genotype ability for in vitro propagation.

Experimental Design and Planting
The experimental design of the trial was a staggered-start design [49], established in two different French locations mainly contrasted by their soil and climatic conditions and approximately 290 km apart. The first site was located at the INRAE GCIE experimental unit of Estrées-Mons (Somme, 49° 53′ N, 3° 00′ E) and characterized by a deep loam soil (Orthic Luvisol, WRB) while the second one was located at the INRAE GBFOr experimental unit in Orléans (Loiret, 47° 49′ N, 1° 54′ E) and offered a sandy soil (Dystric Cambisol, WRB). For each location, the field trial consisted in a staggered-start design where two groups of genotypes were established in two subsequent years. These different years were planted side-by-side according to two adjacent plots in the field (one plot per year). The first group of genotypes was established in 2014 and coded G1, while the second group of genotypes was established in 2015 and coded G2 (Fig. 1a). The first group G1 was observed for a 4-year period and the second group G2 for a 3-year period. In Estrées-Mons, 159 replicated genotypes, including the parents, were distributed across each plot, with most genotypes being common between the two establishment years (Fig. 1a). In Orléans, 106 replicated genotypes, including the parents, were established in each plot, with in addition a set of genotypes common to the two establishment years (Fig. 1a). In both locations, the staggered-start designs were composed of two incomplete randomized block designs, with genotypes being replicated in four of the five blocks of each on average. They were organized with single plants equally spaced within and between rows 1 m apart, which led to a plant density of 1 plant per m 2 . To ensure an equal competition effect between plants along the trial, a border row was planted on each side of the plots that contained repeated genotypes of the population. After transplanting, the field trial was watered and drip irrigation was installed, in both years and locations, thus supplying sufficient water for the overall duration of the trial. No fertilizers were applied and weed management was carried out using a hoe whenever necessary. Plants that did not survive the first year were replaced by plants of the same age that were kept in a nursery near the trial.
The strength of the staggered-start design made it possible to decompose the "year" effect into "age" and "climate" effects ( Fig. 1b). For example, when considering the "year" 2017 in one location, G1 and G2 were contrasted by  the "age" effect, because G1 plants were 3 years old while G2 plants were only 2 years old. Thus, the age effect was decomposed based on the year considered, in which a given climatic condition occurred. Plants with the same "age" of 3-year-old genotypes were then considered, for instance, with G1 observed in 2017 and G2 in 2018: G1 and G2 were now distinguished by the "climate" effect, because G1 plants were evaluated during the year 2017 and G2 plants during 2018. Therefore, each climatic condition that occurred in each of the 2 years was considered using plants of the same "age." The effect of the soil, another environmental factor, was not substantial as the two plots of the staggered-start design were adjacent in the field (in a given location). All these decompositions were applied in each location to all biomass production and composition traits in order to assess the corresponding genetic parameters, genetic and phenotypic correlations. It can be noted that the analysis of these two kinds of data required two different statistical models as described later in the statistical analysis section.

Climatic Conditions According to Plant Cycle Determined Per Year and Location
As reported previously, the staggered-start design was established in two different locations, each one characterized by specific soil and climatic conditions. Weather variables were calculated according to five periods of the plant cycle determined on Malepartus which is one parent of the population. These periods were determined for each year in which the population grew, by considering each location separately (Tables S2a and S2b). Due to the destruction of plants for the floral transition stage determination by apex dissections, the Malepartus clone was established in 2014 in a specific trial not far from the population trial, at a plant density of 2 plants/m 2 [53]. The pre-emergence phase was coded as Period 1 and lasted from senescence in the fall of year N − 1 to shoot emergence in year N (and accordingly, this period was not available for the first year). Period 2 lasted from shoot emergence to floral transition, i.e., which corresponded to the differentiation of the terminal vegetative apex into a panicle and was determined as the date of the first stem that showed the transition to floral stage. This date was determined according to Hou et al. [54]. Period 3 lasted from floral transition to panicle heading and Period 4 from panicle heading to senescence. Period 5 slightly contrasted with the other periods that were identically determined in both locations: in Estrées-Mons, it started with senescence and ended on the date of the first winter frost, while it ended on the harvest date in Orléans. This was due to the fact that winter frost always occurred after senescence in Estrées-Mons and therefore shortened the plant cycle, while no winter frost occurred in Orléans for all the years observed.
Based on the widely used BBCH scale or "Biologische Bundesanstalt, Bundessortenamt, CHemische Industrie" scale [55], it can be noted that the shoot emergence, panicle heading, and senescence dates corresponded to 09, 51, and 99 stages, respectively. Moreover, the floral transition was determined as the date of the first floral transition appearance. It corresponded to the widely used "double ridge stage" in winter wheat [56]. It was based on extensive apex observations made in Malepartus in 2018 [54] and on assessments for the other years. As for maize [57], the use of a degree-day model implying constant degree-day durations must predict the floral transition and flowering: the ratio between the growing degree-days (GDD) for the emergenceto-floral transition period and the GDD for the emergenceto-flowering period was calculated for Malepartus in 2018 and reached 0.17, this ratio was then applied for the estimations of the other years. A base temperature of 8 °C was used for the estimation of GDD, as refined for Malepartus [54]." In each location, weather variables were initially recorded daily in a meteorological station located at a short distance from the trial and downloaded from the CLIMATIK database (INRAE AGROCLIM, 2020). The weather variables were as follows: air minimum temperature (MinT), air maximum temperature (MaxT), air mean temperature (MeanT), soil minimum temperature (MinTs), precipitation (Prec_M), mean humidity (MeanH), vapor-pressure deficit (VPD), maximum vapor-pressure deficit (VPD_M), Penman potential evapotranspiration (ETPP), photosynthetically active radiation (PAR), cumulated growing degree-days (CGDD). Their mean was then summarized for each year in each location according to the five periods of the plant cycle previously defined (Tables S3a and S3b). In addition, cumulated precipitation (Prec_S) was also determined.

Phenotyping of Biomass Production and Composition Traits
Five biomass production traits were phenotyped. The following morphological traits related to biomass production were evaluated. Canopy height (CH_cm) was phenotyped in October using a measuring rod from the ground to the ligule of the last ligulated leaf. Plant maximum height (HMax_cm) was measured in October from the ground to the highest panicle of the plant (if the plant flowered), or to the horizontal flexion of the highest leaf (if the plant did not flower). The plant stem number (PSNb) was determined in October on each plant, by considering all stems with at least one ligulated leaf. Plant circumference (C50_cm) corresponded to the measurement, just before harvest, of the plant circumference 50 cm from the ground. The aboveground biomass yield (ABM_tDMha) or aboveground dry matter (DM) yield, expressed in tDM/ha, was calculated based on the fresh weight and dry weight of one plant (dry matter per ha) after the winter harvest late February. It can be noted that the traits related to flowering time were published by Hou et al. [54].
Six biomass composition traits were measured from a representative sample of approximately 500 g of fresh matter for each plant. Each sample was dried at 55 °C for 4 days in a well-ventilated oven and ground using a hammer mill grinder for it to pass through a 1-mm grid. The following traits were considered, according to the total dry matter content (%DM) or cell wall content (%CW): neutral detergent fiber, which corresponded to the cell wall content (coded NDF_%DM), acid detergent fiber (ADF_%DM), cellulose (CL_%DM; CL_%CW), hemicelluloses (HEM_%DM; HEM_%CW), acid detergent lignin (ADL_%DM; ADL_%CW), and ash (Ash_%DM).
In each location, all the traits were phenotyped over a 4-year period, from 2014 for G1 and 2015 for G2. Each trait was named according to a miscanthus ontology, developed at the INRAE BioEcoAgro research unit of Estrées-Mons, thanks to the GnpIS multispecies integrative information system [58] from the INRAE URGI of Versailles (https:// urgi. versa illes. inra. fr/ ephes is/ ephes is/ ontol ogypo rtal. do).
As too many samples were available (about 4500), a highthroughput method based on near-infrared spectroscopy had to first be developed for the measurement of biomass composition traits to reduce the cost of the analyses by using a fraction of the dry powder samples (called calibration sample). NIR spectra were measured using an Antaris II -Thermo spectrometer and expressed in absorbance with a wave number range between 4000 cm −1 and 10,000 cm −1 and a resolution of 4 cm −1 . In order to determine a representative calibration sample, the spectra of all samples of a given year were first analyzed using a principal component analysis (PCA). About 10% of the genotypes were then selected based on their position on the two first principal components of the PCA in order to best reflect the whole population for the year under consideration. Every year, the calibration sample was then biochemically analyzed in order to determine the six biomass composition traits NDF, ADF, CL, HEM, ADL, and Ash in the Laboratoire Agronomique de Normandie (LANO), using an adaptation [59] of the Van Soest method [60]. Considering all sample calibrations of all years together, a partial least square (PLS) regression was carried out on each biomass composition trait (used as a single Y quantitative response variable) while the near-infrared absorption spectra were used as an X matrix of predictors. The validation of the PLS predictions was then carried out based on independent samples. The calibration and validation parameters showed good predictions of all biomass composition traits by near-infrared absorption spectra with correlations above 0.89 between measures and predictions ( Table 1). Once the predictors were obtained by PLS, the biomass composition traits of all available samples were finally predicted.

Statistical Analysis
Phenotyping data were analyzed for each location separately. Linear mixed models [61] were carried out for the unbalanced datasets: for the age effect modeling according to each year and for the climate effect modeling according to each age. For the staggered-start design in Estrées-Mons, 159 genotypes were included, with 82 in common between the G1 and G2 groups of genotypes. As for the staggered-start design in Orléans, 106 genotypes were considered, with 59 in common between G1 and G2.
All production and composition traits were normally distributed and no data transformations were thus required prior Table 1 Characteristics of NIRS calibrations developed for composition traits in miscanthus. Traits considered were neutral detergent fiber (NDF_%DM), acid detergent fiber (ADF_%DM), cellulose (CL_%DM; CL_%CW), hemicelluloses (HEM_%DM; HEM_%CW), acid detergent lignin (ADL_%DM; ADL_%CW), and ash con-tent (Ash_%DM). Traits were expressed according to the percentage of dry matter (%DM) or the percentage of cell wall (%CW) of the plants. The standard error of prediction (SEP) is also displayed, with the number (n) of samples used for the calibration and validation steps and the associated correlations (r)  Fig. 1b). It initially consisted of a full statistical model which took into account block and spatial effects using the "breedR" package implemented in the R software [62]. It accounted for the environment heterogeneity within each trial among the blocks and diagnosed spatially distributed patterns of remaining residual variations using variograms. The remlf90 function of the package, based on restricted maximum likelihood estimations, was used to estimate all variance components of the progeny.
where Y ijkl represents the phenotypic value measured on plant k of genotype i at age j in block l; µ is the overall mean; α i is the random effect of genotype i; β j is the fixed effect of age j; (αβ) ij is the random interaction between genotype i and age j; δ l is the effect of block l and ε ijkl is the random residual for plant k of genotype i at age j in block l. An autoregressive spatial component was included in the model based on x and y coordinates in each plot in order to decompose the residual ε ijkl into a spatially dependent parameter, θ ikl , for plant k of genotype i in block l, and an independent remaining residual [63]. This full model was compared with a sub-model with no decomposition of the residual term into spatially dependent and independent effects. It was also compared with another sub-model with no block effect but with the decomposition of the residual term into spatially dependent and independent effects. Based on the Akaike information criterion (AIC) [64] of each of the three models, the last one which dropped the block effect but took into account the spatial effect was found to have a lower AIC (i.e., better performance) in all datasets (Fig. S1) and was retained afterwards. In addition, each model detailed hereafter was initially tested on balanced datasets by considering the 82 common genotypes between G1 and G2 in Estrées-Mons. Then, each model was tested on unbalanced datasets using all genotypes, which were common and non-common between G1 and G2. As models applied to balanced and unbalanced datasets yielded very similar results, all the analyses were finally carried out on unbalanced datasets, thus including more genotypes. The Model 2 finally used was written as follows and used for 2016-year, 2017-year, or 2018-year data subsets (Fig. 1b): where all the terms are described as previously without taking into account the block effect.
To account for the climate effect modeling (Model 3), the following linear mixed model was applied to each studied trait, considering each location and each age of the progeny (age 1, age 2, or age 3 data subsets in Fig. 1b): where each term is similar to Model 2 except that age effect β j is replaced by climate effect γ l of year l, and the interaction between genotype i and age j (αβ) ij by the interaction between genotype i and climate effect of year l (α'γ) il . Note that the terms α' i , ε' ikl , and θ' ik were different from the terms given by the previous Model 2 because the corresponding effects were estimated based on different datasets (Fig. 1b).
The ratio of the genetic variance to the genotype × age interaction variance and the ratio of the genetic variance to the genotype × climate interaction variance were respectively computed based on corresponding variance components estimated using Model 2 and Model 3. For each trait, these ratios made it possible to determine the part of genetic determinism mainly due to stable genetic effect or genotype × age and genotype × climate interaction effects.
According to the two previous models, individual plant broad-sense heritability sensu lato values ( H 2 sl ) were estimated using the following formulas: For the age effect modeling according to each year (Model 2), where 2 is the variance attributed to the genotype, 2 is the variance of the genotype × age interaction and 2 is the residual variance.
For the climate effect modeling according to each age (Model 3), Here, each term is similar to the previous formula, except σ 2 α'y which is the variance of the genotype × climate interaction.
According to the two previous models again, progenymean broad-sense heritability values [40,65] were also estimated for the age effect modeling (Model 2) and the climate effect modeling (Model 3), respectively, as: where J is the number of ages considered, L is the number of years with associated climates, and K is the mean number of replicates per genotype. Genetic and phenotypic correlations between traits were then assessed, considering a given year with the age effect modeling (Model 2) or a given age with the climate effect modeling (Model 3). Genetic correlations were estimated using the equation described in Howe et al. [66] which was coded "Eq. 9" in their paper. To calculate phenotypic correlations with individual values, Pearson correlation coefficients were computed using R package "stats" and visualized using the "corrplot" R package [67]. In order to characterize the climatic conditions over the years of the experiment, an initial principal component analysis (PCA) considered the weather variables determined for the five periods of the growing season for each year and location: weather variables were considered as variables and the individuals were composed of the combinations of period × year. Each principal component corresponded to a linear combination of the weather variables. Then, a subsequent Hierarchical Clustering on Principal Components (HCPC) was carried out. For this hierarchical classification, the method used the coordinates of the individuals (combinations of period × year) on the principal components, which distinguished the combinations of period × year on their climatic conditions. PCA and HCPC were performed by using the "FactoMineR" R package [68].

Biomass Composition Traits Were Generally Less Heritable than Biomass Production Traits with Substantial Variations According to Climatic Conditions and Ages for All Traits
For each trait, the means of individual plant broad-sense heritability ( H 2 sl ) and the means of progeny-mean broadsense heritability ( H 2 Pi ) were computed based on heritability values determined using each age effect for each year, or each climate effect for each age, with Model 2 or Model 3, respectively. Each location was considered separately (Table S1). This made it possible to know how much a trait was heritable, before analyzing their evolution while considering the age effect throughout the years or the climate effect throughout the ages. In Estrées-Mons, the biomass production traits were moderately to highly heritable, with means of individual plant broad-sense heritability ranging from 0.42 to 0.62 for estimates of Model 2 and from 0.43 to 0.52 for estimates of Model 3 ( Table 2). Biomass composition traits were generally less heritable, ranging from 0.29 to 0.44 for estimates of Model 2 and from 0.26 to 0.47 of Model 3. In Orléans, the ranges were similar to those in Estrées-Mons, for estimates of both models and trait types ( Table 2).
By considering the means of progeny-mean broad-sense heritability ( H 2 Pi ), which were more meaningful at clonal level, traits were highly heritable for both models and both locations, with the corresponding means ranging from 0.60 to 0.88 (Table S1). According to those means, the progenymean broad-sense heritability was around 1.5 to 2.5 higher than individual plant broad-sense heritability, most often with a higher difference for composition traits.

Lower Genetic Variances Along with Larger Interaction Variances Explained the Lower Heritability Values for Biomass Production Traits at Young Plant Ages
According to Model 2 which considered the age effect for successive years, lower genetic variances were observed for most of the production traits in 2016 than for the subsequent years in both locations (Fig. 2a). This year corresponded to younger plants that were 2 years old and 1 year old according to G1 and G2, respectively. The corresponding genotype × age interaction variance (i.e., age effect) was larger for these young plant ages (Fig. 2a) as highlighted by the lower ratio of the genetic variance to genotype × age interaction variance (Table 2). This ratio was even lower than 1 for aboveground biomass yield and lower than 2 for plant stem number and plant circumference in Orléans, respectively meaning that the extent of the genotype × age interaction variance was higher than the genetic variance or corresponded to more than 50% of it ( Fig. 2a; Table 2). Accordingly, the resulting individual plant broad-sense heritability values ( H 2 sl ) of production traits were generally lower in 2016 while they increased throughout the years. Canopy height, plant maximum height, and aboveground biomass yield continuously increased until they reached a maximal heritability level above 0.50 in 2018 for both locations ( Fig. 2a; Table S1). However, plant stem number in both locations and plant circumference in Estrées-Mons reached a maximum level one year earlier, in 2017, with values of 0.45, 0.56, and 0.49, for each case respectively ( Fig. 2a; Table S1). When comparing the individual plant broad-sense heritability values obtained in 2018 to those obtained in 2016 in Estrées-Mons, the increase varied from 10 to 30% except for the plant stem number. In Orléans, this increase ranged from 30 to 80% (Table S1). Biomass composition traits were not so highly subjected to genotype × age interaction: consequently, the heritability values were higher in 2018 for Orléans, mainly because of an increase in genetic variance ( Fig. 2b; Fig. S2a; Table 2; Table S1).
Model 3 was used to model the climate effect for successive ages: lower genetic variances were highlighted for a major part of production traits at age 1 than at older plant Table 2 For each trait in each year, Model 2 was used to determine the genetic variance ( 2 α ), the genotype × age interaction variance ( 2 αβ ), and the associated ratio was then computed and dis-  ages in both locations (Fig. 3a). For this age, the genotype × climate interaction variance (i.e., climate effect) was evaluated according to G1 plants observed in 2015 and G2 plants in 2016: the associated ratio of genetic variance to genotype × climate interaction variance was often lower than for other ages, meaning that young plants were the most subjected to climate variations (Table 2). Thus, the related individual plant broad-sense heritability values ( H 2 sl ) were generally higher at age 2 or age 3. They increased to reach a maximum level above 0.50 at age 3 for canopy height and plant maximum height in both locations ( Fig. 3a; Table S1). For the plant stem number and aboveground biomass yield, the maximal heritability was reached at age 3 in Orléans, with values of 0.56 and 0.53, respectively. In Estrées-Mons, these traits showed the highest values, 0.47 and 0.50 respectively, at age 2. When comparing the individual plant broad-sense heritability values obtained at age 3 to those obtained at age 1, the increase ranged from 20 to 100% in each location, except for the plant stem number in Estrées-Mons (Table S1).

Most of the Production Traits Showed Large Interaction Variances According to Age and Climate in Both locations While All Composition Traits Showed Large Interaction Variances According to Climate in Orléans
The ratio of the genetic variance to genotype × age interaction variance, determined according to Model 2, was relatively low for most of the production traits evaluated from Pi ) heritability values assessed for each year according to the age effect and considering locations separately. Each heritability value depended on the associated genetic variance, genotype × age interaction variance, and residual variance. The genotype × age interaction variance between G1 and G2 allowed to know the age effect for each year. For example, considering 2016, G1 plants were 2 years old and G2 plants were 1 year old. These heritability values (right y-axis) and variance components (left y-axis) were assessed for biomass production (a) and composition (b) traits. See "Materials and Methods" section for trait name. Blanks in the plots are due to nonavailable data 2016 to 2018 in both locations. For the plant stem number, plant circumference and aboveground biomass yield, this ratio ranged from 1.9 to 6.1 in Estrées-Mons and from 0.9 to 3.6 in Orléans ( Table 2). The lower the ratio values, the more genotype × age interaction variance was associated to the genetic determinism of a trait, which meant that genotype × age interaction variances were substantial for these traits and conditions ( Fig. 2a; Table 2).
The ratio of the genetic variance to genotype × climate interaction variance, determined using Model 3, was also generally low for the production traits evaluated from age 1 to age 3 and for both locations. In Estrées-Mons, it ranged from 2.2 to 5.0 when considering canopy height and plant maximum height at age 1 and age 2, while it ranged from 1.9 to 5.5 for plant stem number, plant circumference and aboveground biomass yield when considering all ages. In Orléans, the ranges of this ratio were relatively similar concerning the same traits (Table 2).
In contrast to Estrées-Mons, the biomass composition traits evaluated in Orléans showed large interaction variances according to the climate effect ( Fig. 3b; Fig. S2b), as highlighted by the ratio ranging from 1.0 to 5.7 (Table 2). This ratio was especially low in 2017, which contributed to a decrease in heritability values compared to 2018 ( Fig. 3b; Fig. S2b; Table 2; Table S1).
Accordingly, most of the production traits showed large interaction variances according to age and climate in both locations, while all composition traits showed large interaction variances according to climate in Orléans. Therefore, the corresponding heritability values were lowered. Pi ) heritability values assessed for each age according to the climate effect and considering locations separately. Each heritability value depended on the associated genetic variance, genotype × climate interaction variance, and residual variance. The genotype × climate interaction variance between G1 and G2 allowed to know the climate effect for each age. For example, plants of age 1 were evaluated for G1 in 2015 and G2 in 2016. These heritability values (right y-axis) and variance components (left y-axis) were assessed for biomass production (a) and composition (b) traits. See "Materials and Methods" section for trait name. Blanks in the plots are due to nonavailable data

Years Reflected Differences in the Climatic Conditions that Occurred During the Growing Season
For each location, a principal component analysis (PCA) was conducted to identify the weather variables that describe the climatic conditions over time. In Estrées-Mons, the first two principal components (PCs) accounted for 87.5% of the total variation in the correlation matrix, with PC1 explaining 75.2% and PC2 explaining 12.3%. In Orléans, those values were respectively 91.9%, 78.8%, and 13.1% (Fig. 4a, c).
For both locations, eight weather variables had high positive loadings (> 0.9) on the right-hand position of the PC1 first principal component (Fig. 4a, c). Those indicators were related to temperature (MinT, MeanT, MaxT for the air, and MinTs for the soil), Penman potential evapotranspiration (ETPP), photosynthetically active radiation (PAR), and vapor-pressure deficit (VPD and maximal VPD_M). They were also highly correlated with each other. MeanH had a high negative loading (< − 0.9), on the left-hand position of the PC1 first principal component (Fig. 4a, c). For both PCA, the PC2 second principal component was explained by the sum of precipitation calculated for each period and each year (Prec_S) as well as by the corresponding mean precipitation (Prec_M). Apart from irrigation (not taken into account in the calculation of the indicators), there was no real effect of climate as precipitation explained only 12 to 13% of the variability in climatic conditions.
As for Estrées-Mons, the clustering of the coordinates on the first two PCA components of the individuals (i.e., periods according to year) divided them into six distinct clusters (Fig. 4b). to include 4 years of a single period, which corresponded to Period 5 (from senescence to the date of the first frost day of the year) and Period 2 (from emergence to floral transition), respectively. Cluster 1 included the three available years for Period 1 before emergence. The year 2018 was very warm for Period 2, Period 4 and Period 5, which contrasted with all other years for each of these periods, especially in 2017 (Fig. 4b). This was particularly observed for two weather variables related to air and soil temperatures (Table S3a). Therefore, these contrasted climatic conditions between those years, for most of the periods related to the yield-building phase, partially accounted for the higher genotype × climate interaction variances observed for plant stem number, plant circumference, and aboveground biomass yield.
In Orléans, the years of Period 1 and Period 5 were spread over two clusters (Fig. 4d). There was indeed an overlap in the definition of the two periods: Period 5 of year N was defined from senescence in year N to winter harvest in year N + 1 while Period 1 in N + 1 from senescence in year N to emergence in year N + 1. Three clusters included the 4 years of a single period. The climatic conditions that occurred in 2015 were relatively rainy before floral transition and after panicle heading, and very warm and dry during the floral transition (i.e., Period 3) (Table S3b). Such conditions during the floral transition were highly contrasted with the year 2016, for which lower temperatures and higher precipitation were observed than for other years. For Period 5, the year 2016 also contrasted with years 2015 and 2017, with lower precipitation and temperatures (Fig. 4d, Table S3b). Thus, the climatic variations between 2016 and other years were likely to affect plant growth differently and contributed to explaining the genotype × climate interaction, which was especially highlighted for plant stem number, plant circumference and aboveground biomass yield at ages 1 and 2.

Except for Hemicelluloses, Strong Correlations Were Positive Within Biomass Production Traits and Biomass Composition Traits, while They Were Intermediate Between Both Types of Traits.
The genetic and phenotypic correlations were calculated in both locations, for each year-using age effect modeling (Model 2)-and each age-using climate effect modeling (Model 3). As they were relatively similar in all cases, they were illustrated in two cases: in the first case with year 2017 in Estrées-Mons (Fig. 5), the correlations corresponded to the relationships between the traits according to the age effect modeling (between age 3 for G1 and age 2 for G2) and in the second case, with age 3 in Orléans (Fig. 6), they The biomass production traits were positively correlated with each other, whether at a genetic or at a phenotypic level (Figs. 5 and 6). The strongest correlations were observed between both plant heights, canopy height and maximum height (above 0.6), between the canopy height and aboveground biomass yield (above 0.7), and above all, between the plant circumference (C50_cm) and aboveground biomass yield (above 0.9).
The biomass composition traits were positively correlated with each other, except hemicelluloses and ash (Figs. 5 and 6). Relatively strong correlations (above 0.49) were mainly highlighted, except for lignin (%CW) which was less correlated with NDF, cellulose in %DM, and cellulose in %CW (from 0.29 to 0.49). The correlations were sometimes lower when the biomass composition traits were expressed in %CW.
Most of the biomass production traits were moderately and positively correlated with the biomass composition traits, the threshold of 0.5 never being exceeded for phenotypic correlations. The relationships were sometimes stronger at a genetic level than at a phenotypic level, especially with NDF, ADF, and cellulose (%DM and %CW), for which they were sometimes just above the threshold of 0.5. Excepting ash content, the hemicellulose content (in %DM or %CW) was negatively correlated with all the other traits, including biomass production traits.

Discussion
The genetic parameters of biomass production and composition traits estimated for the M. sinensis population, highlighted intermediate to high individual plant broad-sense heritability values ( H 2 sl ) according to the age effect per year and the climate effect per age, for each location considered. These values mainly increased over age, depending on genetic, genotype × age, and genotype × climate interaction variances as well as residual variances. Contrasted climatic conditions between the years and identified throughout the plant cycle for each year made it possible to explain the high genotype × climate interaction variance highlighted for some traits. Finally, genetic and phenotypic correlations between biomass production and composition traits were assessed according to specific years and ages. Strong and positive correlations were observed within biomass production traits and within biomass composition traits, except for Genetic and phenotypic correlations (respectively above and below the diagonal) calculated in Orléans for the age 3. Biomass production and composition traits were considered. See "Materials and Methods" section for trait name. Displayed correlations are significant with a p-value < 0.05, while non-significant correlations are in blanks hemicelluloses, while intermediate correlations were highlighted between these two types of traits.
Therefore, three main points will be discussed in this section: (1) the contrast between young and old plant ages for genetic parameters thanks to the staggered-start design and its consequences on the miscanthus breeding cycle, (2) higher heritability values for biomass production traits than biomass composition traits, and (3) the intermediate correlations observed between biomass production and composition traits.

The Contrast Between Young and Old Plant Ages for Genetic Parameters Thanks to the Staggered-Start Design and Its Consequences on Breeding Cycle Duration
Using our staggered-start design, the genotypic and environmental variability was assessed based on genetic variance, genotype × age, and genotype × climatic condition interaction variances and residual variance. For most traits in both locations, the genetic variance was generally lower for the first year of evaluation according to the age effect (Model 2) or for the first age of evaluation according to the climatic condition effect (Model 3). This first age corresponded to young plants. The genetic variance often increased throughout the years and ages: this can be due to the fact that, while genotypes grow, their differences increase, resulting in higher genetic variability for each trait studied. According to each model, genotype × age and genotype × climate interaction variances were highlighted and were larger for young plant ages than for older plant ages, especially for biomass production traits. This was illustrated by the lower ratios of genetic variance to genotype × age interaction variance and the lower ratios of genetic variance to genotype × climate interaction variance. Moreover, the competition effects decreased with age while the genetic effects increased ( Fig  S1). However, most of the production traits, especially plant stem number, plant circumference, and aboveground biomass yield, also showed large genotype × age interaction variances for successive years and large genotype × climate interaction variances for successive ages, in both locations. The population was indeed composed of genotypes showing different growth strategies: for instance, some may have more small stems, while others fewer large stems. Besides that, their strategies can vary according to the ages and climatic conditions. The presence of greater genotype × age interaction variances for successive years means that we observed an age effect and that this age effect varied among genotypes. Moreover, these variations were more pronounced for successive years. This was observed for plant stem number, plant circumference, and aboveground biomass yield. This means that the values of these variables did not increase in the same manner for the genotypes in the population. The presence of large genotype × climate interaction variances for successive ages implies that differences among progeny were amplified by the more pronounced climatic conditions in 2017-2018. We indeed found that the year 2018 was very warm for most of the growing season (Period 2, Period 4, and Period 5), which contrasted with all other years for each of these periods, especially in 2017. Biomass composition traits showed large genotype × climate interactions for ages 2 and 3 studied in Orléans. These interaction variances rarely overtook the genetic variances in our study, meaning that genetic determinism is mainly led by stable genetic effects for each trait. Nevertheless, those interactions revealed that some genotypes evaluated in the population responded differently: (1) for different plant ages (i.e., the age effect) in the same year, with corresponding climatic condition or (2) in different climatic conditions (i.e., the climate effect) for the same age. Segura et al. [50] dissected apple tree architecture into ontogenetic and environmental effects: by using a staggered-start design, they related genotype × age interactions to ontogenetic effects and genotype × year interactions to climate effects. We found results in accordance with their approach. Dong et al. [47] considered the genotype × year interaction to assess the broad-sense heritability of miscanthus production traits and found a substantial effect for plant height, plant circumference, and culm-associated traits in one of their mapping populations. Kar et al. [69] evaluated the broad-sense heritability of miscane production traits (i.e., hybrids between sugarcane and miscanthus) and found substantial genotype × year variances, especially for the stem length and number of stems. Our study is consistent with their findings, even the genotype × year interaction variances presented in these last studies were not decomposed into age and climate effects as in the present study.
Intermediate to high individual plant broad-sense heritability values and high progeny-mean broad-sense heritability values were found for biomass production and composition traits in both locations. The low genetic variances found for the first year or first age evaluated, in combination with large genotype × age interactions or large genotype × climate interactions, tend to decrease associated heritability values. Thus, trait heritability values mostly increase throughout years and ages. These observations throughout years or ages were made possible thanks to the staggered-start design established in each location. Each single year was evaluated with its corresponding climate (i.e., climatic conditions that occurred during the year considered), and the age effect was modeled based on the different ages between G1 and G2 groups. Moreover, each age was evaluated and the climate effect was modeled according to the different climatic conditions between the year considered for G1 and the year considered for G2 (G1 and G2 having the same age). Our results are consistent with the majority of previous studies: Slavov et al. [41,42] evidenced an increase in broad-sense heritability values for most biomass production and composition traits, by studying an M. sinensis population during 3 years. Clark et al. [45] evaluated accessions from different miscanthus species in six locations: an increase in broad-sense heritability values was observed for most biomass production traits, when considering broad-sense heritability for all locations between the second and third year. Gifford et al. [46] and Dong et al. [47] also reported higher broad-sense heritability values from year 2 to year 3 after establishment, when considering biomass production traits in M. sinensis biparental populations. However, Van der Weijde et al. [48] found that the broad-sense heritability values of biomass composition traits decreased between the second and third year after the establishment on an M. sinensis biparental population. The evolution of heritability estimates is specific to each population and environment tested [37]. As a perennial crop, miscanthus is subjected to changes in conditions from one year to the next, depending on each environment considered. However, such changes highlighted between years can be due to plant age or climatic conditions: to our knowledge, the present study is the first to allow the assessment of the broad-sense heritability of traits throughout years (i.e., each year with its corresponding climate when the age effect was taken into account) or throughout ages (i.e., each age when the climate effect was taken into account). Broad-sense heritability values are important to consider in a plant-breeding program, in order to start initial selections and provide genetic progress [40]. Because higher broad-sense heritability values were highlighted throughout years and ages, progress through selection would be greater as plants got older, especially for biomass production traits. Accordingly, this increases the duration of the selection process.
This study highlights the advantage of the staggered-start design in decomposing the year effect into climate and age effects to efficiently assess the genetic parameters in the perennial Miscanthus sinensis. In particular, it shows that the genotype × age and genotype × climate interaction variances are larger for young plant ages than for old plant ages, especially for biomass production traits. In addition, genetic variance as well as broad-sense heritability values increases across age spectrums. Such information must therefore be considered as important to miscanthus breeding.

Higher Heritability Values Were Observed for Biomass Production Compared to Biomass Composition Traits
Broad-sense heritability values and their means over years or ages were mainly higher for biomass production traits than biomass composition traits. For individual plant broad-sense heritability means over years or ages, biomass production traits were moderately to highly heritable, ranging from 0.42 to 0.62, while biomass composition traits were moderately heritable and ranged from 0.26 to 0.47. According to the progeny-mean broad-sense heritability means, both types of traits were highly heritable and the trends between them were similar to those of individual plant broad-sense heritability: means ranged from 0.67 to 0.88 for biomass production traits and from 0.60 to 0.80 for biomass composition traits. Contrary to our study, previous studies generally showed that biomass composition traits were more heritable than biomass production traits. According to their M. sinensis population evaluated over the course of 3 years, Slavov et al. [41] found higher broad-sense heritability means for biomass composition traits than biomass production traits, with means of 0.67 and 0.60, respectively. Petit et al. [70] studied the genetic variability of morphological, flowering, and biomass quality traits in hemp (Cannabis sativa L.) in three locations during one season: biomass composition traits were generally more heritable than biomass production traits.
Thus, the lower broad-sense heritability values observed for the biomass composition traits in our study are probably due to lower variability caused by genetic differences among the progeny (i.e., genetic variance) than for biomass production traits. In addition, it appears that biomass composition traits are frequently more affected by environmental effects (especially when focusing on the residual variance compared to other variances).
According to the staggered-start design, the means of heritability values over years and ages are higher for biomass production traits than for composition traits. This shows that the potential genetic improvement of biomass production traits is greater in the present study. Biomass composition traits are more affected by multiple-environment effects, which hamper the expression of the genetic effect in both locations studied. Thus, the improvement of biomass composition traits has to be carefully considered in breeding programs, in order to meet the specific requirements of each bioconversion routes.

Intermediate Correlations Were Highlighted Between Biomass Production and Composition Traits
The genetic and phenotypic correlations showed that biomass production and composition traits were moderately correlated in both locations. Even positive and strong correlations were highlighted within the biomass production traits and within most of the biomass composition traits, for each year according to the age effect (Model 2) or each age according to the climate effect (Model 3), the correlations between both types of traits were intermediate. In accordance with our study, Slavov et al. [42] reported intermediate correlations between most of the biomass production and composition traits, considering each evaluated year of their M. sinensis population. Petit et al. [70] also found intermediate correlations between biomass production and composition traits in hemp. While it would be possible to easily improve biomass production traits or most of the biomass composition traits separately, the improvement of both types of traits together would be challenging due to their intermediate correlations. In addition, hemicelluloses would be difficult to improve with the other biomass production and composition traits, due to their negative correlations (in %DM and %CW). This was also reported by Slavov et al. [42] and Van der Weijde et al. [48]. This would interfere with downstream processing of biomass, as hemicellulose concentrations influence the conversion processes, including conversion into bioethanol or polymer composites: in miscanthus, Brancourt-Hulmel et al. [71] showed that the best composite performances were obtained with stems that showed the lowest concentrations of hemicellulose in association with higher concentrations of cell wall, lignin, p-coumaric acid, and cellulose. In contrast, Belmokhtar et al. [59] found that more digestible genotypes contained higher amounts of hemicellulose and lower amounts of cellulose and lignin.
Morphological traits contribute to the improvement of biomass production in order to enhance yield, regardless of the age or climate effect considered. However, the improvement of biomass yield for newly bred plants is not concomitant of the improvement of their biomass hemicellulose content. The improvement of biomass production traits does not generate a substantial increase in the other desirable components, due to the intermediate genetic correlations observed. Thus, the adjustment of biomass composition traits regarding age and climatic conditions must be carefully considered by miscanthus breeders, in order to select a biomass that is suited to the end-use targeted.

Conclusion
In this study, intermediate to high individual plant broadsense heritability and high progeny-mean broad-sense heritability values were assessed for biomass production and composition traits. These heritability values changed in relation to the variance components estimated for each year according to the age effect modeling (Model 2) or each age according to the climate effect modeling (Model 3), when considering each staggered-start design per location. Lower genetic variances associated with substantial genotype × age and genotype × climate variances were observed for young plant ages, which tended to decrease heritability values, especially for biomass production traits. Such large interaction variances, compared to genetic variances, were highlighted according to the ratio of genetic variance to genotype × age interaction variance and the ratio of genetic variance to genotype × climate interaction variance. Heritability values generally increased through successive ages and years for most of the traits, which suggested an improved genetic progress when plants get older and implied an increase in the duration of the selection process. However, these values were sometimes lowered by large genotype × age interaction variances for successive years, especially in the case of the plant stem number, plant circumference, and aboveground biomass yield in both locations. Similarly, large genotype × climate interactions for successive ages were found for most biomass production traits in both locations and biomass composition traits in Orléans.
Overall, biomass production traits were more heritable than biomass composition traits, which means that the response to selection would be better for biomass production traits in our population. Also, even high genetic and phenotypic correlations were highlighted among biomass production and composition traits separately, intermediate correlations were diagnosed between biomass production and composition. In contrast, hemicelluloses and ash contents were negatively correlated with other traits. Thus, miscanthus breeders have to consider these parameters in order to design breeding programs that follow the requirements of each bioconversion route. The next step of this study will be QTL detection for biomass production and composition traits, while considering the age and climate effects evaluated thanks to the staggered design in each location.