Fodder Yield, Quality and Growth of Chia (Salvia hispanica L.) as Affected by Sowing Density and Top-Dressing Nitrogen Fertilization

Abstract

Chia (Salvia hispanica L.) seeds are considered a superfood, but research on vegetative biomass uses is scarce. We conducted a two-year trial under the hypothesis that short-day flowering black chia would provide quality forage at high latitudes and tested sowing density (D1 = 125, D2 = 25, D3 = 8 and D4 = 4 plants m⁻²) and top-dressing N fertilization (N = 0, N = 20 kg ha⁻¹) at three stages (EV = early vegetative, LV = late vegetative and EF = early flowering). Forage yield reached 9.0 and 5.64 t ha^-1 of total and 2.07 and 1.56 t ha⁻¹ of leaf dry mass at EF in 2013 and 2014, respectively. High plant density corresponded to higher biomass but to lower plant relative growth rate, stem diameter, branching, proportion of leaves and lipid and crude protein content. Crude protein declined from 18% at EV to about 8% at EF and was lower in D1 than in D4. Lipids ranged between 33.83 g kg⁻¹ (D4, EV) and 17.34 g kg⁻¹ (D1, EF) (p < 0.005). N topdressing affected forage quality but not yield. Alpha-linolenic acid was the most abundant fatty acid (FA) (608 g kg⁻¹ FA at EV). The long vegetative growth of short-day flowering chia in southern Europe is favorable to fodder production, and management should be optimized by seeking balance between growth rate and stand density while optimizing vegetative stage growth for the highest forage quality.

1. Introduction

Widening the range of crops by using neglected or under-utilized species has many advantages [1]. Exploiting their multifunctionality may be an important tool to attain global sustainability and to raise the nutritional and nutraceutical potential of crops. Salvia hispanica L., commonly known as chia, is an annual macrothermal plant from the genus Salvia of the Lamiaceae family. Due to the high proportion of alpha-linolenic (ALA) fatty acid (FA) (~60% of total FA) in the seeds, chia is considered one of the richest botanical sources of omega-3 [2]. The seeds are also a source of proteins [3] and antioxidants [4]. The potential of chia for multiple functions and ecological services is slowly emerging from research. Leaf metabolites [4] indicate possible roles in health, medicine and the control of plant pathogens. The ecological role of chia myxodiaspores has been highlighted by work in References [5,6,7]; the mucilage capsule extruded by germinating achenes has been shown to create a microenvironment relevant for the water relations of the small seedling, row-soil stability and plant-microorganism relations [6]. Chia mucilage has been shown to interact with soil chemicals [5], to stabilize soil aggregates [6,7] and to retain soil at the surface of hydrated chia achenes against flowing water [6].

Chia is becoming increasingly popular in the U.S.A., Australia and Europe, where it has a strong presence in new product launches and innovation; novel uses may be proposed on the basis of properties relevant for improving the rheological and nutraceutical properties of gluten-free pasta [8] and for amending soil [6]. Chia has also been proposed as a superfood; a large body of literature encourages the use of chia in animal nutrition [9]. Chia seeds have been used as an omega-3 boosting supplement in poultry science [3,10,11,12]. They have been shown to improve both milk [13] and meat [14] lipid profile in ruminants. Research has also been conducted on other species, and this shows a higher proportion of omega 3 in the meat of rabbit [15,16], pork [17], fish [18] and edible insects [19]. Improving the lipid profile of animal products is a very important target for the industry since many products are perceived by consumers as unhealthy due to the association between consumption and increased incidence of cardiovascular disease, diabetes and some types of cancer [20,21,22]. In ruminants, nutrition is a fast, natural and low-cost way to sharply improve FA profile of milk [23,24]. Lactating goats supplemented with chia seeds were shown to produce milk with a lower content of saturated fatty acid and a higher proportion of health-promoting conjugated linoleic acid, while the atherogenicity index decreased by as much as 25% [13]. Until now, the use of chia in animal nutrition has been limited to the consumption of seeds. However, potential for forage use was shown [25] since high forage quality of vegetative parts was found, especially at early vegetative stages, when a high digestibility was coupled with a proportion of alpha-linolenic acid slightly above 60%. Chia leaves contain several nutraceuticals, among which are flavonoids such as acetyl vitexin, acetyl orientin, apigenin and luteolin glycosides [4]. The crop has been experimented with in Europe, in both the north [25] and south of Italy [26], in Greece [27] and in Germany [28]. At European latitudes, however, seed production is currently limited by thermal and photoperiod requirements; domesticated chia varieties are short-day flowering since floral initiation requires a photoperiod threshold of 11.8 h and a base development temperature of 10 °C [29]. In Europe, low temperatures prevent sowing before late spring, and lower initiation only begins at the onset of autumn, when low temperatures slow down grain ripening and plants are killed by frost before seed maturation [26]. Novel long-day flowering genotypes have been developed [30,31] but are not freely available to European farmers. Switching the cultivation from seed to fodder production could be an interesting option at European latitudes where cooler and longer growing seasons are favorable to herbage biomass accumulation. Despite the great interest in chia cultivation, literature on the effect of agronomic management on chia yield and quality is scarce, and even less is available on biomass production. Even when oriented towards seed production, research has been mainly focused on assessing the impact of location, sowing dates and seed source [2,32,33,34,35,36]. Very few works addressed the effect of sowing density alone [29] or in combination with fertilization strategies [27,28,37]. Agronomic management protocols still need to be established, especially for novel uses such as forage production and growing chia outside the area of origin. Our aim is to test the effect of sowing density in combination with top-dressing N fertilization on a short-day flowering commercial black chia grown in southern Italy. In a preliminary report aimed at investigating seed yield, biomass and leaf area [26], we found that seed yield was low because of late-flowering at a high latitude.

In this work, we used a two-year agronomic trial to test the hypothesis that late-flowering allows vegetative biomass of good forage value and omega-3 boosting potential to accumulate in Mediterranean Europe and to assess the effect of density and nitrogen top-dressing on growth dynamics and forage quality. To this end, we performed growth analysis and investigated the architecture, forage yield and gross properties from the early vegetative stage to the onset of flowering.

2. Materials and Methods

2.1. Experimental Setup

