Stevia rebaudiana (Bertoni) as a Multifunctional and Sustainable Crop for the Mediterranean Climate

Abstract

Stevia rebaudiana (Bertoni) is a promising medicinal and aromatic plant for Mediterranean agroecosystems given its positive agronomic attributes and interesting quality features. It has both food and pharmaceutical applications, since its leaves contain sweet-tasting steviol glycosides (SVglys) and bioactive compounds, such as phenolics, flavonoids, and vitamins. We evaluated the agronomic and qualitative performances of nine stevia genotypes cultivated, in open field conditions, for two consecutive years under the Mediterranean climate of central Italy. Growth, biomass production, and accumulation of bioactive compounds (SVglys, polyphenols, and their related antioxidant activities) were evaluated, considering the effect of harvest time and crop age (first and second year of cultivation). The results showed high variability among genotypes in terms of both morpho-productive and phytochemical characteristics. In general, greater leaf dry yields, polyphenol accumulation, and antioxidant activities were found in the second year of cultivation, harvesting the plants in full vegetative growth. On the other hand, total SVglys leaf content reached the highest values in the first year when plants were at the beginning of the reproductive phase. On the other hand, although the SVglys profile (Rubusoside, Dulcoside A, Stevioside, Rebaudioside A, C, D, E, and M) remained stable over harvest times, it differed significantly depending on the crop age and genotype. Our findings provide useful information on the influence of crop age and harvest time in defining quanti-qualitative traits in stevia, with PL, SL, BR5, and SW30 being the best performing genotypes and thus suitable for breeding programs. Our study highlighted that stevia, in the tested environment, represents a promising semi-perennial crop which offers new solutions in terms of cropping system diversification and marketing opportunities.

1. Introduction

In many Mediterranean areas, high evapotranspiration rates, increased precipitation variability, and intense summer drought are the main environmental constraints to agricultural management and crop production. Under such conditions, sustainable agricultural practices and the introduction of new crops that are able to diversify cropping systems and to mitigate climate changes are very promising strategies for addressing cropping system sustainability.

Compared to continuous cropping systems, crop diversification involves better utilization of land resources, lower risks from pests and diseases, and greater yield stability. In this context, medicinal and aromatic plants (MAPs) can be included in low-input productive systems: they can enhance the multifunctionality traits of the agricultural sector, produce safe final products, respect the environment, and promote rural areas. Various MAPs find optimal growth conditions in Mediterranean environments and, among these, perennial and semi-perennial species can contribute to agro-ecosystem services.

The cultivation of perennial crops reduces soil erosion, minimizes nutrient leaching, sequesters more C in soils, protects water resources, creates a better pest tolerance, and provides a continuous habitat for wildlife [1,2,3]. In addition, perennial species require less use of farm equipment than annual crops as well as fewer fertilizers and herbicides [4]. Among perennial MAPs, Stevia rebaudiana (Bertoni) is grown successfully in a wide range of agro-ecological environments, from semi-humid, subtropical to temperate zones. Thanks to its extreme versatility, it can be grown as a pluriannual crop (three to five years) in temperate to warm climates and, as an annual crop, in colder regions [5,6,7]. Unlike other traditional sweetener crops, such as sugar beet and sugar cane, stevia cultivation reduces agronomic inputs, including nutrients, water, and energy [8].

Besides these agronomic benefits, stevia has an exceptional phytocomplex composition that can be exploited in several applications, from the pharmaceutical to cosmetic and nutraceutical industry. Metabolic disorders such as type-II diabetes and obesity, associated with an excessive sugar consumption, are becoming more prevalent and stevia is the perfect sugar substitute in foods and drinks, given that its leaves contain non-calorie high-intensive sweetener compounds, namely steviol glycosides (SVglys) [9,10,11]. Its leaves also contain a complex mixture of triterpenoids, sterols, essential oils, phenols, and flavonoids, with functional and health-promoting properties [12,13,14,15,16,17]. Thanks to its positive agronomic and phytochemical characteristics, stevia could thus offer new solutions in terms of cropping system diversification and marketing opportunities.

Stevia cultivation is widespread all over the world and it has been introduced as a commercial crop in several countries [8,18]. Experimental cultivations are increasing in Europe aimed at producing higher-performing and yielding stevia genotypes [19,20,21,22,23]. Access to genotypes that can adapt to different environmental conditions is important for the selection of those characterized by both high levels of secondary metabolites and biomass yield.

To improve the competitiveness of stevia production in the Mediterranean region, it is important to identify higher-performing genotypes in terms of yield and quality, resource use efficiency, and resistance/tolerance to a wide range of biotic/abiotic stress combinations. Therefore, the aim of this study was to determine the best performing genotypes suitable for developing site-specific recommendations for stevia cultivation under the Mediterranean climate of central Italy. Nine stevia genotypes and two harvest times were compared for two consecutive years, exploring differences in plant growth traits (growth cycle, biometric characteristics, leaf yield) and in desired compounds (SVglys content and profile, total phenol and flavonoid accumulation, and antioxidant activity).

2. Materials and Methods

2.1. Experimental Design and Plant Materials

A field-plot trial was carried out during two growing seasons (2018 and 2019) at the experimental Centre of the Department of Agriculture, Food and Environment (DAFE, University of Pisa), located in San Piero a Grado, Pisa, central Italy (43°40′ N; 10°19′ E, 5 m above sea level). A completely randomized block experimental design, with four replications, was adopted, with the plot as the experimental unit and selected plants within the plot as the observational unit. The plot size was 3 m × 1.8 m (width × length), and each plot consisted of six rows of six plants (36 plant plot⁻¹), with a plant density of 6.67 plants m⁻² by adopting an inter-row and intra-row spacing of 0.5 × 0.3 m. Seedbed preparation included moldboard plowing (30 cm depth), disk harrowing, and use of a cultivator. Pre-planting phosphorus and potassium fertilizations were performed at a rate of 100 kg ha⁻¹ of P₂O₅ by triple superphosphate and 80 kg ha⁻¹ of K₂O by potassium sulfate. Nitrogen (as ammonium nitrate) was applied at a rate of 40 kg N ha⁻¹ after transplanting. Nine genotypes of different origins, belonging to DAFE’s germplasm collection were used for the field-plot trial (

Table 1. Origin of S. rebaudiana genotypes.

Genotype ID	Origin
PL	Israel
BR16	Brazil
RG	Italy
SL	Israel
NU	Italy
CO	Israel
BR5	Brazil
SW30	Italy
BR1	Brazil

). Plants for each genotype were first reproduced by micropropagation and then multiplied by stem cuttings to avoid plant genetic variability in terms of both morphological and phytochemical characteristics. Subsequently, they were kept under controlled conditions at the DAFE greenhouse until open field transplanting between mid-May and early June 2018.

