Effect of Planting Density in Two Thistle Species Used for Vegetable Rennet Production in Marginal Mediterranean Areas

: In Mediterranean cropping systems, it is important to utilize marginal lands for the cultivation of non-food crops. Spontaneous thistle genera such as Cynara cardunculus L. and Onopordum tauricum Willd. are native to southern Europe. Previous research has focused on their spontaneous growth in the environment or ability to supply biofuel. The aim of this study was to identify the effect of two different planting densities on the ﬂower differentiation, vegetable rennet production and weed control of O. tauricum Willd. and C. cardunculus L. in rainfed unfertilized Mediterranean cropping systems. The results showed that plant density inﬂuenced most of the biomass parameters for O. tauricum Willd. but no signiﬁcant differences were observed for C. cardunculus L. Higher densities of both species were beneﬁcial for weed control. We recommend the use of C. cardunculus L. and O. tauricum Willd. as suitable species for vegetable rennet production in Mediterranean organic cropping systems.


Introduction
In Mediterranean areas, to increase food productivity, it is important to utilize some marginal land for the cultivation of non-food crops (including livestock breeding, forestry and cultivation), thus reducing competition for land with food crops [1]. It is also important to use surplus and marginal lands to promote regional economic structures, provide alternative sources of employment in rural areas, reduce CO 2 emissions, and decrease vulnerability to short-term changes in weather [2].
In the Mediterranean region, it is recognized that the water availability is one of the most limiting factors for plant growth and productivity under rainfed conditions [1]. Water availability strongly affects the choice of crops in Mediterranean environments [3] where only few options are suitable. In addition, the choice and the introduction of crops should consider the existing cropping systems and plant needs at a local level [4]. In marginal areas, it is necessary to identify the most promising species in terms of adaptation capacity, required inputs, and yields [4]. Most of these species belong to the Asteraceae (Compositae), which is a large and widespread family of vascular plants represented by 1600 genera and about 23,000 species worldwide [5]. In Mediterranean cultivation systems, native spontaneous thistle genera such as Cynara cardunculus L. and Onopordum tauricum Willd. require little additional irrigation because they are well adapted to the xerothermic conditions of southern Europe, which is characterized by dry and hot summers [6,7]. Their adaptation to southern European environmental conditions and their high yield suggests that these species could be useful for industrial [3,[8][9][10], technological food [11] and pharmacological [12][13][14] use. Cynara cardunculus L., also known as Spanish thistle, artichoke or cardoon, is a perennial crop of Mediterranean origin [15,16] which is well adapted to environmental conditions in the region owing to the positive balance of the phases of its growth cycle with the climatic trends [17,18]. In addition, its long vertical tap-root allows the plant to access a large volume of soil so that it requires low cultivation inputs [19,20]. This valuable crop shows high yields and drought tolerance, as well as benefitting soil properties, erodibility, and biological and landscape diversity [21,22]. It is a robust plant with a characteristic rosette of large spiny leaves and branched flowering stems [6,23]. Cynara cardunculus L is a multipurpose crop that can be utilized as a raw material in paper pulp industry, as forage in wintertime, but most importantly as solid and/or liquid biofuel in the bio-energy sector [3]. Extracts of plant proteases have been used for centuries as milk coagulants in cheesemaking, especially for the curdling of raw ovine and caprine milks. The flowers of wild cardoon are traditionally used as a milk-clotting agent to produce cheeses in countries such as Portugal, Spain, Morocco, and Italy [24]. The use of vegetable rennet is geographically circumscribed and limited to a few products, but the scarcity and high price of conventional animal rennet has led to a growing interest in the vegetable sources of milk-clotting enzymes [11].
Onopordum L. is a genus of about 60 species of thistles belonging to the family Asteraceae, which are native to Europe (mainly the Mediterranean region), northern Africa, the Canary Islands, the Caucasus, and southwest and central Asia [5,25]. Onopordum tauricum Willd. is native to southern Europe and southwest Asia. It is a branched, robust biennial or annual herb with stems up to 2 m with dense erect glandular hairs, capitula solitary, terminal, on short peduncles. It is readily distinguished by its green, sticky, glandular leaves, stems, and bracts [26]. These parameters distinguish O. tauricum from two similar species, O. acanthium and O. illyricum. In several locations, for example, in Colorado, putative hybrids between O. acanthium and O. tauricum Willd. have been observed where the two species cohabit [26]. Most studies have focused on O. acanthium L. and O. illyricum L. The former was used traditionally as a bactericide, cardiotonic, haemostatic, diuretic, and antitumor agent, as well as to treat nervousness, and inflammation of the bladder and the respiratory and urinary systems [27]. Pharmacological studies show that O. acanthium L. possesses antibacterial, antioxidant, anticancer, anti-inflammatory, analgesic, antipyretic, hypotensive, antiepileptic, wound healing, xanthine oxidase, and ACE inhibitory effects [27].
As shown by the work on C. cardunculus L. in central and southern Italy [4,28,29], research on thistle species with regard to planting density has mainly examined the chemical composition of biomass components (stalks, leaves and heads) in different locations with fertilizer inputs over several years. The present study is part of the European project (PRIMA 2018-Horizon 2020) "Valorization of thistle-curdled CHEESEs in MEDiterranean marginal areas", acronym VEGGIE-MED-CHEESES (https://veggiemedcheeses.com/).
Its aim was to evaluate the effect of different plant densities on the flower differentiation, vegetable rennet production and weed control of O. tauricum Willd. compared with C. cardunculus L. in rainfed unfertilized Mediterranean cropping systems to assess the potential of these species, given the lack of agronomic information in the scientific literature.

Experimental Site
The experimental site is located at the Pasquale Rosati experimental farm of the Marche Polytechnic University in Agugliano, Italy (43 • 32 N, 13 • 22 E, at an altitude of 100 m above sea level and a slope gradient of 10%) ( Figure 1). The site has a silty-clay soil classified as Calcaric Gleyic Cambisol [30].
The weather data were provided by ASSAM (Agenzia Servizi al Settore Agroalimentare delle Marche, Osimo Stazione, Ancona, Italy). The climate of the site is Mediterranean (Table 1).  Table 1. Weather monthly data comparison for temperature (Tavg, Tmax and Tmin ( • C)) and precipitation (mm) in the period 1998-2019 and in 2020 at the experimental site.

Month
Tavg The annual average temperature was 17 • C, and the monthly mean temperature increased from March (10 • C) to August (26 • C) and then decreased. The mean minimum temperature of the coldest month was 3.8 • C (January) and the mean maximum temperature of the warmest month was 29.4 • C (August).
The annual rainfall was 359 mm vs. 564 mm of the historical period, distributed mainly from March to May (43%). The winter (January and February) recorded 21.4 mm of precipitation (6% of the precipitation observed during the period), in strong contrast to the historical average 1998-2016 (100 mm).
Soil characterization was performed immediately before transplanting (Table 2) on a total of three soil samples that were taken at a depth of 0-30 cm. Table 2. Soil properties in the 0-30 cm layer in the experimental plots in 2020.

Experimental Design and Crop Management
The soil in the experimental site was ploughed along the maximum slope using a mouldboard (with two ploughs) to a depth of 40 cm in autumn. No fertilization or chemical-synthetical plant protection measures were applied. The dates of all agronomic practices are reported in Table 3. Table 3. Agronomic management practices using during the experimental period (dd/mm/yyyy).
The O. tauricum Willd. seeds were collected in the Sibillini Mountains of the Marche region while the C. cardunculus L. seeds came from the Spanish region of Murcia.
O. tauricum Willd. and C. cardunculus L. seeds were sown between 1 and 10 February 2020 in a greenhouse with an average temperature of 20 • C. The seeds were planted in biodegradable pots (type: Ø 16) in a substratum consisting of a mix of topsoil and peat (50:50 ratio) and periodically moistened. To characterize the phenological development of the crop we adopted the Bundesanstalt, Bundessortenamt, CHemischeIndustrie (BBCH) scale [31] during the full crop biological cycle. When three or four ellipsoidal leaves were already visible (phenological phase BBCH 13 and BBCH 14) the field transplantation was performed.

Phenological Phases
According the BBCH scale, we conducted phenological monitoring once a week from the day of transplanting until the flower heads were in full bloom as growth stages from BBCH 59 to BBCH 69. This monitoring was essential for planning the emergency irrigation and organizing the harvesting of the flower heads.

Plant Height and Number of Branches
To evaluate the influence of different plant density on the crop biometric characteristics, we measured the plant height in the full flowering phase (BBCH 69), when the main stem and the branches reached their maximum development. In addition, in the same phenological phase, we counted the primary ramifications emerging from the main stem for each plant.

Number and Weight of Flower Heads
To determine the number of flower heads per square meter we counted the flower heads on each branch (in the apical portion and along the branches themselves). Once the flower heads were harvested, the fresh weight was measured to determine the flower weight (g m −2 ). Furthermore, each flower head was dissected, to determine the weight of the flower tubules (g and g m −2 ) which are the potential source of vegetable coagulant. Other authors [32] have reported the average weight of receptacles, bracts and pappi, which totals about 80-85% of the average weight of the flower heads; the tubules are about 15-20% of the average weight.

Number of Harvesting Days
To evaluate whether the different planting density had an effect on the differentiation of flower heads, we measured the number of harvesting days (the temporal phase of the plant's growth cycle in which the flower heads were completely open and could be cut at the stalk level and stored for laboratory use (phenological phases between BBCH 61 and BBCH 69).

Weed Biomass Determination
To evaluate the effect of plant density on weed development, during the maximum thistle leaf development stages (BBCH 45-BBCH 49), the weed biomass was cut at 3 cm from the soil surface in the inter-row in three 0.25 m 2 test areas. It was then weighed fresh and subsequently oven dried at 105 • C for 48 h to determine its dry weight expressed in kg ha −1 .

Statistical Analysis
All statistical analyses were performed with R statistical software [33]. Good analysis is based on a model that is biologically relevant and uses realistic assumptions [34]. To select the model that could fit the crop parameters better, we built two different models by using the "stats" [35] and "nlme" [36] R packages: (1) A one factor linear model was built by using the "lm" function in which the species was considered the main factor. (2) A full factorial linear model was also created by using the "lm" function in which the species and density were considered as factors.
In order to evaluate the best model, we made an 'posteriori' selection, based on statistics which put a penalty on 'complexity', such as the Akaike Information Criterion (AIC) [37] the Bayesian Information Criterion (BIC) [38] and likelihood ratio tests (LRTs) [39]. Based on the previous model analyses, the full factorial model better fitted the crop parameters.
Before performing the analysis of variance (ANOVA), we checked that the full factorial model met the three assumptions of the ANOVA. The normality distribution of the model residuals was checked both graphically (QQ-plot) and by performing the Shapiro-Wilk normality test. The homoscedasticity was checked using the Levene test. The final ANOVA assumption was satisfied by the experimental design and random sampling.
When all three ANOVA assumptions were met, we applied the ANOVA to the model. When the ANOVA showed a significant difference between the factors (p-value < 0.05), we performed an estimated marginal means post-hoc analysis by using the "emmeans" function with the Bonferroni adjustment of the "emmeans" R package [40].

Effect of Species and Plant Density
The results of the ANOVA (Table 4) showed significant differences associated with thistle species (S), plant density (D) and their interaction (S × D) for almost all measured variables. Table 4. Analysis of variance (ANOVA) applied to the mixed model. The thistle species showed a significative difference at p < 0.001 for all analyzed variables except phenological growth stage (PGS). For this last parameter, C. cardunculus L. and O. tauricum Willd. showed a similar phenological development. The response to plant density was significantly different (p < 0.001) only for the average number of flower heads per plant (FHP). Less significant differences from the statistical point of view (p < 0.01) were found for the development of PGS and plant height (HPT). These differences were repeated in the same proportions even in the interaction of the two factors, S × D. Regarding biometric variables ( Table 5), neither of the compared densities (D1 and D2) for C. cardunculus L. showed a significant effect. Table 5. Estimated biometric marginal mean differences calculated with the Bonferroni adjustment between the compared treatments.
Regarding the weight of flower heads and tubules (Table 6), plant density did not affect their weight on C. cardunculus L., and there were no significant differences per unit area. In contrast, O. tauricum Willd. showed a significant difference in the weight of flower heads per square metre (higher weight in D1 compared with D2), although the other variables did not vary significantly with density. Table 6. Estimated biometric marginal mean differences calculated with the Bonferroni adjustment between the density treatments relative to weights of flower heads and tubules.

TS PD WHF (g) WHFS (g m −2 ) WTF (g) WTFS (g m −2 )
Cc  We measured the quantity of weeds in the first ten days of June to assess biomass, i.e., a time when the plants started to develop flower heads along the stem and the primary and higher order branches. Both in C. cardunculus L., and in O. tauricum Willd. plant density had a significant effect on the weed biomass, with the biomass being lower in the D2 treatment (Table 7).

Linear Regression Analysis
We analyzed the correlation between the variables in relation to the single plant and to the combination of species and plant density.
As concerns biometric variables, there was a high correlation between plant height and number of branches (R 2 = 0.81) as well as between the number of heads produced by each plant and plant height (R 2 = 0.70) for C. cardunculus L. at D1 (Table 8). At D2 (Table 8), these correlations remained high. Compared to D1, the correlation between the number of flower heads and plant height increased by almost 20% (R 2 = 0.86) and, between plant height and number of primary branches, it decreased by almost 20% (R 2 = 0.66). At D2, the highest correlation was between the number of flower heads and the number of primary branches (R 2 = 0.94).
Onopordum tauricum Willd. at D1 (Table 9) had a high correlation between the number of flower heads and the number of primary branches (R 2 = 0.93). The same behavior was observed at D2 (Table 9) but the coefficient of determination values was lower, by about 20% (R 2 = 0.73). For the parameters referring to the unit per area (Table 10), the number of flower heads was highly correlated with the weight of the flower tubules and the weight of flower heads at the two plant densities for C. cardunculus L. The four correlations mentioned had values of R 2 > 0.98 on average.    The correlations were much lower at the two plant densities for O. tauricum Willd. The value of R 2 between the number of flower heads and their weight was higher in D2 (0.54) than in D1 (0.43), by about 20%. The situation was reversed in the case of the correlation between the number of flower heads and the weight of flower tubules, which was higher in D1 (R 2 = 0.54).

Discussion
In the described experimental context, it was observed that the morphometric characteristics were not affected by either of the two plant densities in Cynara cardunculus L. maintained similar values in terms of the number of flower heads produced, number of branches, average plant height, phenology, and harvest days in both densities. In Egypt, a similar height was reported for C. cardunculus L. at a plant density of 1 m × 1 m to that we observed in this study [41]. In the Basilicata region in southern Italy [9], a higher average plant height and number of flower heads per plant in the first year of transplanting at a density of 1 plant per m −2 with pre-fertilization and more rain was observed, compared to our experimental site. Also, in Sardinia, Italy, in the first growing season with input fertilization, the plant heights and number of flower heads on milk thistles of local genotypes were higher than in our experimental test [4]. In Sicily, southern Italy, however, Raccuia and Melilli [10] reported lower plant heights for most local genotypes than those recorded on our site with a plant density of 1 m × 0.5 m during the first year of transplantation, although the production of flower heads per plant was similar.
Several cultivated and wild genotypes were examined in Central Italy that showed higher plant heights than in our case study, but these were situated in cultivation systems with higher fertilization inputs and plant densities; on average, the number of flower heads was the same [32]. Interestingly, we found that the parameters that affected the plant-based rennet extract such as the quantity of flower heads and tubules produced per unit of land were maintained without significant differences in the two plant densities we tested. The main advantage for this species was that in the D2 density, which had plants that were closer together, there was a smaller quantity of weeds which reduced the competition for vital nutrients and water.
Onopordum tauricum Willd. showed different results from morphological parameters between the two plant densities. In D2, the parameter values related to the production of flower heads, branches, and height were lower because there was stronger competition with the neighboring plants and less tendency to develop reciprocally. However, even for this species, the thistles out-competed the weeds, which improved weed control by reducing the number of interventions. Another advantage was that with smaller plants, although they were significantly behind in terms of phenology compared to plants at D1 spacing, the flowering period was concentrated into a few weeks compared with larger plants that had longer flowering duration.
Given that there were no significant differences in flowering parameters within the same tested herbaceous species at different densities, we examined which species between C. cardunculus L. and O. tauricum Willd. was better in terms of the production of flower tubules per unit area as these are the most interesting part of the flowers for the production of vegetable rennet. C. cardunculus L. produced about 800 g m −2 (average of the values of D1 and D2) while O. tauricum Willd. produced about 3500 g m −2 (average of D1 and D2).
To make rennet, some authors have suggested that 2 g of dried flowers from C. cardunculus L. diluted in 50 mL of distilled water [42,43] can be used as a crude enzymatic extract. A concentration of 0.5% and 1% aqueous extracts was tested in cows' and goats' milk, respectively. In general, an increase in the cheese yield was observed in goats' milk rather than in cows' milk (202.6 g L −1 and 150.5 g L −1 ). Tubular flowers of O. tauricum Willd. need to be macerated in demineralized water (1:10 w/v) for 24 h at 4 • C. The aqueous crude extract is freeze-dried and reconstituted in demineralized water 1:10 w/v [11]. Further research, not published yet, shows that 20 mL of reconstituted extract is necessary to clot 5 L of ewes' milk with the addition of CaCl 2 solution (10 mM), obtaining a yield of around 50%.
According to these data, the O. tauricum Willd. could be used as a milk coagulant in cheesemaking, to design a valuable food product that utilizes marginal lands and stand ardizes the technological properties of the plant. However, further studies are necessary to better undestand the no specific proteolitic activity which determine the sensory properties of the obtained cheese. In the future it could be interesting to test both species in Mediterranean organic cropping systems also for bioenergy intents [44,45] to promote their diffusion in marginal areas otherwise underused.

Conclusions
This study confirmed the potential of Cynara cardunculus L. and Onopordum tauricum Willd. as herbaceous crops for vegetable rennet production in a Mediterranean environment in a niche cultivation system with low inputs, in terms of biomass parameters.
The compared plant densities did not affect the crop parameters of C. cardunculus L., which maintained similar values in terms of the number of flower heads produced, number of branches, average plant height, phenology, and harvest days in both tested densities. O. tauricum Willd., in contrast, showed different results between the two plant densities, with a lower production of flower heads in D2.
Other studies [6] suggest these species are competitive crops in their native Mediterranean environments, mostly because of their low water requirement. In this study, emergency irrigation was necessary in the weeks following transplantation because of a particularly hot and dry spring compared to the historical meteorological series.
Comparing the food technology aspects of vegetable rennet production in the literature suggests that, because of the similar levels of production of dry matter and flower tubules, the different technological properties in the two species did not allow us to establish which species was more advantageous. Therefore, we recommend the potential adoption of both C. cardunculus L. and O. tauricum Willd. as suitable species for vegetable rennet production in Mediterranean organic cropping systems.
The research activities that are still ongoing will focused on the characterization of thistle aqueous crude extracts and the compatibility of these extracts in the qualitative production of local cheeses.