The experiment was conducted in 2013 and 2014 at Masserie Saraceno (Atella, PZ, Southern Italy, Lat. N 40°51′37.59″, Long. E 15°38′49.43″) on a Luvi-vertic Phaeozem [38] loam soil with the following characteristics: sand (50–2000 μm) 43.6%, silt (2–50 μm) 34.2% and clay < 2 μm) 22.1%. The soil was amended in June 2013 with 25 t ha⁻¹ of a solid fraction of biogas digested materials with the following characteristics: dry matter 8.5%, carbon 20.4 kg t⁻¹, nitrogen (N) 2.8 kg t⁻¹, ammonium (N-NH₄) 0.6 kg t⁻¹, P₂O₅ 1.4 kg t⁻¹ and K₂O 2.5 kg t⁻¹. After amendment, the soil chemical parameters were: pH 6.8, N 1.9 g kg⁻¹, P₂O₅ 50.3 mg kg⁻¹ and K₂O 1430 mg kg⁻¹. The site has an average annual rainfall of 678 mm, concentrated mainly during October to May. The precipitation and temperature during the experiment are reported in Figure 1.

Black chia (Salvia hispanica L.) (Eichenhain, Hofgeismar, Germany) was sown on 21 June 2013 and 26 June 2014 to a 1 cm depth and grown with a non-limiting water supply. Every 2 days, drip irrigation provided 100% of ET0 corresponding to the evaporative demand of the atmosphere [39] measured with a TE-ETG atmometer (Tecnoel, Rome, Italy). Drip irrigation amounted to 711.3 m³ ha⁻¹ in 2013 and 1193.0 m³ ha⁻¹ in 2014. Soil tillage consisted of ploughing at 35 cm and accurate seedbed preparation with chain and rotary harrows, given the small seed size (1000 seed weight = 1.26 g). The crop was hand-weeded. Top-dress nitrogen fertilization and sowing density were tested with the following treatments:

• (i). N top-dress fertilization at two levels, T0 and T20, respectively corresponding to 0 and 20 kg N ha⁻¹; the T20 treatment was fertilized with ammonium nitrate (13.5% NH₄⁺ and 13.5% NO₃⁻) applied 7–8 weeks after sowing;

• (ii). sowing density at four levels (D1 = 125, D2 = 25, D3 = 8 and D4 = 4 plants m⁻²) with a distance of 0.8 m interrow.

The treatments were tested in a randomized complete block split-plot design with three replications, where N top-dressing was assigned to the main plots and sowing density to the split plots. The split plots were 5 × 8 m².

2.2. Plant Measurements

Samples were taken at the early vegetative (EV, 7 weeks after sowing), late vegetative (LV, 10 weeks after sowing), early flowering (EF, 10% of stems had flowers) stages by harvesting plants of 100 cm on the row at randomly chosen positions per triplicate plot.

The fresh weight of the leaves and stems was determined within 5 min of cutting. The samples were divided in two subsamples. One was frozen at −20 °C at the field premises. The other was dried at 70 °C until it was a constant weight to determine dry biomass. The shoot biomass was calculated as the sum of the leaf and stem mass. The stem diameter, plant height and number of primary branches were measured at all growth stages (EV, LV and EF) in 2013, while, in 2014, measurements were made at LV and EF. The number of leaves was measured in 2013.

The ratio of leaf to stem biomass was calculated as:

Leaf to steam ratio (LSR) = leaf dry biomass/shoot dry biomass (g g⁻¹))

The relative growth rate was calculated for plants of T0 for two periods: from EV to LV and from LV to EF. It was computed as:

Relative growth rate (RGR) = (mean ln shoot biomass t2—mean ln shoot biomass t1)/(t2 − t1) g g⁻¹ d⁻¹ plant⁻¹ [40].

2.3. Chemical Analyses

Chemical analyses were performed on the D1 and D4 plant densities. The samples were dried, ground to pass through a 1-mm screen and analyzed for the following determinations: the total N content by the Kjeldahl method, the ash content by ignition to 550 °C and ether extract by the Soxhlet method were determined as described in Reference [41]. Acid detergent fiber (2ADF) and neutral detergent fiber (NDF) were determined as described by Reference [42] expressed exclusive of residual ash. Lignin was determined by solubilization of cellulose with sulphuric acid, as described in Reference [43]. Lipid extraction was performed on freeze-dried samples according to Reference [44]. Fatty acids (FA) were analyzed as their methyl esters (FAME). The analysis was carried out by gas chromatography using a Varian Star 3400 CX GC (Varian-Agilent, Milan, Italy) equipped with a SLB^®-IL111 Capillary GC Column (100 m × 0.25 mm × 0.20 μm) (Sigma-Aldrich, Milan, Italy). The separation was carried out at 90/240 °C with helium as the carrier gas and using a flame ionization detector (FID) at 300 °C. FAMEs were identified by the comparison of retention times with FAME standard mixture under the same conditions (Supelco 37 Component FAME Mix analytical standard, Sigma-Aldrich, Milan, Italy).

2.4. Statistical Analysis

The experiment was designed as a split-plot randomized complete block design replicated three times, with the top-dressing N levels (0–20 Kg ha⁻¹) randomly assigned to the main plots within each of the three blocks and the sowing density (D1:D4) assigned to the sub-plot. The effects of plant density and nitrogen fertilization were assessed by a mixed effect model estimated by restricted maximum likelihood [45]. We fitted a mixed effect model, testing the fixed effect of plant density, crop stage, fertilization and year and of their interaction, as well as a nested random effect (Nitrogen level, block or growth, stage or year) [45]. As a consequence of including random terms, a correlation structure between observations belonging to the same group is incorporated in the model; this is called induced correlation [46]. The analysis of variance was followed by Tukey’s post-hoc means comparison. Regression analysis was conducted to analyze the relationship between plant traits. Statistical analysis was performed by R software (R Foundation for Statistical Computing, Vienna, Austria) [47].

3. Results

3.1. Crop Growth and Forage Yield

The results of the ANOVA for biomass are reported in Table 1, and the results of the post-hoc mean comparison for total dry mass (TBDW), leaf dry mass (LBDW) and stem dry mass (SBDW) are shown in Figure 2.