Weather parameters (daily minimum, T_min, maximum, T_max, and mean temperatures, T_mean, and cumulative rainfall) were recorded by an automated weather station nearby the experimental site, from the beginning of vegetative plant development (June) to the full plant flowering (October) in each year of cultivation. Throughout the experiment, total rainfalls were 164.2 and 233.2 mm, in 2018 and 2019, respectively, with a mean average temperature of 19.0 °C and 16.4 °C in 2018 and 2019, respectively. Mean average T_max and T_min temperatures did not differ notably between 2018 and 2019 (

Table 2. Monthly cumulative (mm) rainfall and mean temperatures (T_mean, T_max, and T_min) throughout the entire experimental period (May 2018–September 2019).

Year	Month	Tmax (°C)	Tmin (°C)	Tmean (°C)	Rainfall (mm)
2018	May	21.7	15.0	18.4	87.4
	June	25.6	17.6	21.6	5.0
	July	28.0	20.3	24.1	49.6
	August	29.3	20.4	24.8	25.2
	September	26.1	17.2	21.6	19.4
	October	22.7	14.1	18.4	65.0
	November	16.8	9.9	13.4	111.6
	December	13.3	5.5	9.4	51.0
2019	January	11.4	1.9	6.6	41.6
	February	13.9	4.1	9.0	58.6
	March	15.6	6.6	11.1	3.0
	April	17.6	9.0	13.3	112.2
	May	18.6	11.6	15.1	87.6
	June	26.3	17.8	22.1	2.8
	July	28.6	20.2	24.4	88.4
	August	28.9	20.2	24.6	3.6
	September	25.5	17.3	21.4	72.6

).

Before beginning the experiment, soil samples were collected at 0 to 30 cm depth in order to evaluate the physical and chemical characteristics of the soil (

Table 3. Physical and chemical characteristics of the soil at the experimental site.

Characteristics	Values
Clay (<0.002 mm, %)	14.5
Silt (0.05–0.002 mm, %)	38.9
Sand (2–0.05 mm, %)	46.5
pH (H2O 1:2.5 soil:water suspension; McLean method)	7.6
N tot (Kjeldahl method, g kg−1)	1.4
S.O. (Walkley–Black method, g kg−1)	3.1
Available phosphorus (Olsen method, mg kg−1)	9.9
Exchangeable potassium (Thomas method, mg kg−1)	211.5
CE (mS cm−1)	0.4
CSC (Method BaCl2, pH 8.1, meq 100 g−1)	20.1

). The soil was loam with a sub-alkaline reaction, a good content of total nitrogen, organic matter, and exchangeable potassium, and with a low level of available phosphorus and salinity.

2.2. Crop Sampling and Agronomic Measurements

Phenological observations, from the beginning of vegetative plant development until full plant flowering, were performed in the ratooning crop, starting from the regrowth of new shoots, per each genotype and in each plot. For every observation, six plants per genotype were randomly selected in the field, with 54 total plants observed. The length of both the vegetative and reproductive phases was evaluated using accumulated thermal time (°C/day) and the accumulated growing degree days (GDDs) were calculated daily according to the equation (1) presented by McMaster and Wilhelm [24],

GDD = Σ [(T_max + T_min)/2 − T_b]

(1)

where T_max is the daily maximum air temperature, T_min is the daily minimum air temperature, and T_b is the base temperature, for which 10 °C was used [22].

For each growing season, two destructive samplings (namely HT1 and HT2) were carried out at the following phenological stages (according to the BBCH scale by Le Bihan et al. [25]): HT1–Stage 48 = about 80% of final leaf biomass is developed; and HT2–Stage 55 = 50% of apex leaves are differentiated and present inflorescence, but flower buds are still closed.

HT1 corresponded to the full vegetative growth which was reached between 27 July–2 August in the 1st year and between 15–17 July in the 2nd year, depending on the genotype. HT2 corresponded to the beginning of flowering, occurred in early September in the 1st year (10–13 September 2018) and between the end of August and the beginning of September in the 2nd year (from 20 August to 2–3 September 2019). Samplings were manually performed, collecting six plants per genotype (1 sample = 1 plant). After each harvest, plant height (cm), basal stem diameter (mm), branching (n_. stem/plant), fresh and dry weight of leaves and stems (g plant⁻¹), harvest index (HI), and specific leaf weight (LSW) were measured. Leaves of each sample were air-dried in a ventilated oven from 30 to 40 °C until constant weight. Dry leaves were ground to a fine powder by a laboratory mill (Grindomix GM 200, Retsch, Pedrengo (BG), Italy) and stored until the subsequent analyses. The harvest index (HI), which represents the plant’s efficiency at producing leaves, was calculated as the ratio between leaf dry yield (g plant⁻¹) and total aboveground biomass yield (leaves plus stems). Specific leaf weight (SLW) was calculated as the ratio between dry leaf weight and leaf area (mg cm⁻²). The leaf area was measured by collecting and subjecting to a color scan, two apical leaves, two middle leaves, and two basal leaves. The leaf images were then analyzed with ImageJ (Fiji Particle Analysis plug-in).

2.3. Phytochemical Analysis

2.3.1. Chemicals

