Characterization of phenology, physiology, morphology and biomass traits across a broad Euro‐Mediterranean ecotypic panel of the lignocellulosic feedstock Arundo donax

Giant reed (Arundo donax L.) is a perennial rhizomatous grass, which has attracted great attention as a potential lignocellulosic feedstock for bioethanol production due to high biomass yield in marginal land areas, high polysaccharide content and low inhibitor levels in microbial fermentations. However, little is known about the trait variation that is available across a broad ecotypic panel of A. donax nor the traits that contribute most significantly to yield and growth in drought prone environments. A collection of 82 ecotypes of A. donax sampled across the Mediterranean basin was planted in a common garden experimental field in Savigliano, Italy. We analysed the collection using 367 clumps representing replicate plantings of 82 ecotypes for variation in 21 traits important for biomass accumulation and to identify the particular set of ecotypes with the most promising potential for biomass production. We measured morpho‐physiological, phenological and biomass traits and analysed causal relationships between traits and productivity characteristics assessed at leaf and canopy levels. The results identified differences among the 82 ecotypes for all studied traits: those showing the highest level of variability included stomatal resistance, stem density (StN), stem dry mass (StDM) and total biomass production (TotDM). Multiple regression analysis revealed that leaf area index, StDM, StN, number of nodes per stem, stem height and diameter were the most significant predictors of TotDM and the most important early selection criteria for bioenergy production from A. donax. These traits were used in a hierarchical cluster analysis to identify groups of similar ecotypes, and a selection was made of promising ecotypes for multiyear and multisite testing for biomass production. Heritability estimates were significant for all traits. The potential of this ecotype collection as a resource for studies of germplasm diversity and for the analysis of traits underpinning high productivity of A. donax is highlighted.