We found that TBDW and the individual components LBDM and SBDW, as well as their fresh biomass counterparts (TBFW, LBFW and SBFW, respectively), were significantly affected by sowing density and by growth stage. There was also a significant interaction between year and growth stage, while fertilization had no significant effect on forage yield.

Chia forage yield increased until the beginning of flowering (EF), when dry biomass ranged between 9.27 (D3 N20) and 7.57 (D4 N20) t ha⁻¹ in 2013 and between 7 (D1 N0) and 4.6 (D3 N20) t ha⁻¹ in 2014. About 44 days after sowing (EV), dry biomass ranged between 0.70 t ha⁻¹ (D3 20N) and 1.80 t ha⁻¹ (D1 0N) in 2013 and was much less in 2014, ranging from 0.19 t ha⁻¹ (D1 N0) and 1.05 t ha⁻¹ (D4 N20) respectively. In 2014, the biomass values were lower than in 2013 at both the EV and LV stages, and the difference was reduced but still significant at the EF stage. Leaf-to-stem ratio (LSR) was the only parameter affected by a four-way interaction that included the effect of N fertilization. LSR from the EV to EF stages decreased on average from 1.08 to 0.32 in 2013 and from 1.29 to 0.35 in 2014. Only at EV were differences between sowing densities significant (p < 0.005), with values ranging between 0.70 in D1 and 1.44 in D4 in 2013. In 2014, LSR was 0.77 in D1 and 1.71 in D4.

Top-dressed D4 plants showed a significantly higher LSR at LV (1.64 vs 0.86; p < 0.005). LSR decreased with total biomass, following a power relationship (y = 0.8805x −0.479, R² = 0.90, p < 0.05).

Different letters over bars indicate significant differences among means by Tukey’s multiple comparison test at p < 0.05).

The relative growth rate (RGR, Figure 3) increased with decreasing sowing density, although differences were not statistically significant between D1 and D2 and, in the period from EV to LV, between D3 and D4. The values of RGR were also significantly higher in the period from LV to EF in D4 only. In 2013, RGR was higher than in 2014 but only for the D4 density (2.80 vs. 1.74 g g⁻¹ day⁻¹ plant⁻¹ in 2013 and 2014, respectively).

3.2. Chia Architecture

Table 2 reports ANOVA results for plant architectural parameters. Sowing density and growth stage alone or in interaction significantly affected architecture, while fertilization did not.

The results of post hoc mean comparison are depicted in Figure 4. In both 2013 and 2014, stem diameter and the number of primary branches (Figure 4a,b,e,f) were affected by plant density (p-value < 0.0001), but the values of the two highest densities (D1 and D2) could not be significantly discriminated, while their differences from the lower densities (D3 and D4) became more marked as the plant cycle proceeded. In fact, while at low densities, values increased throughout the plant cycle; in the D1 treatment, there were no significant differences between growth stages.

Crop height (Figure 4c,g) increased with phenological stage in both years, and, in 2013, there was a significant interaction (p < 0.005) between growth stage and sowing density since, at EF, only height was significantly lower in D1 than in D3 and D4, while no effect of density was detected in other cases (Figure 4c). Height was linearly related to total plant fresh (y = 14.89x − 252.15; R2 = 0.8; p < 0.05) and dry (y = 2.25x − 47.71; R² = 0.8; p < 0.05) mass, as well as to the fresh mass of leaves (y = 2.78x + 8.69; R2 = 0.79; p < 0.05). The total number of leaves per plant was measured only in 2013 (Figure 4g), when a significant interaction was found between sowing density and growth stage. The number of leaves increased from D1 to D4 in LV and more so in EF, but the gap between densities was amplified with growth stage.

3.3. Forage Quality

Table 3 reports the ANOVA results for forage gross properties, and the results of post-hoc mean comparisons are depicted in Figure 5, Figure 6 and Figure 7.

The quality markedly changed during growth, declining from EV to EF, and was affected by sowing density, the year alone or interactions between factors.

Figure 5 depicts the main effects of the experimental treatments on forage quality, while Figure 6 and Figure 7 show the effects of the interactions with sowing density and N-fertilization, respectively. Crude protein content was significantly higher in D4 (14.18%) compared to D1 (12.28%) (Figure 5a; p < 0.005). Differences between years were significant for crude protein and lignin content; namely, in 2014, there was a higher protein content while lignin decreased (Figure 5b,c; p < 0.05). Ash content was only influenced by the growth stage declining from EV to EF (Figure 5d; p < 0.05). For lipids, lignin and NDF, there was a significant interaction between sowing density and growth stage; lipids decreased with increasing density (Figure 6a; p < 0.005), while fiber content increased (Figure 6b,c, p < 0.005).

Fertilization as a main effect was significant for crude protein, lignin and acid detergent fiber (ADF) (Table 3). There were significant interactions between fertilization and growth stage for crude protein content (Figure 7a; p < 0.005) and for lignin (Figure 7b; p < 0.005). Crude protein declined from 18% EV to about 8% at EF, and this loss was higher in unfertilized plots at both the LV and EF stages. ADF was affected by a significant interaction between fertilization, growth stage and sowing density (Figure 7c; p < 0.005).

The results of ANOVA for the forage fatty acids (FA) profile are reported in Table 4. The FA composition of chia lipids is dominated by polyunsaturated FA, of which ALA was by far the most abundant, with an average of 54%, followed by linoleic acid (11%) (average across growth stages, years, fertilizers and density levels).

The most abundant saturated fatty acid (SFA) was palmitic acid, with an average proportion of 117.70 g kg⁻¹ FA, while stearic acid constituted only 30.41 g kg⁻¹ FA. Monounsaturated oleic acid was, on average, 28.83 g kg⁻¹ FA. In many cases, individual FA proportions changed during growth and between years. The results of the post-hoc comparison are depicted in Figure 8 and Figure 9.

Saturated FA palmitic (C16) and stearic (C18) increased during growth, with a significant interaction between year and growth stage (Figure 8a,b; p < 0.05). Oleic acid increased during growth, but differences between years were not significant (Figure 8c; p > 0.05). Only the linoleic acid proportion increased slightly but significantly with sowing density (Figure 9c; p < 0.05). Linoleic acid was also affected by a significant interaction between growth stage and year (Figure 9a; p < 0.05), while, in 2013, the linoleic acid proportion increased steadily during growth. In 2014, a slight decrease at LV was recorded.