Common Steviol Glycosides Standards Kit (Rubusoside, Dulcoside A, Stevioside, Rebaudiosides A, C, D, E and M) was purchased from Chromadex (LGC Standards S.r.L., Milan, Italy). Ferrous sulphate, 2,2-diphenyl-1-picrylhydrazyl (DPPH), gallic acid monohydrate (3,4,5-trihydroxybenzoic acid), 2,4,6-tri(2-pyridyl)-triazine (TPTZ), trizma acetate, Folin-Ciocalteu reagent, Trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid), sodium carbonate, and ferric chloride were obtained from Sigma-Aldrich Chemical Co. (Milan, Italy). All chemicals used in the present study, including solvents, were of analytical grade.

2.3.2. Sample Extraction

A total of 0.1 g of leaf powder per sample was extracted with 10 mL of 70% (v/v) EtOH, and sonicated for 30 min at 60 °C. At the end of the sonication, the extracts were centrifuged (3500 rpm for 10 min) and filtered with a syringe filter (Ø 0.45 μm) to remove any suspended material. The extracts obtained were stored at 4 °C until subsequent analyses.

2.3.3. Steviol Glycosides Determination

The extraction procedure and the determination of steviol glycosides (SVglys) were carried out following Zimmerman et al. [26]. Steviol glycosides were analyzed using a Jasco PU980 HPLC system (JASCO Benelux B.V., Utrecht, Netherlands) coupled with a UV-visible wavelength detector. A hydrophilic column (Luna HILIC 200A, Phenomenex Inc., Torrance, CA, USA), 5 μm, 250 mm × 4.6 mm (Phenomenex Inc., Torrance, CA, USA) in combination with the corresponding pre-column (4 × 3.0 mm) was used. UV detection was carried out at 205 nm at room temperature with a flow rate of 0.68 mL/min and a run time of 20 min. Separation was achieved in acetonitrile/water (80:20) as isocratic mobile phase at pH 3.6 regulated with acetic acid. Chromatograms were acquired online, and data were collected using a Jasco interface (Hercules 2000 Interface Chromatography). Steviol glycosides were quantified using authentic standards, through calibration curves (0.005–1.00 g L⁻¹), obtained from standard mixtures containing Rubusoside (Rub), Dulcoside A (Dulc A), Stevioside (Stev), and Rebaudiosides A, C, D, E, and M (Reb A, C, D, E, and M).

2.3.4. Analysis of Total Phenols and Flavonoids

Total phenols were determined using the Folin-Ciocalteu method according to Dewanto et al. [27] and expressed as gallic acid equivalents (mg GA g⁻¹ dry leaf). This method involved the reduction of Folin-Ciocalteu reagent by phenolic compounds, with a blue complex formation determined at 765 nm by UV-Vis spectrophotometer (Varian Cary 1E, Palo Alto, CA, USA). Total flavonoids were determined using the aluminum trichloride method according to Jia et al. [28] and expressed as catechin equivalents (mg CE g⁻¹ dry leaf). The flavonoids-aluminum reaction created a pink complex formation measured at 510 nm using a UV-Vis spectrophotometer (Varian Cary 1E, Palo Alto, CA, USA).

2.3.5. Ferric Reducing Antioxidant Power (FRAP) Assay and Free Radical-Scavenging Assay

The determination of ferric reducing antioxidant power and the free radical-scavenging activity (FRAP and DPPH assay) of stevia leaf extracts followed Tavarini et al. [29]. The FRAP method is based on the ability of the antioxidant compounds to reduce Fe³⁺ to Fe²⁺ which, in the presence of TPTZ (2,4,6-tris(2-pyridyl)-s-triazine) led to a blue complex formation (Fe²⁺-TPTZ), measured at 593 nm using a UV-Vis spectrophotometer (Varian Cary 1E, Palo Alto, CA, USA). The DPPH assay is based on the reducing activity of the antioxidant molecules against the 1,1-diphenyl-2-picryl-hydrazil (DPPH) radical which was characterized by a purple red color. The extent of the disappearance of DPPH is directly proportional to the amount of antioxidant present in the reaction measured at 517 nm using a UV-Vis spectrophotometer (Varian Cary 1E, Palo Alto, CA, USA). Total antioxidant activity and free radical-scavenging capacity were expressed as trolox equivalents (mmol TE g⁻¹ dry leaf).

2.4. Statistical Analyses

All data were subjected to analysis of variance (ANOVA) using GraphPad Prism v. 8.0.2 (GraphPad Software, Inc., La Jolla, CA, USA). A three-way ANOVA analysis was conducted to assess the effect of genotype (G), crop age (CA), harvest time (HT), and their reciprocal interactions, on the agronomic characteristics and on total phenols and flavonoids, total SVglys and antioxidant activities (FRAP and DPPH). Means were separated on the basis of the least significant difference (LSD) only when the ANOVA F test showed significance at 0.05 or 0.01 probability level.

Hierarchical cluster (HC) and principal component (PC) analyses were performed on (i) the total content of phenols, flavonoids, and steviol glycosides, and (ii) on the individual steviol glycoside (Dulcoside A; Rebaudiosides A, C, D, E and M; Rubusoside, Stevioside) concentrations with JMP^® Pro 13.2.1 (SAS Institute Inc., Cary, NC, USA). As unsupervised methods, the groups of samples obtained with hierarchical cluster analysis (HCA) and PCA can be observed even when there are no reference samples that can be used as a training set to establish the model. For both observation groups, the hierarchical cluster analysis (HCA) was conducted on the normalized average values, with Ward’s algorithm, using Euclidean distances as a measure of similarity among the samples. In addition, principal component analyses (PCA) were carried out in order to reduce the dimensionality of the multivariate data of the matrix, whilst preserving most of the variance [30]. For the PCA of the total content of phenols, flavonoids, and steviol glycosides, a 35 × 3 (35 samples, 3 metabolite contents, 105 total data) dimensional matrix was used. The score plot obtained was defined by a PC1 and a PC2 covering 78.3 and 21.3% of the variance, respectively, for a total explained variance of 99.6%. For the PCA of the total content of individual steviol glycoside (Dulcoside A; Rebaudiosides A, C, E, M, and D; Rubusoside, Stevioside) concentrations, a 35 × 8 (35 samples, 8 steviol glycosides, 280 total data) dimensional matrix was used. The score plot obtained was defined by a PC1 and a PC2 covering 68.2 and 19.3% of the variance, respectively, for a total explained variance of 87.5%.