. A. donax lignocellulosic biomass feedstock can be readily used as solid biofuel in direct combustion for heat generation, or for bioethanol production from alcoholic fermentation of lignocellulose material (Pilu et al., 2012;Silva, Schirmer, Maeda, Barcelos, & Pereira, 2015). Based on a comprehensive utilization indexes obtained from an evaluation system that takes into consideration numerous factors deemed essential for the conversion of lignocellulosic biomass to bioethanol, a recent study revealed that A. donax is a highly promising lignocellulosic feedstock for producing bioethanol compared to corn stalks, switchgrass, pennisetum and silvergrass (Kou, Song, Zhang, & Tan, 2017). Moreover, A. donax has been used as a potential feedstock for thermo-chemical conversion to biochar, bio-oil and gas for energetic purposes (Saikia, Chutia, Kataki, & Pant, 2015). A. donax cropping system is important not only because it provides farmers with an innovative opportunity but also for mitigating greenhouse gas (GHG) emissions, controlling soil erosion and phytoremediating contaminated soils (Cosentino, Copani, Scalici, Scordia, & Testa, 2015;Fagnano, Impagliazzo, Mori, & Fiorentino, 2015;Nsanganwimana, Marchand, Douay, & Mench, 2014;Zucaro et al., 2015). Furthermore, A. donax is a resilient and very productive species. It has been shown that the photosynthetic capacity of A. donax in full sunlight is high compared to other C3 species, and comparable to C4 bioenergy grasses (Rossa, Tüffers, Naidoo, & von Willert, 1998). Because of these advantages, A. donax is expected to play a major role in the provision of perennial lignocellulosic biomass across much of Europe. However, successful commercialization of renewable biofuel production depends on the sustainability and biomass yield of lignocellulosic crops.
Field trials of A. donax in central Italy showed an average yield of above-ground biomass of about 40 Mg dry matter (DM) ha −1 yr −1 , with a peak of about 50 Mg DM/ha at the end of the second and third years of growth over a production cycle of 10 years (Angelini et al., 2009). Similarly, in a multi-year study in Alabama, A. donax achieved an average yield of 35.5 Mg DM ha −1 yr −1 (Ping, Bransby, & van Santen, 2014). However, biomass yield can be doubled at the end of the same growth cycles in northern Italy under fertirrigated conditions (Borin et al., 2013). Biomass yields in these experiments using small plots may be an overestimate of the yields achievable in competitive commercial planting. While increasing yield from the first to the third one-year coppice rotation was often observed in A. donax (Borin et al., 2013;Cosentino, Copani, Patanè, Mantineo, & D'Agosta, 2008;Nassi o Di Nasso, Roncucci, & Bonari, 2013), different conclusions for the yield trend from the third year onward can been drawn (Lewandowski et al., 2003). Although the reported biomass yield for A. donax is high under irrigated conditions, the large-scale plantations of this crop remain a major limitation due to the high cost of propagation (Ceotto & Di Candilo, 2010). Biomass yields may also be considerably lower on marginal land. For example, on a sandy loam with limited nutrient availability, A. donax required 3 years to accumulate 20 Mg DM/ha, and build a rhizome mass of 16 Mg DM/ha (Nassi o Di Nasso et al., 2013).
Although A. donax can flower and produce seed under some conditions, as yet there is no evidence that the seeds are viable (Barney & Ditomaso, 2008;Haddadchi, Gross, & Fatemi, 2013;Lewandowski et al., 2003;Mariani et al., 2010;Pilu et al., 2012). Thus, the natural variability in existing populations may occur following spontaneous mutations selected by natural pressure (Cosentino, Copani, D'Agosta, Sanzone, & Mantineo, 2006). Despite the low genetic diversity revealed in a few studies using molecular markers (Ahmad, Liow, Spencer, & Jasieniuk, 2008;Hardion, Verlaque, Baumel, Juin, & Vila, 2012;Hardion, Verlaque, Saltonstall, Leriche, & Vila, 2014;Mariani et al., 2010;Pilu et al., 2014), heritable phenotypic differences have been described in growth-related traits (Cosentino et al., 2006;Pilu et al., 2014). For this reason, ecotype selection remains important for genetic improvement of A. donax (Corno et al., 2014;Pilu et al., 2013). Further, it is important that traits which contribute to harvestable biomass yield in grasses such as A. donax are well understood in order to drive genetic gain for accelerated domestication as discussed for other biomass crops . A thorough understanding of the interactions between traits, both positive and negative correlations, is required because selection for one trait can have profound effects on another (Farrar et al., 2012).
Diversity in growth-related traits has been reported using ecotypes collected from regions within the Mediterranean basin (Cosentino et al., 2006;Sánchez et al., 2015) including some studies of ecotypes from a wider range of provenances (Amaducci & Perego, 2015;Pilu et al., 2014); however, a comprehensive study of different parameters that may interact with yield (i.e. phenology, reed morphology, leaf functional traits and physiology) in a large panel of A. donax ecotypes collected from diverse environments has not been previously reported.
This study aimed: (a) to characterize a broad Euro-Mediterranean panel of 82 A. donax ecotypes by assaying 21 biomass, morphology, physiology and phenology traits in a common garden experiment and to understand how these may correlate to each other; (b) to identify traits that may have contributed to differences in biomass production in field conditions; and (c) to select potential high yielding clusters through a multivariate approach. This study design therefore also serves as a key resource for parametrizing biomass production models of A. donax, an emerging lignocellulosic crop that would ultimately be an important contributor to any substantive EU renewable energy strategy that is based on lignocellulosic biomass supply.

| Collection of plant material
From 2006 to 2010, eighty-two Euro-Mediterranean ecotypes of A. donax were collected for use in the study. To avoid clonal ramets collection, sites were at least 4 km and at most 474 km apart (mean distance = 75 km) (Figure 1). The local geography of collection sites included riverbanks and riverbeds of the Mediterranean basin, uncultivated fields and marginal areas. At each site rhizomes were collected from a single clump that was considered to be a distinct ecotype.

| Experimental site
In May 2011, rhizomes were planted in a common garden experiment at the experimental field of the "Alasia Franco Vivai," Savigliano in northern Italy (SAV, 44°36′N, 07°38′ E, 340 m above sea level). The soil is alluvial with sandy loam texture. The regional climate is fully humid and warm temperate with warm summers, according to the Köppen-Geiger climate classification system (Kottek, Grieser, Beck, Rudolf, & Rubel, 2006). The average annual precipitation and temperature was 848.2 mm and 10.5°C, respectively. During this study, spring was wetter and colder, and summer was drier than the 5-year seasonal average ( Figure 2). During the sampling period and specifically from 19 July to 23 August, the weather was extremely dry and warm, with a daily average maximum temperature of 29.4°C and 12 mm of accumulated rainfall ( Figure 2) as assessed by a weather station located 1 km from the common garden experimental field.

| Experimental design
Each ecotype was planted in a 9 m 2 plot (3 × 3 m) using five homogeneous rhizomes planted in 2 × 2 m spacing with one rhizome in the centre of that space (quincunx system of planting). Plots were harvested at the end of each F I G U R E 1 Map showing the collection sites of the 82 ecotypes of the bioenergy grass Arundo donax. Solid circles indicate the locations where samples were collected growing season. In 2011 and 2012, the site was fertilized and irrigated by furrowing every 20 days while no irrigation water or fertilizer were applied during the 2013 growing season. During the 3 years, no serious incidents of insects or diseases were observed, and no pesticides or fungicides were applied.

| Phenotypic data collection
Arundo donax ecotypes start the maturity phase by the third one-year coppice rotation (Angelini et al., 2009), making phenotypic comparisons more meaningful from the third growth year onwards. Data were collected from July to November 2013, and traits and their abbreviations are listed in Table 1. For each ecotype, the main plot was divided into four 1.5 × 1.5 m subplots and three of the subplots were randomly selected as biological replicates for all measurements. Empty subplots where the original rhizomes died after establishment were excluded from the selection. No centre rhizome pieces died therefore the contributions to the subplots were always the same. Consequently, at the beginning of the growing season, the mean number of available planted clumps per ecotype was 4.51 (Supporting information Table S1).

| Phenology
Flowering phenological data were recorded in autumn 2013, from 19 September [day of the year (DOY) 262] to 30 October (DOY 303). The timing of five discrete phenological stages was assessed daily by visual inspection following a new scoring system. In detail, Phase 0 (Ph0)-absence of flowers; Phase 1 (Ph1)-onset of flowering with the appearance of the first small inflorescence at the singlestem level; Phase 2 (Ph2)-the first inflorescence had reached at least 30 cm length; Phase 3 (Ph3)-the first inflorescence was completely open; and according to Cosentino et al. (2006), Phase 4 (Ph4)-50% of the reeds showed completely open inflorescences.

| Physiology
Mid-day stomatal resistance (Rs, s/cm) was measured on July 10th (Rs july ) and August 22nd (Rs august ) with a dynamic diffusion porometer (AP4; Delta-T devices, Cambridge, UK). Data were taken from the abaxial mid-leaf position of three healthy, fully mature and fully expanded leaves (5th from the top) from three randomly chosen reeds per ecotype between 11:00 and 13:00 hr. Before each session of measurement, the porometer was calibrated to field temperature and humidity.
The same leaf previously used for stomatal resistance measurement (Rs august ) was collected, oven-dried at 75°C to constant weight, and stored for subsequent carbon isotope discrimination (Δ, ‰) analysis. Leaves were ground into a homogenous fine powder. Samples were combusted in an elemental analyzer (Thermo Flash EA 1112 Series, Bremen, DE); CO 2 and N 2 were separated by chromatography and directly injected into a continuous-flow Isotope Ratio Mass Spectrometer (Thermo-Finnigan Delta XP, Bremen, DE) for isotope analysis. Peach leaf standards (NIST 1547) were run every six samples. 13 C/ 12 C isotope ratio (δ 13 C) was calculated according to the definitions described by Farquhar, Ehleringer, and Hubick (1989):  where R sample /R standard were referred to a Pee Dee Belemnite standard.
And Δ was calculated according to Farquhar et al. (1989): where the δ 13 C atm was assumed −8‰.

| Canopy and leaf morphology
Leaf area index (LAI, m 2 /m 2 ) was indirectly calculated in each subplot by measuring incoming photosynthetic active radiation (PAR) below the canopy and in an open area using the AccuPAR LP-80 ceptometer (Decagon Devices, Pullman, WA, USA). Three mature leaves of each ecotype were collected from the middle part of the main axis of the longest reed (stem) in each subplot (Hardion et al., 2012). Leaves were immediately transported to the laboratory and photographed using a digital camera. Before image analysis, pictures were corrected for geometric distortions caused by the camera lens using filters available in Adobe Photoshop Elements 12. Leaf length (LfL, cm), leaf width (LfW, cm) and leaf area (LfA, cm 2 ) were then measured using the LAMINA software (Bylesjö et al., 2008). The ratio between LfL and LfW (LfL_LfW) was calculated and represents an independent shape variable that has been used extensively in leaf morphometric analysis (Aravanopoulos, T A B L E 1 Abbreviation, units and description of phenological, physiological, morphological and biomass and growth traits

Traits Abbreviation Units Description
Phenological traits  Sun et al., 2014). All leaves were then oven-dried at 60°C for 72 hr to constant weight, and the specific leaf area (SLA, cm 2 /g) was determined as the ratio between LfA and leaf dry weight.

| Growth and biomass production
To assess growth-related traits and biomass production, all ecotypes were harvested at the end of the growing period in November 2013, by cutting at 5 cm aboveground level. In each subplot, the weight of the above-ground fresh biomass was measured and the number of reeds counted to calculate stem density (StN, n/m). Five stems were randomly selected, and the following measurements were performed: height of the stem (H, cm) from the base-node to the top-node without the inflorescence (if present); number of nodes per stem (StNd, n); diameter of the stem in three positions above the ground: (a) 20 cm (Db, mm), (b) at half of H (Dm, mm) and (c) 20 cm below the stem's apex (Da, mm); and fresh weight (Cosentino et al., 2006). Six sub-samples of whole stems per ecotype, two for each selected subplot, were weighed and oven-dried at 75°C until a constant weight was obtained to determine the stem dry mass (StDM, g DM) and the total biomass production (TotDM, Kg DM/m).

| Statistical analyses
Statistical analyses of the data were performed using R software (Version 2.15.2, A Language and Environment Copyright, 2012). Original values were tested for transformation using the Box-Cox procedure (Venables & Ripley, 2002) to ensure homoscedasticity and normality of residuals. Besides the determination of means and ranges, the coefficient of variation (CV), calculated as the ratio of the trait's standard deviation to its mean value and reported as a percentage (%), was used to compare the variability of traits independently of the units of measure. To test differences among ecotypes and replicates of each ecotype, variables were analysed using the analysis of variance (ANOVA) model, as follows: where Y ijk is the individual value, μ is the general mean, E i is the effect of the i th ecotype (random), R j is the effect of the j th replicate (fixed) and ε ijk is the residual error term. With regard to the stomatal resistance, measured at two time-points, the following model was used to test time (T) and ecotype by time interaction (E × T) effects: In this case, T was the time effect considered as random, and E × T was the interaction effect considered as random.
Because of the absence of a replicated block design, R could not be corrected when significant to eliminate micro-environmental variations, and thus, this source of variation was removed to make the ANOVA models simpler. When the analysis of variance showed significant E × T interaction (P ≤ 0.05), the change in ranking of ecotypes was evaluated by the Spearman rank coefficient (ρ) on ecotypic means.
Pearson's correlation coefficient (r) was used to assess the relationship among all traits, using mean values for each ecotype. All statistical tests were considered significant at P ≤ 0.05.
Multiple linear regression analysis was used to analyse traits affecting the amount of dry biomass production and to select a sub-set of variables to be used for cluster analysis. Hence, the analysis was performed taking into consideration TotDM as dependent variable and all other traits (except flowering stages for which only non-replicated measurements were obtainable) as independent variables. As the residuals have to be normally distributed, the Box-Cox procedure (Venables & Ripley, 2002) was performed to search for the most appropriate transformation of the dependent variables of the linear regression model.
Hierarchical cluster analysis was performed to group ecotypes with similar characteristics on a multivariate basis, and the Euclidean distance was used to determine the difference between the groups. The analysis was performed using standardized values of all traits. Ward's minimum variance method was used to construct the clusters, and the actual number of clusters was discerned using "NbClust" package in R (Charrad, Ghazzali, Boiteau, & Niknafs, 2013), which provides most of the popular indices for determining the optimal number of clusters. ANOVA model was used to test the null hypothesis that clusters (Clu) were equal, as follows: where Y ij is the individual value, μ is the general mean, Clu j is the effect of the j th cluster (random) and ε ij is the residual error term. LSD post hoc test was used to compare significant differences among mean cluster values.
To describe ecotypes with their interrelations and at the same time to create a graphical display of the clusters, the "CLUSPLOT" algorithm of Pison, Struyf, and Rousseeuw (1999) was used. Finally, the dimension of the data was reduced by principal component analysis (PCA), then "CLUSPLOT" displayed the bivariate plot of the objects relative to the first two principal components (components 1 and 2) (Pison et al., 1999).

| Broad-sense heritability
Broad-sense trait heritability on an individual basis (H 2 ) was estimated as follows: where σ 2 G and σ 2 ɛ are the genetic and residual variance components, respectively, which were obtained using the restricted maximum-likelihood (REML) procedure (Holland, Nyquist, & Cervantes-Martinez, 2003). H 2 was evaluated for all the studied traits, except for flowering stages for which replicated measurements were not viable.

| Univariate analysis
A significant level of phenotypic variation among the 82 different ecotypes of A. donax was observed in all physiological, canopy, growth and biomass traits measured (Table 2). Flowering phenology of all ecotypes was recorded using a new scoring system developed in this study that covers the onset of flowering (Ph1) to the end of flower development (Ph4). Flowering was evaluated at plot scale, without replication (Table 1); therefore, statistical differences between ecotypes and blocks for flowering phenology trait could not be evaluated (Table 2). Average flowering time of Ph1 and Ph4 at ecotype level were DOY 265 (22 September) and DOY 284 (11 October), respectively, giving an average of 19 days between Ph1 and Ph4. Transition between the first 3 phases was rapid, for example taking an average of 4 days from Ph1 to Ph2 and from Ph2 to Ph3, while an average of 11 days were required to complete the last flowering phase (Figure 3; Table 2). Due T A B L E 2 General means, ecotypic range of variation, coefficient of variation (CV, %) and level of significance of ecotype (E) and replicate (R) effects for phenological, morphological, physiological and growth traits in Arundo donax. Broad-sense heritability (H 2 ) for physiological, morphological, growth and biomass traits. See Table 1  to high phenotypic variability, there were some notable outliers, on the first day of flowering assessment two early flowering ecotypes were in Ph4 (ecotypes "8″ and "35″).
Conversely, ecotypes "14," "16" and "34" did not flower at all. Of the ecotypes that did reach Ph4 the most delayed were "30," "37," "50" and "52" (29 October). A significant ecotypic effect was detected in all traits examined, whereas the block effect appeared to be statistically significant only for Rs july , LAI, TotDM and StN (Table 2). Traits with highest variation among ecotypes were Rs july and Rs august (CV of 28.20% and 27.76%, respectively), ranging from 0.57 to 1.68 s/cm in July, and between 0.8 and 2.69 s/cm in August (Table 2). Notably, Rs august was measured under moderate drought conditions (i.e. 36 days of high temperatures, 21.3°C of mean temperature and maximum daily average temperature of 29.5°C) and low rainfall amounts concentrated in 2 days (5.2 mm in July, 29 and 3.2 mm in August, 8). Accordingly, ANOVA analysis highlighted significant differences between the two time-points (P ≤ 0.001) with ecotypes that reacted differently to different times (P ≤ 0.001 for E × T; Figure 4a,b). In addition, significant changes in the ranking of ecotypes (ρ = 0.363, P ≤ 0.001) were observed. There was an inverse relationship between Rs august/Rsjuly ratio and absolute Rs july, and ecotypes with the lowest Rs july values had higher Rs august values, up to 3 fold than Rs july (Figure 4c). The physiological variation in Δ found within ecotypes of our A. donax collection ranged from 19.43‰ to 22.23‰ (Table 2), with a difference of 2.8‰ observed between the two most divergent ecotypes, "11" and "30." Moreover, the HSD post hoc test revealed significant differences among ecotypes "11," "77" and "82," which showed the lowest Δ values, and ecotypes "49" and "30," which showed the highest Δ values ( Figure 5; Table 2).
Leaf and canopy traits showed medium to high variability (Table 2). LAI ranged from 3.89 to 8.14 m 2 /m 2 and was the most variable of the canopy and leaf morphology group of traits. Values of LfL ranged from 48.25 to 65.30 cm, and a CV of 6.27% accounting for less than half of the variability observed for LAI. Likewise, LfW was also less variable than LAI, but CV was higher than LfL (8.56%). Leaf shape, measured by LfL_LfW, ranged between 8.59 and 11.62 indicating leaves were approximately 10 times longer than wide.
At the end of the growing season, differences in TotDM among ecotypes were significant, with means ranging from 5.26 for ecotype "87" to 15.21 kg DM/m 2 for ecotype "49" (Table 2 and Supporting information Table S1). In addition, TotDM represented one of the most variable traits with a CV of 23.18% that was twice as large as for H (10.26%) and for diameter measured at different positions (approximately 10% for Db, Dm and Da). Similarly, high ecotype differences were observed for StDM and StN, with a slightly higher variation for the latter (22.16% vs. 25.30%, respectively).

| Phenotypic correlations
The linear relationships among traits were studied using Pearson's correlation analysis ( Figure 6). Significant positive r values among flowering phases showed that when the flowering process began, all the other stages followed through a cascade process. Moreover, the date of the latest flowering stage, Ph4, were moderately related to LAI (r = 0.705, P ≤ 0.001) and inversely related to H (r = -0.631, P ≤ 0.001). All flowering phases were significantly and positively correlated with Δ, with Ph4 and Δ (r = 0.489, P ≤ 0.001) showing the highest correlation.
Regarding stomatal resistance, there were few weak and negative relationships between Rs july and Rs august and other traits, such as flowering time, Δ, LAI, Db and Dm, whereas the two traits positively correlated with each other.
LfL and LfW were positively correlated and interestingly, leaf morphology was strongly related to StN, with a higher correlation with LfW (r = −0.610, P ≤ 0.001) than LfL (r = −0.144, P = ns).
Among growth traits, as might be expected the strongest correlations were detected between StDM and the diameter of the reed, especially Db (r = 0.851, P ≤ 0.001). TotDM, on the other hand, did not correlate with diameter of the single stem, but significant correlations were found with StN (r = 0.526, P ≤ 0.001) and H (r = 0.500, P ≤ 0.001; Figure 6).  Table 1 Table 1 for trait abbreviations. Significant correlations are coloured either in red (positive) or blue (negative) hues, while correlations that were not significant are indicated by "×." F I G U R E 4 Phenotypic variation for stomatal resistance (Rs) in a broad Euro-Mediterranean panel of 82 Arundo donax ecotypes grown in northern Italy. Ecotypic means (black dots) ± SE for Rs measured in (a) July, 10 (Rs july ) and in (b) August, 22 (Rs august ) are shown together with the collection mean (horizontal dashed line). (c) Pairwise trait plot between Rs july and Rs august /Rs july ratio. Unique identifiers of each ecotype are written below the dot, except identifiers "1," "18," "30," "36," "42," "47," "49," "55," "61," "76," "82" that are written above the dot, to avoid labels overlapping. See Table 1 for trait abbreviations. In (a) and (b) ecotypes were ordered by increasing Rs means. Statistically significant differences among ecotypes, evaluated with HSD post hoc test (P ≤ 0.05), are shown above the graph (b); however, no statistical significance means are found in (a). In (b), for clarification, for each ecotype only the first and the last letters of the group are shown

| Cluster analysis
Hierarchical cluster analysis was used to group ecotypes with similar characteristics on a multivariate basis. To reduce the number of variables to be used in the analysis, a sub-set of informative traits, that is LAI, StDM, StN, StNd, H, Db and Dm, were selected by a multiple regression analysis (adjusted determination coefficient, R 2 = 0.753,  (Table 3). Ward's hierarchical clustering method indicated that the 82 ecotypes could be divided into five groups (Figures 7 and 8; Table 4). In addition, these clusters were significantly different for all the traits studied (  (Table 5). Clu.1 consisted of 17 ecotypes and was located in the centre of the cluster biplot (Figure 8). This cluster contained ecotypes classified as late in flowering (Ph4 at the end of October), medium for growth performance (biomass and morphological traits), preserved low stomatal resistance in the two data points, and that showed the highest values of Δ (Table 4). However, the biplot indicated Clu.1 was an intermediate group with strong overlapping regions with Clu.3 and Clu.4, that is areas where memberships are ambiguous (Figure 8). The latter two clusters, comprising 25 and 19 ecotypes, respectively, were more divergent for TotDM and StNd. Clu.3 included early-flowering ecotypes (Ph4 at the end of September) that had big leaves (LfL = 59.33 cm, LfA = 222.71 cm 2 ), a high biomass yield (10.90 kg DM/m 2 ) and a low Δ (20.72 ‰). Moreover, ecotypes in this cluster were characterized by dense (331.59 g) and knotty reeds (an average of approximately 39 nodes in each reed). On the other hand, Clu.4 comprised accessions that produced few and small stems, reflecting low biomass production (6.83 kg DM/m 2 ) ( Table 4). Clu.2 and Clu.5 were two distinct groups of ecotypes that are more differentiated for reed dimensions and StN and, as a consequence, for LAI. Clu.2 was a small cluster of 6 early-flowering and high performing ecotypes that form large and tall reeds (Db and H were 24.27 mm and 493.36 cm, respectively) in clumps that were not particularly dense (StN =~35 per m 2 ), together constituting light canopies (LAI = 5.11 m 2 /m 2 ). Clu.5 included 15 lateflowering ecotypes characterized by small reeds (Db = 18.03 mm and H = 395.11 cm) and dense clumps (StN = 55 per m 2 , LAI = 6.25 m 2 /m 2 ). However, both Clu.2 and Clu.5 produced similarly high amounts of biomass that differed significantly only from Clu.4. Finally, Clu.2 and Clu.5 showed similar Δ (20.97 ‰).

| DISCUSSION
In this study, we present a comparative analysis of a large number of A. donax ecotypes collected from different Euro-Mediterranean regions throughout most of the species' native range. At the end of the third one-year coppice rotation, when the species reached its maturity phase, distinct differences expressed in the biomass related characteristics in field conditions were demonstrated using 82 ecotypes. Despite the reported low genetic diversity as revealed by molecular analysis of clonal populations collected in America and Europe (Ahmad et al., 2008;Hardion et al., 2012Hardion et al., , 2014Mariani et al., 2010;Pilu et al., 2014), clonal selection can take advantage of inheritable and promising phenotypic differences that have been described in traits such as number of culms, culm diameter and height (Cosentino et al., 2006;Pilu et al., 2014).
Our results highlighted significant differences among ecotypes for all phenotypic traits studied. Concerning phenology, our collection included ecotypes with large

F I G U R E 8 Representation of the principal components analysis showing interrelations among a broad Euro-Mediterranean panel of 82
Arundo donax ecotypes grown in northern Italy. The representative clustering plot for all ecotypes showed their distribution into five major clusters (circle: ecotypes into cluster 1; triangle: ecotypes into cluster 2; plus: ecotypes into cluster 3; cross: ecotypes into cluster 4; diamond: ecotypes into cluster 5) FABBRINI ET AL. | 163 variation in the late flowering stage. Flowering date is considered one of the most important selection criteria when breeding perennial grasses for cellulosic ethanol (Jakob, Zhou, & Paterson, 2009). In addition, flowering phenology and height growth patterns were shown to be significantly associated with functional traits such as relative growth rate, leaf mass per area and hollow ratio in herbaceous grassland species (Sun & Frelich, 2011). Previously Cosentino et al. (2006) scored flowering in 39 ecotypes collected in southern Italy (35 in Sicily region and 4 in Calabria region) and at Ph4. Surprisingly the variability observed in Ph4 reported by Cosentino et al. (2006), which showed a 41 days difference between the two most divergent ecotypes at the end of the second growing season, perfectly mirrors the variability reported here (Figure 3). Perhaps the most interesting is the finding of phenotypic correlations that highlighted a relationship between flowering time and some morphological traits such as H and SLA. Although we compared within-species variation, a statistically significant and positive relationship among these traits has also been observed in previous studies (Bolmgren & Cowan, 2008;Sun & Frelich, 2011), where taller herbaceous grassland and perennial species flowered earlier than shorter ones. As expected, positive correlations between flowering phenology phases and SLA were observed. Both early flowering time and reduced SLA are associated with plants highly adapted to arid environments (Kooyers, 2015). In A. donax, vegetative and reproductive phases often end at the T A B L E 4 Mean cluster values and coefficient of variation (CV, %) of phenological, physiological, morphological, and biomass and growth traits in Arundo donax. The number of ecotypes in each cluster (Clu) is enclosed by curly brackets. CV was calculated as the ratio of the standard deviation among ecotypic values in the cluster to the mean cluster value, reported as a percentage within round bracket. For each trait, significant differences among clusters (F test) are indicated by asterisks. Different letters among clusters are related to LSD post hoc test (P ≤ 0.05). See Table 1  Da (cm) *** 11.89 (6.78)a 12.01 (6.85)a 11.02 (9.77)b 11.17 (6.57)b 10.29 (9.27)c ns: non-significant. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001. same period (from early to late fall; Pilu et al., 2012Pilu et al., , 2013; therefore, late flowering together with high SLA may be indicative of improved growth capacity. Likewise, there were strong positive correlations between flowering time and LAI reflecting greater canopy soil coverage for the ecotypes that had longer growing seasons. Ecotypes with high LAI had longer vegetative periods, indicating that late flowering may be advantageous under continuous mild water deficit as it allows stress acclimatization over time, but not under severe drought as it may imply that they could not finish their phenological cycle. Gradual water reduction in the soil is common in natural ecosystems, and one of the earliest responses to this depletion is by stomata (Bogeat-Triboulot et al., 2007;Chaves, Flexas, & Pinheiro, 2009;Harfouche, Meilan, & Altman, 2014;Plomion et al., 2006). There is also strong evidence that stomatal limitation is preponderant at mild levels of abiotic stresses like drought (reviewed by Cornic & Fresneau, 2002). We have shown a significant variation among ecotypes in mid-day abaxial stomatal resistance and a clear ecotype-dependent response to a period of water shortage as highlighted by the E × T interaction analysis. Indeed, ecotypes presenting the lowest values of Rs july (i.e. using more water through transpiration) had much higher Rs august (as they might have exhausted the available water) than others, and thus showed a higher Rs august /Rs july , for which a negative relationship between both parameters is apparent (Figure 4). However, this response varies among ecotypes. Ideally, the perfect ecotype could be one showing low Rs july (allowing high rate of photosynthesis during the season when there is still some water in the soil) and either having largely increased Rs august (saving water when it is much less available) or keeping low Rs august (meaning that it is still able to uptake water probably due to an efficient root system). It would therefore be interesting in further studies to check drought tolerance of ecotypes positioned at different places along the Rs august /Rs july versus Rs july relationship. The fact that Δ is only marginally related to stomatal resistance suggests that water use efficiency in Arundo ecotypes is more dependent on variability in photosynthetic capacity rather than on stomatal physiology (O'Leary, 1988), for which further studies on the variability and H 2 of photosynthetic capacity would be also of interest for selection purposes. We also observed a negative relationship between stomatal activity and LAI, as already seen in various broadleaved and coniferous forest stands under different climates (Granier, Loustau, & Bréda, 2000).
The selection of suitable ecotypes for sustainable lignocellulosic biomass production has to be based on key characteristics, among which carbon isotope composition could be a relevant trait. Mean Δ observed in this study was lower than the one observed by Nackley, Vogt, and Kim (2014) (Δ =~24‰) in a growth chamber pot experiment. In our study, we observed a difference of 2.8‰ between the two most divergent ecotypes in terms of Δ. In a collection of chickpea (Cicer arietinum), ranges of variation of 3.1‰ and 2.7‰ were observed among genotypes subjected to terminal drought condition (Krishnamurthy et al., 2013). In cocksfoot (Dactylis glomerata), tall fescue (Festuca arundinacea) and reed canary grass (Phalaris arundinacea) grown under water and water-limited conditions, an overall difference of 2.5‰ was observed between the two treatments (Mårtensson et al., 2017). Plants of sunflower (Helianthus annuus) subjected to two drought treatments (i.e. progressive water-stress and unvaried stressful condition), showed a significant difference between treatments, with a Δ of 2.86‰ and 2.71‰ in stressed and no-stressed plants, respectively (Adiredjo et al., 2014). As supported by these papers, the range of variation observed in our study is of high interest and could be the result of the high variability in the ecology of native sites where A. donax ecotypes were collected. The small CV observed in Δ, not exceeding 2.5%, was several times smaller than the CV here estimated for other physiological, morphological and growth related traits. However, caution should be taken when interpreting significances of low Δ variations. Indeed, a supposed low CV value of 2.42% that referred to a general mean value of 20.98‰ would imply a large variation of Δ of about 2.8‰. Such a variation is physiologically relevant, suggesting good levels of adaptation for Δ in our A. donax panel. Correlation analysis found weak and negative relationship between Δ and TotDM; moreover, cluster differentiation highlighted that more productive ecotypes showed lower Δ values. The consistently negative associations between Δ and WUE might suggest that crop yield and Δ might also be negatively related (Condon, Richards, Rebetzke, & Farquhar, 2004). However, large variance in flowering time and reed height, as here observed, could affect the strength and direction of correlation between Δ T A B L E 5 Loadings of the sub-set of informative traits in the first two principal components of the PCA, accounting for 69.23% of the ecotype variability in Arundo donax. See Table 1  and total yield (Condon et al., 2004). Positive correlations among flowering phases and Δ were found in our work, and this result reflects the antagonist behaviours between drought escape (i.e. early flowering) and avoidance (i.e. higher WUE, i.e. lower Δ) strategies (Kooyers, 2015). Finally, negative correlations between Rs and Δ here observed resulted as expected results, probably due to lower CO 2 inside leaf concentration in response to greater stomatal resistance (Farquhar et al., 1989). Canopy architecture affects growth performance of tree species used in SRC plantations (Dillen, Marron, Sabatti, Ceulemans, & Bastien, 2009); therefore, LAI was measured to evaluate the variation in canopy density among the large set of ecotypes analysed here. Other leaf traits affecting canopy structure and light capture have also been characterized in the perennial biomass crop genus Miscanthus . The range of LAI values observed in the current study is consistent with recent findings. The average value of 6.11 m 2 /m 2 is similar to that reported by Cosentino, Scordia, Sanzone, Testa, and Copani (2014) in a trial carried out in southern Italy using a local ecotype at the third year of growth under moderate water supply (6.90 m 2 /m 2 ). Likewise, the maximum value of 8.14 m 2 /m 2 is in agreement with the previous recorded LAI of about 9 m 2 /m 2 for a local ecotype at its third year of growth in northern Italy (Ceotto et al., 2013). Finally, the value of 3.89 m 2 /m 2 measured in the ecotype with the lightest canopy is comparable with the one estimated by Nassi o Di Nasso, Roncucci, Triana, Tozzini, and Bonari (2011) in central Italy, but in a stand at the seventh one-year coppice rotation. Notably, except for Cosentino et al. (2014), the comparison of LAI values was carried out during the same period of time, that is the middle of July, because LAI changes through the season and in A. donax it reaches the maximum value in early August and tends to decline until the end of the season (Ceotto et al., 2013). Leaf dimensions, described by LfL and LfW, were in the range of those reported by Hardion et al. (2012) where A. donax plant material from 16 sites across the Mediterranean basin was investigated. Finally, the observed average value of SLA was consistent with data collected in wild populations of A. donax along moisture gradients in southern USA (Watts & Moore, 2011). Interestingly, SLA of A. donax was suggested as an indicator of a poor adaptation to frequent or prolonged submerged conditions (Mommer, Lenssen, Huber, Visser, & De Kroon, 2006).
Due to its high biomass yields, combined with low input requirements, A. donax has been proposed as an herbaceous model lignocellulosic crop for Mediterranean environments (Angelini et al., 2009;Borin et al., 2013;Nassi o Di Nasso et al., 2011. However, the biomass data from A. donax in Italy differ considerably. In a study by Cosentino et al. (2006) in southern Italy, a mean yield of 11 Mg DM/ha at the end of the first one-year coppice rotation and of 22 Mg DM/ha in the second year were recorded, with a peak of 34.2 Mg DM/ha reached by the most productive clone at the end of the second growing season. In central-southern Italy, Fagnano et al. (2015) obtained a mean yield of about 15 Mg DM/ha in a 9-year trial carried out in an experimental site characterized by severe climatic conditions and a clayey soil texture. Whereas in northern Italy, a 12-year field study conducted by Angelini et al. (2009) highlighted an average yield of 38 Mg DM ha −1 yr −1 with a peak of about 50 Mg DM/ha at the end of the second and third years of growth under non-limiting conditions of water and nutrient availability. Similarly, in another experimental site in northern Italy, with optimal water and N inputs, Borin et al. (2013) showed a high biomass yield of about 100 Mg DM/ha at the second and third year of growth, being the highest A. donax biomass yield ever recorded in Italy, although the biomass yields from small-scale experimental plots may be higher than those achievable at commercial scale. This result was supported by our data; the highest yielding ecotypes in our collection were more comparable with the yields of 150 Mg DM/ha reported in wild fields in California and India (Giessow, Casanova, Lecler, MacArthur, & Fleming, 2011;Sharma, Kushwaha, & Gopal, 1998;Spencer, Liow, Chan, Ksander, & Getsinger, 2006). However, our yield estimates are likely to be overestimates, because our results were obtained from single plots in a common garden collection and because we were unable to avoid the bias of increased yields often seen in plants growing at the edges of plots. This is also not unexpected as biomass yields in commercial-scale fields might well be lower than those obtained in small research plots (Ping et al., 2014;Verhulst, Cox, & Govaerts, 2013). For our purpose, the homogeneous plot design of the common garden experiment was successful. Comparing biomass components, the variability of H measured in the current study was perfectly comparable to the range for A. donax of between approximately 340 and 600 cm, including panicle (20-70 cm), described by Hardion et al. (2012). The mean H value reported here was higher than those observed in southern and central Italy of 339 cm and 286 cm in trials at the second and seventh year, respectively (Cosentino et al., 2006;Nassi o Di Nasso et al., 2011). As compared to the experiment carried out in Sicily (Cosentino et al., 2006) Correlation analyses in A. donax indicated that high dry biomass yield was largely due to plant apical growth and stem density (Angelini et al., 2009;Cosentino et al., 2006). The low but statistically significant correlations observed between biomass production and various growth traits (i.e. StDM, StN, StNd and H) highlighted biomass yield as a complex trait made up of multiple simple traits (Clifton-Brown, Robson, & Allison, 2008).
H 2 assessments in this study showed that biomass traits had moderate to high H 2 estimates, whereas morphology and ecophysiology trait H 2 estimates ranged from low to high. Low H 2 estimates were observed primarily in physiology traits (Rs, Δ), which is most likely related to some phenotypic plasticity and reflecting the environment in which the trait was sampled. Notably, the wide change of H 2 estimates for Rs july (H 2 = 0.28) and Rs august (H 2 = 0.62) is to be expected because stomatal resistance is an instantaneously measured physiological trait related to leaf gas exchange (Ackerly et al., 2000). Other studies in A. donax have also observed moderate H 2 in biomass and morphology traits (Cosentino et al., 2006;Pilu et al., 2014). Individual traits with higher H 2 (morphology, biomass) tended to be intercorrelated. Our H 2 results for biomass traits (0.32 ≤ H 2 ≤ 0.39) are slightly lower than other perennial grasses such as switchgrass (0.45 ≤ H 2 ≤ 0.50; Rose, Das, & Taliaferro, 2008;Jiang et al., 2014) and miscanthus (0.34 ≤ H 2 ≤ 0.55; Slavov et al., 2013;Gifford, Chae, Swaminathan, Moose, & Juvik, 2015).
Our multivariate analysis demonstrated that biomass and growth traits: StDM, StN, StNd, H, Db and Dm, as well as LAI, were the most indicative of dry biomass production in our panel of 82 ecotypes. These growth-related traits allowed us to identify similar groups and to select promising individual ecotypes. Cluster analysis highlighted clusters 2, 3 and 5 as the most interesting for biomass yield potential and that the high biomass producing ecotypes adopted different growth strategies. For example, ecotypes in Clu.3 deployed heavy and large leaves in a lighter canopy. Ecotypes in the lower side of the biplot were characterized by a higher number of reeds and higher biomass production compared with individuals localized in the upper part (Figure 8). The two most divergent groups, Clu.2 and Clu.5, contained ecotypes that exhibited significant differences in reed diameter and number and, as a consequence, in LAI. High-biomass ecotypes with large and tall stems organized into less-dense and light clumps were grouped in Clu.2. Conversely ecotypes grouped in Clu.5 had small reeds forming dense clumps. Clusters 2, 3 and 5 correspond to 56% of our collection (46 ecotypes); thus an extensive multi-year multi-location examination of phenotypic variation of these 46 ecotypes among many different types of traits could provide considerable insight into key traits favouring the production of biomass yield in different cropping environments. Finally, the most promising ecotypes belong to these three clusters that showed the best combinations of biomass determinants such as H, StN, StNd and LAI. Although Δ is not a significant biomass determinant, Clusters 2, 3 and 5 showed lower Δ than the other two, and this underlines its importance. In some growing areas in semi-arid environments, biomass potential may not be achieved due to water shortage, and hence water use efficiency may be a more relevant trait. These traits should be taken into consideration for both clonal selection and systems and synthetic biology approaches. While selection for these traits is warranted, targeted genetic improvement of A. donax bioenergy feedstock in these critical morphological and biomass traits is of great importance. However, to maximize the efficiency of biomass conversion and reduce biomass recalcitrance in A. donax, future breeding efforts should also target cell wall properties optimized for saccharification and fermentation (Scordia, Cosentino, Lee, & Jeffries, 2011. Creating a successful bioenergy program from lignocellulosic feedstocks is a major challenge and will require significant research effort to maximize plant biomass production. Our study provides an essential stepping stone towards this goal. The results presented here demonstrate high variation in phenotypic traits among a large collection of A. donax ecotypes spanning most of the species' natural range. The assembled panel also represents a suitable resource to screen for improved bioconversion efficiency using a similar approach to improve the already favourable comprehensive utilization index for A. donax as discussed by Kou et al. (2017). The study also highlights the availability of different ideotypes (described by Donald (1968) as a biological model which is expected to perform or behave in a predictable manner within a defined environment) of A. donax that achieve high yield through different combinations of traits and it will be interesting to see if these ideotypes differ in their ability to generate high yield across a range of environments. Extensive assessments of traits underscored the consistency observed in phenotypic traits and to our knowledge present the most comprehensive assessments of biomass traits in A. donax to date. Efforts are under way to exploit genomic tools in A. donax to identify loci responsible for the observed trait variation.