Finally, the concentration of ALA shows a trend opposite to the other FAs; it decreased during growth from a maximum at EV of 608 g kg⁻¹ FA to a minimum value of 403 g kg⁻¹ FA at EF (average across sowing density, fertilization and year) (Figure 9b; p < 0.05). ALA was slightly lower in 2014 (522.4 g kg⁻¹ FA) compared to 2013 (558.9 g kg⁻¹ FA) (Figure 9d; p < 0.05). In no cases were FA concentrations influenced by fertilization. The ratio of omega-6 to omega-3 increased with growth but was always <0.3 and did not differ between densities (Figure 9e; p < 0.05).

4. Discussion

4.1. Crop Growth and Forage Yield

Literature on the effects of agronomic practices on chia biomass production is quite scarce [29], especially regarding the effect of sowing density and fertilization [27,37]. This is partially due to the fact that chia is mainly grown as a seed crop, and only recently has interest in fodder production been raised [4,9,25]. The chia shoot biomass values found in our experiment are higher than the values obtained in northwestern Argentina [34] across different locations and seeding dates (between 0.49 and 4.67 t ha⁻¹) and as affected by geographic location [35] (between 1.12 and 4.18 t ha⁻¹). This can be ascribed to the long vegetative growth we recorded using a short-day flowering genotype at a high latitude. Our values are in agreement with those reported for the Mediterranean environment of southern Europe, where dry shoot biomass values ranging between 4.48 and 15.36 t ha⁻¹, as a function of sowing rate and fertilization, are reported in Greece [27]; these are higher than the values obtained for the Mediterranean environment of Chile [29], with values between 2.05 and 5.15 t ha⁻¹ for early and late sowing dates, respectively. A direct comparison of our data on the dry matter yield of chia at early flowering with yields of forage crops is difficult due to differences in harvest stage and to biomass values coming from multiple cuts for some species. Nevertheless, as an example of irrigated annual crops in southern Italy, Reference [48] reports an experiment in which different species of clover (Trifolium alexandrinum L., T. incarnatum L., crimson; T. resupinatum L. and T. squarrosum L.) were compared, and their average yield was 4.36 t ha⁻¹ of dry matter taken when 10–15% of the stems were flowering. By contrast, Reference [49] report yields ranging from 5.5 t ha⁻¹ for Italian ryegrass (Lolium multiflorum Lam.), to 32.66 t ha⁻¹ for silage maize (Zea mays L.), harvested at around 2/3 milkline in kernel and therefore at a physiological stage corresponding to a high potential biomass accumulation.

In our trial, sowing density proved to be a relevant factor affecting fodder yield. Experimental evidence on the role of sowing density in chia growth and development [27,37] reports that plant density significantly affects crop growth. Compared to our data, much less biomass was reported [37] for the tropical environment of Ghana (0.23 t ha⁻¹ at a density of 4 plants m⁻²). Differences in sampling stage (end of grain filling vs. early flowering in our experiment) and photoperiod are probably the main reasons why plant growth was much lower in their case and why the number of leaves was less than half of ours. In Greece [27], it was shown that a reduction of the inter-row space from 0.60 to 0.40 m corresponds to a three-fold biomass increase (from 4.48 to 14.75 t ha⁻¹).

High biomass values can be achieved by increasing sowing density, but our data demonstrate that the biomass yield gain we obtained in D1 was counterbalanced by a forage quality loss. Density effects on branching pattern and leaf yield and proportion (e.g., parameters related to forage quality) need to be considered. Concerning the role of fertilization, we found no significant effect on forage biomass; this is consistent with References [4,27]. In the literature, data on chia nitrogen fertilization are inconclusive so far; doses of nitrogen fertilizer from 0 to 125 kg ha⁻¹ were tested in a pot experiment [48], and vegetative growth at 30 and 60 days after sowing was found to be highest at the maximum dose. The same paper cites technical papers on chia management in which recommendations on nitrogen rates span quite a large range (15 to 100 kg ha⁻¹ of N). Other experiments show little or no effect on plant biometrics from N [26], and authors [26,50,51] argue that fertilization can actually reduce plant density by increasing lodging rates and can hence reduce total biomass on a surface basis. The same comment was made in a paper reporting that heavy fertilization coupled with long vegetative cycles in the Sichuan Basin, China [52] leads to lodging rates as high as 98% by stimulating vegetative growth. Nitrogen can enhance lodging also through a lower lignification of plant tissues. The effect of the water regime on lodging is controversial [52].

Our data on RGR (Figure 3), span one order of magnitude, showing the large effect of sowing density on plant size and growth dynamics. The values at high densities (D1 and D2) are of the same order of those found in a study in Chile [53,54] in which maximum values of 0.15 g g⁻¹d⁻¹ plant⁻¹ are reported for chia grown at 50 plants m⁻², which corresponds to an intermediate density between our D1 and D2. In that research, an early increase and a later decrease in RGR were found, and this is commented on by assuming that the plant’s basal leaves contribute less to photosynthesis as time passes. Comparison with our data can only be made in a wide sense since we calculated RGR for two long periods, both of which were during vegetative growth, but the explanation in Reference [54] is consistent with our finding of a significant interaction between sowing density and plant stage for RGR. The values of RGR go from a slight non-significant reduction with time in dense plant standings (D1), where it is likely that basal leaves become shaded during growth, to an increase in the second period at lower densities that becomes significant only at 4 plants m⁻². At this density, the wide spacing between plants is likely to allow basal leaves to keep receiving solar radiation and contributing to photosynthesis later than in dense plants, and they may do so even more with time due to the very strong increase in the number of leaves that we recorded in D4 (Figure 4).