3. Results

3.1. Growth, Biometric, and Productive Measurements

Figure 1 shows the accumulation of thermal time, expressed as growing degrees days (GDDs), required by each genotype to pass from the vegetative phase to reproductive one in the second year after transplanting. Significant differences among genotypes were observed: RG, SL, CO, BR1, and NU (1414.4 GDD °C d⁻¹) accumulated fewer GDDs to develop the first flower buds than the others (BR16 and SW30 = 1779.15 °C d⁻¹, PL, and BR5 = 1575.2 °C d⁻¹).

Genotype (G), crop age (CA), harvest time (HT) and their reciprocal interactions significantly affected biometric and productive traits, except for (i) branching, which did not vary depending on CA, (ii) basal stem diameter, in relation to G × HT interaction and (iii) plant height and SLW which were not affected by CA×HT (

Table 4. Results of three-factorial analysis of variance (ANOVA) for agronomic characteristics of stevia genotypes tested in the open-field experimental trial.

Factor	Factor Level	Height (cm)	Basal Stem Diameter (mm)	Branching (n. Stem Plant−1)	HI (%)	SLW (mg cm−2)	Leaf Dry Yield (g Plant−1)
Main Effects
Genotype (G)	PL	47.9 ± 3.1 de	7.5 ± 0.4 bc	18.0 ± 2.3 b	65.1 ± 2.5 cd	11.3 ± 0.3 ab	43.7 ± 3.6 bc
	BR16	48.5 ± 2.4 d	5.1 ± 0.2 e	27.0 ± 3.3 a	52.8 ± 2.7 g	9.4 ± 0.5 e	26.9 ± 3.2 f
	RG	55.3 ± 4.4 ab	7.6 ± 0.3 b	24.0 ± 2.8 a	56.0 ± 3.4 f	9.9 ± 0.4 de	36.9 ± 4.3 cde
	SL	44.5 ± 3.7 e	6.8 ± 0.4 cd	15.0 ± 1.7 b	73.5 ± 2.8 a	11.5 ± 0.4 a	29.1 ± 4.3 ef
	NU	58.2 ± 4.7 a	7.7 ± 0.3 b	19.0 ± 1.7 b	55.5 ± 3.4 f	9.5 ± 0.3 e	55.7 ± 5.9 a
	CO	56.6 ± 3.9 a	6.9 ± 0.4 cd	15.0 ± 1.3 b	64.0 ± 2.4 d	10.2 ± 0.3 cd	41.5 ± 3.0 bcd
	BR5	52.9 ± 2.7 bc	7.2 ± 0.4 bcd	17.0 ± 1.3 b	67.1 ± 2.8 b	10.6 ± 0.4 c	44.3 ± 3.8 bc
	SW30	47.3 ± 2.2 de	8.8 ± 0.2 a	15.0 ± 1.6 b	62.2 ± 2.7 e	10.8 ± 0.5 bc	33.9 ± 2.1 def
	BR1	50.0 ± 4.5 cd	6.6 ± 0.2 d	15.0 ± 1.3 b	65.8 ± 2.2 bc	10.6 ± 0.5 c	48.5 ± 5.4 ab
Crop age (CA)	2018	48.8 ± 1.6 b	7.5 ± 0.2 a	19.0 ± 1.1 a	56.0 ± 0.8 b	9.3 ± 0.1 b	35.4 ± 2.9 b
	2019	53.9 ± 1.8 a	6.8 ± 0.2 b	18.0 ± 1.1 a	67.9 ± 1.6 a	11.4 ± 0.2 a	44.6 ± 1.9 a
Harvest time (HT)	July	39.2 ± 0.8 b	6.6 ± 0.2 b	17.0 ± 1.0 b	69.3 ± 1.5 a	10.8 ± 0.2 a	34.3 ± 2.8 b
	September	63.1 ± 1.2 a	7.7 ± 0.2 a	19.0 ± 1.1 a	55.4 ± 0.9 b	10.0 ± 0.2 b	45.8 ± 2.0 a
Significance	Main effects
	G	***	***	***	***	***	***
	CA	***	*	ns	***	***	***
	HT	***	***	***	***	***	***
	G × CA	***	***	***	***	***	***
	G × HT	***	ns	***	***	***	***
	CA × HT	ns	***	***	***	ns	***
	G × CA × HT	***	***	***	***	***	***

Values followed by identical letters are not significantly different for p < 0.05, according to the LSD post-hoc test. The significance of variability factors according to the F-test: ns, not significant; *, significant at p ≤ 0.05; ***, significant at p ≤ 0.001 level. HI, harvest index; SLW, specific leaf weight.

).

Taking into account the genetic effect (G), plant height ranged from 47.31 to 58.19 cm, reaching the highest values in RG and NU genotypes and the lowest values in SL, PL, and SW30. On the other hand, SW30 was characterized by the greatest stem basal diameter (8.80 mm), while BR16 showed the lowest value (5.06 mm), although along with RG, BR16 had the greatest branching than the other genotypes (PL, SL, NU, CO, BR5, SW30, and BR1) (Table 4).

Regarding leaf dry yield per plant, NU and BR1 were the most productive, while BR16 appeared to be the least productive genotype. The other genotypes showed medium leaf dry yields. Harvest index ranged from 73.51% (SL) to 52.79% (BR16), reflecting the plant architecture, in terms of height and branching, and leaf yield. PL and SL exhibited the highest specific leaf weight (SLW), while RG and NU showed the lowest weights.

In terms of crop age, all biometric and productive measurements showed significant increases in the second year of cultivation, except for the stem basal diameter, which was significantly smaller in 2-year-old plants, and branching, which did not vary in the two years of cultivation (Table 4).

Additionally, harvest time played a key role in defining the agronomic responses of stevia genotypes. In general, maximum plant height, stem basal diameter, branching and leaf dry yield increased from the point of full vegetative growth to the beginning of flowering (Table 4). On the other hand, the highest values of specific leaf weight and harvest index were recorded during the full vegetative growth, with a significant decrease in the subsequent harvest.

3.2. Phytochemical Evaluation and Antioxidant Activities

Table 5. Results of three-way ANOVA for secondary metabolites (total phenols, flavonoids, and steviol glycosides) and antioxidant activities of stevia leaf extracts.

Factor	Factor Level	Total Phenols (mg GAE g−1 DW)	Total Flavonoids (mg CE g−1 DW)	FRAP (mmol TE g−1 DW)	DPPH (mmol TE g−1 DW)	Total SVglys (g 100 g−1 DW)
Main Effects
Genotype (G)	PL	58.2 ± 5.4 e	54.5 ± 3.9 d	0.34 ± 0.02 e	0.26 ± 0.03 d	25.2 ± 2.1 b
	BR16	41.6 ± 4.7 f	41.3 ± 4.6 e	0.29 ± 0.03 f	0.22 ± 0.02 e	16.2 ± 1.5 d
	RG	37.3 ± 3.5 g	33.3 ± 3.2 f	0.28 ± 0.02 f	0.20 ± 0.02 f	15.5 ± 1.9 de
	SL	71.3 ± 7.0 a	74.1 ± 5.2 bc	0.47 ± 0.02 bc	0.36 ± 0.04 a	24.6 ± 1.4 b
	NU	38.9 ± 3.7 fg	39.9 ± 2.9 e	0.27 ± 0.03 f	0.23 ± 0.02 e	12.2 ± 0.5 e
	CO	63.7 ± 6.6 cd	74.9 ± 7.4 b	0.45 ± 0.03 c	0.32 ± 0.04 c	19.7 ± 1.4 c
	BR5	66.1 ± 5.5 bc	81.2 ± 8.4 a	0.49 ± 0.03 ab	0.34 ± 0.04 b	23.9 ± 2.6 b
	SW30	61.5 ± 5.9 d	70.9 ± 6.9 c	0.42 ± 0.04 d	0.31 ± 0.04 c	29.3 ± 2.7 a
	BR1	67.9 ± 6.5 ab	84.1 ± 7.3 a	0.51 ± 0.02 a	0.31 ± 0.04 c	19.5 ± 1.6 c
Crop age (CA)	2018	32.2 ± 1.3 b	37.3 ± 1.8 b	0.32 ± 0.01 b	0.17 ± 0.003 b	24.3 ± 1.3 a
	2019	78.3 ± 2.1 a	83.8 ± 3.2 a	0.46 ± 0.02 a	0.39 ± 0.01 a	17.1 ± 0.6 b
Harvest time (HT)	July	58.7 ± 2.8 a	64.2 ± 3.6 a	0.47 ± 0.02 a	0.31 ± 0.02 a	18.5 ± 0.6 b
	September	53.2 ± 3.7 b	58.3 ± 4.0 b	0.32 ± 0.01 b	0.26 ± 0.01 b	22.5 ± 1.3 a
Significance	Main effects
	G	***	***	***	***	***
	CA	***	***	***	***	***
	HT	***	***	***	***	***
	G × CA	***	***	***	***	***
	G × HT	***	***	***	***	**
	CA × HT	***	***	***	***	***
	G × CA × HT	***	***	***	***	***

Values followed by identical letters are not significantly different for p < 0.05, according to the LSD post-hoc test. The significance of variability factors according to the F-test: ns, not significant; **, significant at p ≤ 0.01; ***, significant at p ≤ 0.001 level. FRAP, ferric reducing antioxidant power; DPPH, 1,1-diphenyl-2-picrylhydrazyl; Total SVglys, total steviol glycosides.

shows data on the secondary metabolites and antioxidant activities, depending on genotype, crop age and harvest time. Total phenols and flavonoids, total SVglys, as well as in vitro antioxidant activities (measured by FRAP and DPPH assays) were significantly affected by all variability factors (G, CA, and HT) and their interactions. Regarding the effect of genotype, SL showed the highest values of both total phenols and DPPH. Conversely, the highest total flavonoids and FRAP capacity were observed for BR5 and BR1. Finally, SW30 showed the highest SVglys total content.

Interestingly, compared to the other genotypes, BR16, RG, and NU exhibited the lowest values of all secondary metabolites (total phenols, total flavonoids and total SVglys), as well as the lowest in vitro antioxidant activities (Table 5). Crop age (CA) positively influenced the phytochemical characteristics of stevia leaf extracts, with an increase in their values in the plants in the second year after transplanting. Only the total SVglys showed a significantly decrease passing from the first to the second year of cultivation. Taking into account harvest time, total phenols, total flavonoids, and the antioxidant activities decreased from the vegetative phase (July) to the beginning of the reproductive one (September). On the other hand, an opposite trend, with a significant increase from the 1st to the 2nd harvest time, was observed for total SVglys, confirming that the development of the first flower buds corresponds to the time of major accumulation of SVglys in stevia plants.

Table 6. Results of three-way ANOVA for steviol glycoside profile.