Based on field research [55], spring wheat (T. aestivum L.) is defined a fast-growing species (0.085 g g⁻¹ day⁻¹) compared to the slow growing Aegilops comosa Sm. in Sibth. and Sm. var. comosa (0.068 g g⁻¹ day⁻¹). For faba bean (Vicia faba L.), pea (Pisum sativum L.) and lentil (Lens culinaris Medik), the minimum relative growth ranged between 0.10 and 0.20, 0.26 or 0.25 g g⁻¹ day⁻¹, respectively [56]. Chia values are therefore in range with those of the major broadleaf crops, but they are shown to be greatly affected by stand density. The ecological implications of high RGR are intuitive: fast growth results in a competitive advantage due to the rapid accumulation of space; a fast growth cycle, for instance, is essential for ruderals [57,58]. On the opposite side, slow growers have been classified as “stress tolerators” [58]. As reviewed in Reference [59], a low RGR is beneficial in case of limited resources since a low demand may imply a slower resource depletion rate. From an agronomic perspective, the growth rate has implications for competitive ability, which might be critical for stand establishment. This is relevant for weed management in chia stands, where chemical weed control is still under evaluation [60], and a compromise between growth rate and stand density in relation to competition with weeds needs to be found with further research.

Plant morphological plasticity as a response to sowing density is relatively low in species with a strong apical dominance (e.g., sunflower), while other crops, such as soybean or wheat [61], can modify their space-filling capacity by changing the branching pattern and tiller density. In our trial, chia showed a strong response to sowing density, with a three-fold variation of the number of primary branches from D1 to D3 and D4 and with a more than two-fold variation in stem thickening. The effect of plant density on the number of primary branches was mostly observed at the EF stage for D2, D3 and D4, while in the high-density stand (D1) branching pattern remained almost unaltered during the growth cycle, showing that this is a limiting density for canopy expansion. With the exception of the D1 group, plant morphology changed during growth in our experiment: branching patterns and most of the stem thickening occurred during the last third of the growth cycle (LV to EF). The chia plants reached a remarkable size (1.6 m height), higher than the maximum height of 1 m reported for the south area of Santa Fe province, Argentina (33°14′ S, 61°2′ W) [62] and higher than that reported in Ghana [37]. A greater height is reached in chia when, due to the photoperiod, the reproductive stage is delayed; Reference [29] argues that, in Mediterranean sites, during the great majority of the growth cycle, photo-assimilates are used for biomass production until autumn’s shorter days trigger the partition of dry matter to the seeds. In Reference [29], height is negatively associated to seed yield and to the harvest index but positively related to shoot biomass and leaf production. Crop height was not significantly affected by plant density, but low heights were measured in D1 at EF in 2013, and this height reduction coupled with the consistently low number of primary branches suggests that D1 might be a threshold density above which plant growth is reduced. Under high plant density, leaf photosynthesis can be decreased due to shading; Reference [63] observed the effect of stand density on the photosynthetic performances of soybean and showed that photosynthesis was linearly related to light interception during the first part of the growth cycle but that, after canopy closure, light interception did not decline at a rate proportional to the loss of leaf area, indicating that abscission of leaves was not involved in light interception. All three of the morphological parameters measured in our trial had higher values compared to those reported in Reference [37] at 4 plants m⁻² (corresponding to our D4 treatment) in the tropical environment of Ghana. The number of branches reported by Yeboah et al. (2014) is only 22 at 4 plants m⁻²; the maximum height was 97 cm, and stem diameter was below 0.6 cm. This difference can be possibly related to our longer growth cycle. In their case, physiological maturity was reached 66 days after planting, while, in our case, the onset of flowering was recorded around 100 days after sowing. Our primary branch number was similar to that reported in China by Reference [52], which measured chia morphology as affected by population densities and by sowing dates. Unlike us, they found no effect of stand density on chia morphology. In their trial, the growth cycle lasted up to 273 days, being about 2.6 times longer than ours, and their plants were about 1.6 times taller than ours, with a maximum height of 250 cm. According to those authors, these very long growth cycles caused overgrowth and very high lodging rates (up to 98.52%). They report that the upper stem restored an upright growth habit after falling down, but this could have caused alteration in morphological development; for instance, branching was significantly lower in an early-sown stand that suffered more from lodging. The authors also emphasized that early spring sowing, in addition to lodging, was also more affected by summer high temperatures and drought, and they concluded that all densities considered in their study were probably too high, causing competition for light and nutrients and high lodging rates.

4.2. Forage Quality

The assessment of chia fodder quality is important in view of a possible expansion of the crop at European latitudes, where seed production is currently limited by the photoperiod [26,29,64]. We found that the chia whole-plant is rich in beneficial FA and crude protein, especially at the early vegetative stage. As for all forage species, in our data, herbage quality declines during growth. Substantial quality loss of chia with age has been also reported in northern Italy [25] and Greece [27]. Our crude protein values at LV and EF are below those reported by Reference [27] at 61, 91 and 161 DAS (days after sowing), though ADF is slightly lower at all growth stages. Compared to temperate summer legumes and grasses grown under non-limiting water and nutrient supplies [65], chia ADF and NDF at the late vegetative stage were higher than those of alfalfa, similar to those of cow-pea and slightly lower than those of many summer grasses. Chia crude protein concentration at the late vegetative stage was lower than that of summer grasses and legumes (>20% DM), and only at the early stages of growth (EV) did it reach the threshold values of 16.0–17.9% DM suggested by Reference [66] as adequate to meet rumen microbial N needs. The highest quality forage for chia can be obtained by harvesting plants at the EV stage. This does not necessarily mean that the crop needs to be terminated early, since regrowth after herbage cuts may correspond to setting the ontogenetic stage back. In view of forage production, therefore, the effect of multiple cuttings and regrowth needs to be investigated, as for common sage (Salvia officinalis L.) for which research indicates that high yields of good-quality fresh herb can be obtained by harvesting the plants twice per growing season at a cutting height between 10 and 15 cm above the soil level [67]. A good yield of vegetative biomass may be exploited also for alternative multifunctional uses, such the extraction of nutraceuticals, due to the high content of anti-oxidants [4], or the natural control of plant pathogens, due to the antimicrobial properties of leaf essential oils and other compounds [68]. In addition, dual-purpose uses might be envisaged, wherein an early cut for high-quality forage would be coupled with a later harvest of seeds. This multiple use implies that whole-plant behaviors and adaptation to different environments are considered and tested. For instance, a second harvest of seeds may be pursued with common short-day flowering chia types at low latitudes. In European environments, though, long-day flowering types would be needed in order to obtain early flowering and early seed maturation so that grain filling is not impaired by low fall temperatures [26,64].