Figure	Factor Level	Rub	Dulc A	Stev	Reb C	Reb A	Reb E	Reb M	Reb D
Main Effects
Genotype (G)	PL	0.59 ± 0.12 c	1.01 ± 0.04 b	12.05 ± 0.66 b	1.29 ± 0.05 b	7.47 ± 0.53 a	0.93 ± 0.14 bc	0.09 ± 0.001 b	1.36 ± 0.05 bc
	BR16	1.53 ± 0.12 b	1.71 ± 0.22 ab	3.86 ± 0.23 de	1.59 ± 0.19 b	5.15 ± 0.34 b	0.87 ± 0.06 c	0.08 ± 0.001 b	1.43 ± 0.03 ab
	RG	2.77 ± 0.38 a	1.93 ± 0.19 ab	1.19 ± 0.14 ef	5.31 ± 0.33 a	0.88 ± 0.08 c	1.19 ± 0.16 ab	nd	1.37 ± 0.01 bc
	SL	2.47 ± 0.57 a	2.22 ± 0.65 ab	8.83 ± 0.23 c	1.72 ± 0.10 b	7.55 ± 0.20 a	0.89 ± 0.01 c	nd	1.25 ± 0.02 bcd
	NU	2.19 ± 0.53 a	2.08 ± 0.22 ab	0.94 ± 0.09 f	5.61 ± 0.34 a	1.01 ± 0.19 c	1.35 ± 0.27 a	nd	nd
	CO	0.92 ± 0.02 bc	2.36 ± 0.25 a	4.51 ± 0.17 cd	5.70 ± 0.29 a	4.49 ± 0.31 b	0.94 ± 0.15 bc	nd	nd
	BR5	0.92 ± 0.04 bc	2.83 ± 0.52 a	5.71 ± 0.55 cd	7.79 ± 1.19 a	6.10 ± 0.81 ab	1.12 ± 0.25 ab	0.27 ± 0.01 ab	0.74 ± 0.01 d
	SW30	0.84 ± 0.04 c	1.76 ± 0.20 ab	16.02 ± 1.13 a	2.07 ± 0.30 b	4.39 ± 0.46 b	1.38 ± 0.26 a	0.54 ± 0.23 a	1.97 ± 0.19 a
	BR1	0.86 ± 0.14 c	2.31 ± 0.31 a	4.50 ± 0.29 cd	5.42 ± 0.38 a	4.68 ± 0.50 b	0.88 ± 0.13 c	nd	1.60 ± 0.01 a
Crop age (CA)	2018	0.40 ± 0.10 b	2.86 ± 0.28 a	7.17 ± 1.49 a	4.87 ± 0.87 a	5.30 ± 0.75 a	1.40 ± 0.13 a	0.81 ± 0.05 a	1.63 ± 0.18 a
	2019	2.03 ± 0.21 a	1.23 ± 0.06 b	5.54 ± 0.96 b	3.42 ± 0.40 b	3.84 ± 0.53 b	0.61 ± 0.04 b	0.11 ± 0.02 a	1.32 ± 0.10 a
Harvest time (HT)	July	1.50 ± 0.32 a	1.75 ± 0.17 b	5.97 ± 1.13 a	3.59 ± 0.48 b	3.88 ± 0.53 b	0.88 ± 0.11 b	0.11 ± 0.001 a	1.35 ± 0.08 a
	September	1.20 ± 0.24 a	2.19 ± 0.35 a	6.74 ± 1.43 a	4.48 ± 0.85 a	5.05 ± 0.76 a	1.26 ± 0.18 a	0.39 ± 0.02 a	1.69 ± 0.25 a
Significance	Main effects
	G	**	ns	***	**	**	**	***	***
	CA	***	***	*	**	*	***	**	ns
	HT	ns	*	ns	*	**	*	ns	ns
	G × CA	**	ns	**	**	*	ns	***	ns
	G × HT	ns	ns	ns	**	**	ns	ns	ns
	CA × HT	ns	*	**	**	**	**	**	**
	G × CA × HT	***	ns	***	***	**	**	***	ns

Values followed by identical letters are not significantly different for p < 0.05, according to the LSD post-hoc test. The significance of variability factors according to the F-test: ns, not significant; *, significant at p ≤ 0.05; **, significant at p ≤ 0.01; ***, significant at p ≤ 0.001 level. Rub, Rubusoside; Dulc A, Dulcoside A; Stev, Stevioside; and Reb A, C, D, E, and M, Rebaudiosides A, C, D, E, and M.

reports the SVglys profile and related statistical significance. A significant effect of genotype was observed for almost all steviol glycosides, with the exception of Dulcoside A. Interestingly, the identified compounds were not present in all genotypes. This difference was evident for Reb M which was not detected in RG, SL, NU, and CO, and for Reb D which was not detected in NU and CO. The reduction in SVglys content observed in plants at the second year of age was also seen in the most represented compounds such as Stev, Reb C, Reb A, and Reb C. On the other hand, Rubusoside increased with crop age, while no effect was observed for Reb M and Reb D. With regard to harvest time, as already noted for the total content, each steviol glycoside significantly increased from the vegetative phase to the beginning of flowering, except for Rubusoside which remained stable in both harvest times.

HCA and PCA were carried out for total phenols, flavonoids, and SVglys content. The two-way dendrogram of the HCA is reported in Figure 2. The first macro-cluster (red) was grouped by itself, while the second comprised three sub-clusters (green, blue, and yellow). Among the metabolites analyzed, the total contents of phenols and flavonoids were clustered together, while the total SVglys concentration was grouped by itself, thus highlighting a higher degree of dissimilarity based on this parameter. Based on their total SVglys concentration, samples in the yellow cluster, all from 2018, shared the highest SVglys content; samples of the red cluster, all from 2019, were, instead, characterized by intermediate concentrations of SVglys; the lowest concentrations of these metabolites were, instead, common to samples of the blue and green clusters, especially those from 2019. As evidenced by their grouping in a common macro-cluster, the total phenols and total flavonoids showed a common quantitative distribution among the samples. In fact, higher concentrations of both these chemical classes were evidenced for samples of the red cluster, all from 2019; also, intermediate concentrations were evidenced for samples from 2019 in the green cluster; finally, samples belonging to the blue and yellow clusters exhibited the lowest contents of these compounds. With the exception of the green cluster, the samples were homogeneously distributed in the groups based on the harvest time.

Some genotypes showed very similar total abundances of the analyzed compounds between the two harvest times (July and September) within the same year: i.e., CO 2019, NU 2018 and 2019, RG 2018 and 2019, SW30 2018 and 2019, BR5 and BR1 2019 (Figure 2). The score plot of the PCA (Figure 3) confirmed the HCA grouping: with the exception of two samples (PL_2_2019 and BR5_2_2018), only the red samples were grouped in the upper quadrants (PC2 > 0), due to their higher phenols and flavonoids contents.

Among the other groups, the higher content of SVglys detected in the yellow samples meant that they were plotted in the right quadrants (PC1 > 0), whereas all the other samples were in the left quadrants (PC1 < 0). The same multivariate analyses were conducted on the SVglys profile of all the genotypes studied (Figure 4 and Figure 5). The two-way dendrogram (Figure 4) of the HCA on the compounds showed two main groups driving the dissimilarities: Rubusoside, Stevioside, Rebaudiosides A, M, and D for macro-cluster 1; Rebaudioside E, Rebaudioside C, and Dulcoside A for macro-cluster 2. The samples in the dendrogram (Figure 4) were distributed in two macro-clusters: the first comprised the red and green sub-groups; the second comprised the blue and yellow sub-groups, of which the latter only contained one sample (SW30_2_2018). The second macro-cluster was mainly composed of 2018 samples, while the first macro-cluster comprised all the 2019 samples, as well as a few of the 2018 ones (Figure 4).

Compared to the HCA of the total compounds, the dendrogram of the SVglys profile suggested that the differences between the samples were less due to the year and more to the steviol glycosides pool produced by each genotype (i.e., all the PL samples are in the red cluster; all the BR16 samples are in the green cluster). In addition, the SVglys pool of samples RG and NU appeared very similar, thus suggesting the proximity of these two genotypes; the same seemed true for BR5, BR1, and BR16 (Figure 4). Regarding SVglys profile, the influence of the year on the sample distribution was also confirmed by the PCA. The score plot is reported in Figure 5. Samples belonging to the first HCA macro-cluster, thus almost all 2019 samples were mostly distributed in the bottom quadrants (PC2 < 0); in particular, green samples were almost all grouped in the bottom left quadrant (PC1 < 0), while red ones were all positioned in the bottom right quadrant (PC1 > 0). The blue and yellow HCA clusters, so all 2018 samples were instead almost all plotted in the upper quadrants (PC2 > 0); blue samples were chiefly plotted in the left quadrant (PC1 < 0), while the only yellow sample was plotted in the right one (PC1 > 0). The similarity evidenced for genotypes RG-NU and BR5-BR1-BR16 was confirmed by the PCA plot (Figure 5); the former (RG-NU), in particular, were grouped closely in the upper left quadrant in the 2018 genotypes, while the 2019 ones were plotted together in the bottom left one. Sample BR5_2_2018 is plotted quite separately from all other samples, in the upmost area of the upper right quadrant, due to its high Dulcoside A concentration.

4. Discussion

We assessed the agronomic and phytochemical performances of nine S. rebaudiana (Bertoni) genotypes for two consecutive growing seasons. The aim was to select the best productive genetic resources to develop new and improved cultivars for the Mediterranean area. To improve the competitiveness of stevia production as a semi-perennial crop in this environment, it is important to produce higher-performing crops in terms of yield and quality. Therefore, the identification of new varieties/cultivars of S. rebaudiana with long stand duration, consistently high and stable leaf and SVglys yields, as well as a high level of other beneficial substances, is the top priority for the successful introduction of stevia into Mediterranean cropping systems.

Our results showed that some of the genotypes tested were characterized by a longer vegetative period before flowering, as demonstrated by the higher thermal requirements (Figure 1). This behavior promotes both higher leaf yield and steviol glycoside concentration, which in turn determine higher SVglys yield. This is important given the time required for the plant to synthesize and accumulate the steviol glycosides in the leaves. In fact, the maximum accumulation of these compounds, which depends on several environmental, agronomic, and physiological factors, is reached when the plant moves from the budding phase to an initial flowering stage with less than 10% flowers [31]. Therefore, the greater the vegetative period, the greater the accumulation of steviol glycosides in the leaves. When the plant starts to flower, nutrients accumulate in the reproductive organs and, as a result, vegetative growth declines.

Given that the leaves are the commercially important part of stevia, a delay in flowering can enhance vegetative growth and economic yield. In this regard, Ceunen and Geuns [31] observed large amounts of SVglys within the upper leaves during the budding phase and, as more and more buds become flowers, leaves contained lower amounts of SVglys. This is important for the choice of best harvesting time, since the maximum SVglys concentration is reached in the physiological stage of flowering-bud formation [32,33]. Flowering and, in general, plant growth are also affected by the photoperiod and temperature. Stevia is, in fact, a short-day plant with a critical photoperiod of between 12–13 h of day length [34,35]. Regarding air temperature, in our environment vegetative growth was lower when the maximum day temperature was below 10 °C or over 35 °C.

Overall, our results showed that genotype, crop age and harvest time represented key pre-harvest factors for defining the morphological, phenological, and quanti-qualitative traits in stevia.

Regarding differences in plant morphology and canopy architecture, we confirmed previous findings [21,36,37], which reported a very high variability, due to genetic characteristics, has been reported for plant height, basal stem diameter, branching, specific leaf weight, and leaf yield. In order to investigate this great variability among stevia genotypes and populations, some authors [21,22] have observed that increased yields were correlated to a high leaf area index (LAI), which is responsible for a greater light interception and, in turn, for a higher leaf photosynthesis.

We assessed the specific leaf weight (SLW), a leaf thickness index that is positively correlated with leaf photosynthesis: leaves with high SLW values are thicker and generally have a higher chlorophyll density per area unit (µg cm⁻²) and, therefore, greater photosynthetic capacities than thinner leaves [38]. SLW significantly changed depending on genotype, crop age and harvest time. In particular, in 2-year-old plants, an increase in SLW was accompanied by increased leaf dry yields, suggesting that, thanks to a greater SLW leaf photosynthesis may have been enhanced with a consequent increase in crop yield [39].