We found that plant density affects forage quality. Our data indicate that fodder quality improves in low density stands, corresponding to a higher leafiness. Crude protein increases while ADF, NDF and lignin decrease with decreasing sowing density. For the lipid fraction, the effect of sowing density changes with growth stage, but the quality gain of D4 compared to D1 is always remarkable: +44%, +28% and +27% at the EV, LF and EF stages, respectively. However, Reference [27] found no significant effect of plant density on quality parameters. This can be possibly attributed to a narrower range of variability in their study (0.60 vs 0.40 cm row spacing). Reducing density to improve quality obviously implies a yield penalty. Therefore, while it is true that concentrations of desirable compounds in shoots increase at low sowing densities, the yield of such compounds on a field surface area basis may not. For instance, in our data, crude protein contents increase by 15% on average by growth stage from D1 to D4 but shoot biomass yield is reduced by up to 21.5%. As a result, the yield of proteins that may be harvested from D4 is 9% lower than that obtained in D1 (0.48 vs. 0.53 t ha⁻¹ respectively).

Another finding of our study is chia responsiveness to fertilization. We found that that even a low dose of N topdressing increased crude protein at LV, where the value for N20 is 13% higher than that for N0, and more so at EF, where the value for N20 is 30% higher than that of N0. This is in agreement with Reference [27], which showed that manure increased chia herbage crude protein by 13%. Fertilization significantly reduced lignin at LV and ADF in D1 at LV and EF. This confirms results for other plant species [69] for which nitrogen fertilizers are found to interfere with the deposition of structural components of plant cells like lignin and other fiber fractions. Such changes are also linked to a decrease in mechanical strength and an increase in lodging [69]. Fertilization strategies therefore need to take into account the lodging risk [26,52], especially for dense standings where the shading of the lower parts of stems also plays a role, as confirmed by the interaction of N topdressing and sowing density on ADF in our data. In chia, the risk of lodging is likely higher for short-day flowering varieties at high latitudes, where they reach remarkable sizes due to long vegetative growth.

Our data are in general agreement with those published for northern Italy [29]. A maximum of above 60% ALA is the main FA, followed by linoleic acid. The percentage of ALA we found is in line with values reported in the literature as reviewed in Reference [9] and slightly lower than the values reported in Reference [29]. The omega-3 content of chia vegetative parts is remarkable compared to the average values reported for some common forage grasses and legumes [70] and higher than values reported for other omega-3-rich species such as birdsfoot trefoil (Lotus corniculatus L.) and salad burnet (Sanguisorba minor Scop.) [71]. The ALA concentration is also higher than values reported for quinoa (Chenopodium quinoa Willdenow) at early vegetative stages (460 g kg⁻¹ of total FA) [72]. Currently, one of the main goals of the dairy and meat industry is to improve omega-3 content. Many studies demonstrated that even small amounts of chia seeds in animal diets can dramatically improve the omega-3 content of eggs [11,19,73] and of poultry and lamb meat [12,74]. Nevertheless, data on lactating goats fed with chia seeds [13] show that the FA profile of milk is improved but that ALA undergoes a substantial biohydrogenation in the rumen. The effect of including chia forage in the diet of ruminants remains to be ascertained, but it is possible that leaf antioxidants would play a role in protecting polyunsaturated fatty acids (PUFA).

Our data show a low LA/ALA ratio. This is reported to be the major determinant of the total omega-6/omega-3 ratio in food and feed and to be widely used to closely track it [75]. Both the omega-6/omega-3 ratio and its proxy LA/ALA are too high in animal and human diets in western countries, and this is thought to correspond to an increased risk of cancer and cardio-vascular problems to the point that forage with a low omega-6/omega-3 ratio is considered an important contribution to healthy animal food [75]. In our data, the LA/ALA ratio ranges between 0.17 and 0.26 from EV to late EF, in agreement with values reported by Reference [25] at similar growth stages. Such values are lower than those reported by Reference [71] for caraway (Carum carvi L.) (0.89) and alfalfa (Medicago saltiva L.) (0.59) and lower than values found in the omega-3 rich birdsfoot trefoil (0.33). This might be important in the feed industry.

5. Conclusions

Short-day flowering chia grown in southern Italy is suitable for the production of vegetative parts with good potential for high quality forage uses due to good biomass yield, especially at high plant density, and to a higher content of omega 3 and a lower omega 6/omega 3 ratio compared to many omega 3-boosting species.

However, a forage yield gain of 21.5% at 125 plants m⁻² compared to 4 plants m⁻² corresponded to a lower quality due to a severe reduction of the number of primary branches, a reduction of the number of leaves of one order of magnitude and an almost halved proportion of leaf to stem biomass. A lower crude protein and lipid content and a higher fiber content were also found at D1. Forage quality was quite high at early vegetative stages and declined during growth. The effect of N topdressing was not seen on biomass accumulation, but it positively affected quality mostly in interaction with growth stage. We can conclude that chia should be harvested before the beginning of the reproductive stage in order to obtain a high-quality forage and that future research might focus on testing multiple cuts as a strategy to ensure a good seasonal biomass production while improving forage quality by setting back the ontogenetic stage. This may prevent short-day flowering chia varieties to reach a large plant size even when grown during long-days, which, in turn, may help control the related lodging risk. Our findings may also be useful for the design of further research directions on the possible dual-purpose uses of chia with an early cut for high-quality forage followed by a later harvest of seeds. This would require taking crop regrowth, yield dynamics and photoperiod sensitivity into account. A strong dependance of RGR on sowing density also indicates that agronomic research is needed, and specifically optimizing the competitive abilities of chia towards weeds may be pursued through seeking compromises between growth rate and stand density in the fast occupation of space and acquisition of resources.

Animal nutrition trials with forage chia still need to be conducted, but data on fodder yield and quality indicate that commercially available short-day flowering varieties can be grown in Mediterranean Europe to obtain a potential forage rich in omega-3. Perspective forage exploitation would extend the current limits of chia cultivation at northern latitudes and open alleys to the development of innovative leaf-based products for the emerging market of omega-3-rich feed and food products.

Author Contributions