Plants with a high SLW have leaves with a low surface/volume ratio, which is more efficient in terms of water use [40]. This ratio plays an important role in leaf functioning and is related to the strategies for acquiring and using the energy resources [41], as well as a tool to evaluate plant productivity [42] with a view to sustainable agriculture. For a better comprehension of the relationships between morphotype and light interception and photosynthetic activity, future studies are needed in order to improve leaf yield in S. rebaudiana.

Regarding stevia productivity, we found that the least productive genotypes are those with more branching, as reported by Tateo et al. [43]. In addition, our findings suggested a positive relation between plant height and leaf dry yield in agreement with previous studies [36,37]. Finally, the aerial biomass variability among genotypes exhibited high heritability, and the differences can also be partly attributed to crop age, pedo-climatic conditions of the cultivation site and plant development stage, as already observed in previous studies [19,44,45].

In our environment, stevia productivity significantly increased from the first to the second year of cultivation, with a very high winter survival rate. This suggests that stevia could grow as a semi-perennial crop in central Italy. This could bring positive agronomic advantages, such as, reductions in soil erosion and nutrient leaching, C soil sequestration, which in turn have important implications in multifunctional and sustainable agriculture.

Our findings confirmed previous results [5,6] in which in the temperate area of central Italy, the cultivation of stevia is long-term (5–7 years), with a vegetative period from April–May to September-October. In the tested climatic conditions, during winter, above-ground parts of the plant became dry, and there was regrowth through new shoots in the buried rhizome the next spring, thus producing a new crop without replanting. On the other hand, at higher latitudes as in central Europe (such as Belgium and Germany) and Canada [46,47,48], with cooler winters, stevia behaves like an annual crop with annual replanting.

We thus believe that stevia can be cultivated as a promising new crop in the Mediterranean climate of central Italy, though the SVgly content and composition still need to be optimized.

We found that the content and composition in bioactive compounds varied greatly depending on the genotype, crop age and harvest time. Two-year-old plants generally exhibited the highest content of total phenols and flavonoids, as well as the antioxidant activities; conversely, total SVglys content peaked in the establishment year. The sampling carried out at the beginning of the reproductive phase only improved total SVglys, again confirming that, in stevia, plant development is a primary factor in influencing the biosynthetic pathway.

Harvesting when plants were in full vegetative development (July harvest) maximized the content of polyphenol compounds and improved the antioxidant activities. Similar results have been obtained in a previous work [49] highlighting how the identification of the optimal harvest time was able to maximize the bioactive compounds of interest and, consequently, the health-promoting properties of stevia leaves. However, since the biosynthesis of secondary metabolites is a plant defense mechanism against biotic and abiotic stresses, through many physiological, biochemical, and molecular changes in plant metabolism, it is not always easy to identify the optimal time in which the different metabolites are maximally expressed.

The accumulation of polyphenols during the beginning of flowering stage could be related to the ecological roles of these compounds, such as intensifying antifungal defenses and attracting pollinators [32]. In addition, since different interactions among environmental and agronomic factors may occur, it is difficult to select individual stimuli that can influence a single metabolic pathway.

Of the secondary metabolites that have been synthesized and accumulated in stevia leaves, the presence of phenolic compounds is the subject of increasing interest because of their significant practical use for nutritional and medicinal applications. In fact, there are important implications for the growing market of natural stevia products [50], where they are employed as natural preservatives, thanks to their capacity to delay the oxidative degradation of lipids and to improve the shelf life of foods and beverages. They are also involved in the prevention of oxidative stress in humans, thanks to the hydroxyl groups in their molecules which have antioxidant anti-inflammatory properties [51,52].

Our findings provide new knowledge about the dissimilarity among genotypes, which was above all generated by the total SVglys content with respect to total phenols and flavonoids. Previous studies have investigated the effect of farming practices and genotype × environment interaction on stevia crop productivity, Reb A and Stev accumulation and on their reciprocal ratio (Reb A/Stev) [19,21,22,53,54,55]. Our results showed that the differences in the steviol glycoside profile are less due to crop age and more to genetics, as suggested by the dendrogram of the SVglys profile. It is well-known that, while the SVglys content in the leaves can vary depending on environmental and agronomic factors [56,57], the SVglys qualitative profile remains quite stable, indicating the high genotypic determinism of this trait [19].

We thus believe that our findings are a good starting point to screen the best genotypes in terms of SVglys profile, taking into account that breeding programs for stevia are increasingly aimed at varieties with an optimal Reb A/Stev ratio, which is considered a good qualitative measure of sweetness, or with a high content of Reb M and Reb D, characterized by a very sweet taste with no liquorice-like taste. We found that genotypes, such as SL, BR5, PL, and SW30 accumulated significantly high amounts of Reb A, Reb M, and Reb D. Among these genotypes, SL and BR5 were characterized by relatively low content of Stevioside. Conversely, PL and SW30 exhibited high contents of Stevioside, which negatively affect the Reb A/Stev ratio, lowering it below 1.

In our study, we try to identify the stevia ideotype, deriving from the combination of morphological, productive, and phytochemical traits. Consequently, considering not only the SVglys profile, but also leaf yield, growth crop cycle and the composition of their whole phytocomplex (polyphenols, total SVglys and their related antioxidant activities), SL, BR5, PL, and SW30 seemed to be the best performing genotypes, in the given environment, and thus suitable as starting point for future breeding programs.

5. Conclusions

Our findings highlighted that stevia, in the tested environment, represents a promising semi-perennial crop, which can contribute to the diversification of traditional cropping systems, thus increasing their sustainability and, at the same time, generating functional products with high added value.

We revealed a strong dependence of growth, crop yield and quality on genotype, harvest time, and crop age. The high variability among genotypes highlighted the importance of morpho-productive and qualitative traits for the identification of the best stevia ideotype, for our environment and for targeted end uses. The best combination of morphological, productive, and phytochemical traits was observed for SL, BR5, PL, and SW30 genotypes.

In conclusion, identifying genotypes characterized by specific SVglys profile (high Reb A, Reb D, Reb M), together with high content of all bioactive compounds and satisfactory yields, is crucial to establish future breeding programs for this crop.