Conceptualization: M.A., R.R. and R.B.; methodology: R.R., R.B. and G.B.; data curation: R.B., R.R. and R.L.; writing: R.R. and M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We wish to thank “Masserie Saraceno” and “Eurogel” by Donato De Carlo for invaluable assistance with the logistics of the experiment and Carmela Bochicchio for assistance with plant biometrics.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Monthly temperature minima (dashed line) and maxima (solid line) and precipitation (bars) during 2013 and 2014.

Figure 2. Dry biomass as affected by experimental treatments; (a–c) interaction between growth stage and year. Bars represent means averaged over density and fertilization. EV = early vegetative; LV = late vegetative: EF = early flowering. (d–f) main effects of sowing density. Bars represent means averaged over fertilization, stage and year; D1 = 125 plants m−2; D2 = 25 plants m–2; D3 = 8 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means across all combinations of growth stages and years (a–c) or across all levels of density (d–f) by Tukey’s multiple comparison test ((p < 0.005).

Figure 3. Relative growth rate (RGR) for the vegetative growth (EV-LV) and the transition to reproductive stage (LV-EF) as affected by growth stage and sowing density. Bars represent means averaged over years. EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D2 = 25 plants m−2; D3 = 8 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means across all combinations of sowing densities and growth stage transitions by Tukey’s multiple comparison test ((p < 0.005).

Figure 4. Plant architecture as affected by experimental treatments as main effects or in interaction in 2013 (a–d) and 2014 (e–g). Bars indicate means averaged over fertilization (a–f) and over density and fertilization (g). EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D2 = 25 plants m−2; D3 = 8 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means across all combinations of sowing densities and growth stages (a–f) or among means across growth stages (g) by Tukey’s multiple comparison test ((p < 0.005). Different letters over bars indicate significant differences among means by Tukey’s multiple comparison test (p < 0.005).

Figure 5. Chia gross forage quality characteristics affected by the main effects of management experimental factors. Bars indicate means averaged over fertilization, growth stage and year (a); over fertilization, growth stage and density (b,c); or over fertilization, density and year (d). EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means by Tukey’s multiple comparison test (p < 0.005). Lipid content was linearly related to crude protein content (y = 0.1113x + 0.8748; R² = 0.56; p < 0.005), and neutral detergent fiber (NDF) increased with growth stage, with slightly higher values in D1 plots at the EV stage (Figure 6b; p < 0.05), but no significant differences between plant densities were found thereafter.

Figure 6. Interaction between sowing density and growth stage on gross forage quality: (a) lipids, (b) NDF and (c) lignin. Bars indicate means averaged over fertilization and year. EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means across all combinations of growth stage and sowing densities D1 and D4 by Tukey’s multiple comparison test (p < 0.005).

Figure 7. Interaction between nitrogen topdressing treatment and other management experimental factors. Effects on gross forage quality: (a) crude protein, (b) lignin and (c) acid detergent fiber (ADF). Bars indicate means averaged over sowing density and year (a,b) or over years (c). EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D4 = 4 plants m−2. Different letters over bars indicate significant differences among means across all combinations of growth stage and nitrogen treatments (a,b) or all combinations of growth stage, nitrogen treatments and densities D1 and D4 by Tukey’s multiple comparison test (p < 0.005).

Figure 8. Saturated (a,b) and monounsaturated (c) fatty acid profile as affected by the growth stage × year interaction. Bars indicate means averaged over sowing density and nitrogen treatment. EV = early vegetative; LV = late vegetative; EF = early flowering. Different letters over bars indicate significant differences among means across all combinations of growth stage and year by Tukey’s multiple comparison test (p < 0.005).

Figure 9. Polyunsaturated fatty acid profile (a–d) and LA/ALA ratio (e) as affected by experimental treatments as main effects or in interactions. Bars indicate means averaged over sowing density and nitrogen treatment (a,e); over density, year and nitrogen treatment (b); growth stage, year and nitrogen treatment (c); or density, growth stage and nitrogen treatment (d). EV = early vegetative; LV = late vegetative; EF = early flowering. D1 = 125 plants m−2; D4 = 4 plants m−2. LA = Linoleic Acid; ALA = alpha-linolenic acid. Different letters over bars indicate significant differences among means across all combinations of growth stage and year (a,e) or across all levels of growth stage (b), density (c) and year (d) by Tukey’s multiple comparison test ((p < 0.005).

Table 1

Results of ANOVA for plant growth as affected by fertilization, sowing density, growth stage and year and their interaction. TBDW = total biomass dry weight, t ha⁻¹; LBDW = leaf biomass dry weight; SBDW = stem biomass dry weight; TBFW = total fresh biomass dry weight; LBFW = leaf fresh biomass dry weight; SBFW = stem fresh biomass dry weight; LSR = leaf to stem ratio; Num DF = numerator degrees of freedom; Den DF = denominator degrees of freedom.

			p Value
Source	Num DF	DenDF	TBDW	LBDW	SBDW	TBFW	LBFW	SBFW	LSR
(Intercept)	1	72	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Fertilization	1	2	0.5691	0.4533	0.8873	0.6612	0.5363	0.6881	<0.0001
Sowing Density	3	72	<0.0001	0.0004	<0.0001	<0.0001	0.0017	<0.0001	<0.0001
Growth stage	2	8	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Year	1	12	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	0.5526
Fertilization × Sowing Density	3	72	0.396	0.5549	0.1625	0.5454	0.6568	0.4997	<0.0001
Fertilization × Growth stage	2	8	0.926	0.4744	0.9639	0.6295	0.605	0.6426	<0.0001
Sowing Density × Growth stage	6	72	0.8639	0.0956	0.9025	0.7967	0.0768	0.5947	<0.0001
Fertilization × Year	1	12	0.9388	0.4576	0.3229	0.179	0.5622	0.1293	0.0248
Sowing Density × Year	3	72	0.591	0.7113	0.6358	0.1588	0.4773	0.1161	0.1174
Growth stage × Year	2	12	<0.0001	0.0001	<0.0001	<0.0001	<0.0001	<0.0001	0.8082
Fertilization × Sowing Density × Growth stage	6	72	0.7207	0.928	0.3845	0.4465	0.8024	0.3411	0.0011
Fertilization × Sowing Density × Year	3	72	0.2604	0.4949	0.4418	0.3741	0.5947	0.3375	0.0002
Fertilization × Growth stage × Year	2	12	0.7798	0.6735	0.8332	0.9913	0.3935	0.8934	0.0002
Sowing Density × Growth stage × Year	6	72	0.5927	0.9735	0.6482	0.7945	0.9992	0.5641	0.0644
Fertilization × Sowing Density × Growth stage × Year	6	72	0.4847	0.4096	0.5035	0.5073	0.5357	0.5219	0.0009

Table 2

Results of ANOVA for plant architecture as affected by fertilization, sowing density, growth stage and their interaction in 2013 (top) and 2014 (bottom). Num DF = numerator degrees of freedom; Den DF = denominator degrees of freedom.

				p-Value
	NumDF	DenDF	Stem Diameter (cm)	Plant Height (cm)	Number of Primary Branches	Number of Leaves Per Plant
2013
(Intercept)	1	36	<0.0001	<0.0001	<0.0001	<0.0001
Fertilization	1	2	0.4527	0.4228	0.3113	0.3218
Sowing density	3	36	<0.0001	0.2889	<0.0001	<0.0001
Growth stage	2	8	<0.0001	<0.0001	<0.0001	<0.0001
Fertlization × Sowing density	3	36	0.2544	0.6341	0.309	0.1526
Fertilization × Growth stage	2	8	0.7126	0.8113	0.2539	0.6596
Sowing density × Growth stage	6	36	<0.0001	0.0466	<0.0001	<0.0001
Fertilization × Sowing density × Growth stage	6	36	0.8918	0.6608	0.4969	0.1968
2014
(Intercept)	1	24	<0.0001	<0.0001	<0.0001
Fertilization	1	2	0.2724	0.0596	0.2749
Sowing density	3	24	<0.0001	0.1175	<0.0001
Growth stage	1	4	0.0003	<0.0001	0.0111
Fertlization × Sowing density	3	24	0.5208	0.1184	0.5641
Fertilization × Growth stage	1	4	0.4111	0.0834	0.5074
Sowing density × Growth stage	3	24	0.0003	0.4	0.0165
Fertilization × Sowing density × Growth stage	3	24	0.2144	0.3896	0.8758

Table 3

Results of ANOVA for forage gross properties as affected by fertilization, sowing density, growth stage, year and their interaction. Num DF = numerator degrees of freedom; Den DF = denominator degrees of freedom.

			p-Value
Source	NumDF	DenDF	Crude Protein %	Ether Extract %	Ashes %	ADF %	NDF%	Lignin %
(Intercept)	1	23	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Fertilization	1	23	<0.0001	0.2496	0.5276	0.0001	0.2515	0.0029
Sowing Density	1	23	<0.0001	<0.0001	0.2926	<0.0001	0.0179	<0.0001
Growth stage	2	23	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Year	1	23	<0.0001	0.5927	0.9879	0.2147	0.2628	0.0483
Fertilization × Sowing density	1	23	0.0726	0.2124	0.3109	0.0043	0.1164	0.6627
Fertilization × Growth stage	2	23	<0.0001	0.8405	0.7173	0.0051	0.5945	0.023
Sowing density × Growth stage	2	23	0.2891	0.018	0.4815	<0.0001	<0.0001	<0.0001
Fertilization × Year	1	23	0.5703	0.972	0.945	0.9332	0.9337	0.818
Sowing density × Year	1	23	0.419	0.9241	0.9022	0.3903	0.9016	0.4281
Growth stage × Year	2	23	0.9403	0.9998	0.99	0.8832	0.9865	0.8288
Fertilization × Sowing density × Growth stage	2	23	0.1524	0.18	0.6569	0.0374	0.1236	0.1924
Fertilization × Sowing density × Year	1	23	0.3299	0.9161	0.8598	0.6378	0.8718	0.7456
Fertilization × Growth stage × Year	2	23	0.6455	0.9958	0.9885	0.9453	0.6155	0.731
Sowing density × Growth stage × Year	2	23	0.1535	0.9844	0.9977	0.6563	0.6922	0.9074
Fertilization × Sowing density × Growth stage × Year	1	23	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001

Table 4

Results of ANOVA for whole-plant individual fatty acid (FA) proportion (g kg⁻¹ FA) as affected by fertilization, sowing density, growth stage, year and their interaction. Num DF = numerator degrees of freedom; Den DF = denominator degrees of freedom. LA = Linoleic Acid; ALA = alpha-linolenic acid.

			p Value
Source	NumDF	DenDF	Palmitic	Stearic	Oleic	LA	ALA	LA/ALA
(Intercept)	1	23	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Fertilization	1	23	0.9181	0.5195	0.3949	0.4514	0.7868	0.6155
Sowing Density	1	23	0.3095	0.2869	0.8152	0.0055	0.0833	0.4245
Growth stage	2	23	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
Year	1	23	<0.0001	0.0154	<0.0001	<0.0001	<0.0001	<0.0001
Fertilization × Sowing density	1	23	0.4437	0.9141	0.3162	0.9245	0.4084	0.659
Fertilization × Growth stage	2	23	0.8475	0.5878	0.36	0.1928	0.8489	0.5828
Sowing density × Growth stage	2	23	0.1249	0.4909	0.0158	0.0824	0.5052	0.6137
Fertilization × Year	1	23	0.8372	0.9141	0.8152	0.5712	0.9668	0.6747
Sowing density × Year	1	23	1	0.7466	0.9379	1	0.6037	0.8035
Growth stage × Year	2	23	0.0001	0.0099	0.0031	<0.0001	0.1707	<0.0001
Fertilization x Sowing density × Growth stage	2	23	0.3882	0.678	0.4745	0.7158	0.9162	0.8408
Fertilization × Sowing density × Year	1	23	0.8775	1	0.8152	0.5093	0.8514	0.5556
Fertilization × Growth stage × Year	2	23	0.638	0.9622	1	0.921	0.9821	0.9072
Sowing density × Growth stage × Year	2	23	0.6201	0.9397	0.9058	0.4744	0.8384	0.6985
Fertilization × Sowing density × Growth stage × Year	2	23	0.8738	0.9911	0.8014	0.8412	0.9418	0.9087