Microbial Dynamics of a Specialty Italian Raw Ewe’s Milk Cheese Curdled with Extracts from Spontaneous and Cultivated Onopordum tauricum Willd

Milk coagulants prepared by maceration of flowers harvested from both spontaneous and cultivated Onopordum tauricum Willd. and a commercially available coagulant obtained from Cynara cardunculus L. (control) were assayed for small-scale manufacturing of Caciofiore, an Italian specialty raw ewe’s milk cheese produced in a family run dairy farm located in the Marche region (Central Italy). The microbiota of the three thistle-based milk coagulants and their effect on the microbial dynamics of raw milk cheeses during fermentation and maturation (from day 0 up until day 60) were investigated through a combined approach based on viable counting and Illumina DNA sequencing. In both the control and experimental cheeses, despite the slight differences emerged depending on the coagulant used, Lactococcus lactis and Debaryomyces hansenii were the prevalent species among bacteria and fungi, respectively. Moreover, raw ewe’s milk was the main factor affecting the evolution of both the bacterial and fungal microbiota in all cheeses. The overall similarities between control and experimental cheeses herein analyzed supports the exploitation of Onopordum tauricum Willd. as an alternative milk coagulating agent for production of Caciofiore and, more in general, raw ewe’s milk cheeses.


Introduction
Cheese is a complex food, whose characteristics largely depend on basic (milk, salt, starter culture, and rennet) and optional unconventional ingredients (e.g., grated lemon peel, aromatic herbs, etc.), as well as on the parameters applied during cheesemaking, which, in turn, affect the biochemical and microbiological transformations occurring during fermentation and ripening [1,2].
It is estimated that more than a thousand varieties of cheese are produced in the world by combining the aforementioned factors [3]; for approximately 75% of produced cheeses, coagulation is achieved through the addition of rennet to the milk [4].
Historically, the term "rennet" refers to the crude extract from bovine calf stomach containing a complex set of protease enzymes that curdles the casein [5]. However, starting from the 1950s, due to a shortage in the supply of animal rennet, the search for "rennet substitutes" began. Four main types of milk coagulants are currently available, being animal rennet, microbial rennet, fermentation-produced chymosin, and vegetable rennet [6].
Concerning the latter, the ability of some plants to coagulate milk has been known since ancient times; in his treatise De Re Rustica (ca. 50 AD), Lucius Columella was the first to suggest the use of wild thistle flowers, safflower seeds, and fig latex for milk (CE) from spontaneous (st) and cultivated (ct) flower heads were prepared following the procedure outlined by Mozzon et al. [24] and illustrated in Figure 1. Briefly, the purple tubular flowers were separated from the receptacle and macerated in demineralized water (1:10 w/v) for 24 h at 4 • C; the resulting aqueous extracts were cleared by filtration through a muslin cloth and subsequent centrifugation (5000× g). Finally, the crude extracts were freeze-dried (VirTis Advantage benchtop freeze dryer, Steroglass S.r.l., Perugia, Italy), stored at −20 • C, and reconstituted in demineralized water at the time of use.

Preparation of Crude Extracts (CE)
Flower heads from spontaneous O. tauricum Willd. plants were collected in July 2020 in the outer fringes of the Monti Sibillini National Park in the municipality of Visso (MC, Italy) according to Foligni et al. [25], while flower heads from cultivated plants were harvested in July 2020 from the experimental and didactic farm "Pasquale Rosati" (Università Politecnica delle Marche, Agugliano, AN, Italy) according to Zenobi et al. [31]. The crude extracts (CE) from spontaneous (st) and cultivated (ct) flower heads were prepared following the procedure outlined by Mozzon et al. [24] and illustrated in Figure  1. Briefly, the purple tubular flowers were separated from the receptacle and macerated in demineralized water (1:10 w/v) for 24 h at 4 °C; the resulting aqueous extracts were cleared by filtration through a muslin cloth and subsequent centrifugation (5000× g). Finally, the crude extracts were freeze-dried (VirTis Advantage benchtop freeze dryer, Steroglass S.r.l., Perugia, Italy), stored at −20 °C, and reconstituted in demineralized water at the time of use.

Cheesemaking and Sampling
Cheesemaking trials for manufacturing of Caciofiore (CF) were conducted in a family run dairy farm located in Pieve Torina (MC), in the Monti Sibillini National Park (Marche Region, Italy), following a traditional method (Supplementary Figure 1) and using raw Sopravissana breed ewe's milk, without the addition of starter cultures. The experimental plan for the cheesemaking trials is depicted in Figure 2. In the first production round (R1), conducted between April and July 2021, experimental cheeses (coded as CF_Ot_st) were curdled with the reconstituted CE_st, while in the second production round (R2), conducted between September and November 2021, experimental cheeses (coded as CF_Ot_ct) were curdled with the reconstituted CE_ct.

Cheesemaking and Sampling
Cheesemaking trials for manufacturing of Caciofiore (CF) were conducted in a family run dairy farm located in Pieve Torina (MC), in the Monti Sibillini National Park (Marche Region, Italy), following a traditional method (Supplementary Figure S1) and using raw Sopravissana breed ewe's milk, without the addition of starter cultures. The experimental plan for the cheesemaking trials is depicted in Figure 2. In the first production round (R1), conducted between April and July 2021, experimental cheeses (coded as CF_Ot_st) were curdled with the reconstituted CE_st, while in the second production round (R2), conducted between September and November 2021, experimental cheeses (coded as CF_Ot_ct) were curdled with the reconstituted CE_ct. Experimental plan of the cheesemaking trials. Two production rounds were run, the first (R1) between April and July 2021, and the second (R2) between September and November 2021; for each Round, trials were repeated on two different days (Batch 1, B1, and Batch 2, B2, respectively). On the same day, both control (CF_C) and experimental (CF_Ot_st in R1 and CF_Ot_ct in R2) Caciofiore cheeses were produced.
For each round, cheesemaking trials were repeated twice, on two different days, using two different batches of ewe's raw milk (B1 and B2, respectively). For each round, control cheeses (coded as CF_C) were also produced using a commercial vegetable coagulant (Dairen, Dairy and Food, Castel San Pietro Terme, BO, Italy) obtained from Figure 2. Experimental plan of the cheesemaking trials. Two production rounds were run, the first (R1) between April and July 2021, and the second (R2) between September and November 2021; for each Round, trials were repeated on two different days (Batch 1, B1, and Batch 2, B2, respectively). On the same day, both control (CF_C) and experimental (CF_Ot_st in R1 and CF_Ot_ct in R2) Caciofiore cheeses were produced. For each round, cheesemaking trials were repeated twice, on two different days, using two different batches of ewe's raw milk (B1 and B2, respectively). For each round, control cheeses (coded as CF_C) were also produced using a commercial vegetable coagulant (Dairen, Dairy and Food, Castel San Pietro Terme, BO, Italy) obtained from flowers of Cynara cardunculus L. (Cc). For each cheese manufacture, seven cheese wheels, each weighing about 250 g, were produced. Both control (CF_C) and experimental cheeses (CF_Ot_st and CF_Ot_ct) were ripened for 60 days at 12 • C and 70% of relative humidity (RH).

Physico-Chemical Measurements
The pH values of milk and Caciofiore cheese samples collected immediately after molding (day 0) as well as after 2, 5, 15, 30, 45, and 60 days of ripening were determined using a pHmeter Model 300 (Hanna Instruments, Padova, Italy) equipped with a solid electrode. Total titratable acidity (TTA) was determined on 10 g of sample homogenized in 90 mL of distilled water, titrated to a pH value of 8.3 with a 0.1 N solution of NaOH [36]. Percentage (%) TTA of lactic acid equivalents was calculated according to the following formula: where N is the normality of the titrant and 90 is the equivalent weight for lactic acid. Both pH and TTA were reported as mean ± standard deviation of two biological and two technical replicates.

Microbiological Analyses
Microbial counts were determined on commercial (Cc) and experimental coagulants (reconstituted CE_st and CE_ct), raw milk, and Caciofiore cheese samples; 10 g of each cheese sample was homogenized with 90 mL of sterile peptone water (bacteriological peptone, 1 g L −1 , Oxoid, Basingstoke, UK) using a stomacher apparatus (400 Circulator, International PBI, Milan, Italy) for 2 min at 260 rpm [37]. Thistle-based coagulants, raw milk, and cheese homogenates were serially diluted in sterile peptone water. Aliquots of decimal dilutions were inoculated in duplicate on the opportune growth media for the enumeration of the following microbial groups: (i) total mesophilic aerobes on Plate Count Agar (PCA) incubated at 30 • C for 48 h; (ii) presumptive lactobacilli on de Man, Rogosa, Sharp (MRS) agar incubated at 30 • C for 48 h; (iii) presumptive mesophilic and thermophilic lactococci on M17 agar incubated at 22 and 45 • C, respectively, for 48 h; (iv) Enterobacteriaceae on Violet Red Bile Glucose Agar (VRBGA) incubated at 37 • C for 24 h; (v) coliforms and Escherichia coli on Chromogenic Coliform Agar (CCA) incubated at 37 • C for 24 h; (vi) Pseudomonadaceae on Pseudomonas Agar Base (PAB) supplemented with Cetrimide, Fusidic Acid, and Cephaloridine (CFC), incubated at 30 • C for 48 h; (vii) yeasts and molds on Rose Bengal (RB) agar incubated at 25 • C for 72 h. The results of viable counts were expressed as the log of colony-forming units (cfu) per gram or mL of sample and reported as mean value ± standard deviation of two biological and two technical replicates.

DNA Extraction and Sequencing
The total microbial DNA was extracted from the cell pellets obtained by centrifugation of 1.5 mL of each biological replicate (dilutions 10 −1 ) using the E.Z.N.A. soil DNA kit (Omega Bio-tek, Norcross, GA, USA) . For each milk and cheese sample, the DNA extracts obtained from each biological replicate were pooled to reduce the inter-sample variability, as previously elucidated [38]. DNAs were quantified using the QUBIT dsDNA Assay kit (Life Technologies, Milan, Italy), and then diluted to 5 ng µL −1 to perform the metataxonomic amplicon sequencing. The amplification of the V3-V4 region of the 16S rRNA gene was carried out to determine the bacterial biota, in accordance with the conditions described by Klindworth et al. [39], whereas fungi were analyzed by amplifying the D1-D2 region of the 26S rRNA gene [40]. The amplicons were purified, tagged, and pooled following the Illumina standard protocol. Paired-end reads (2 × 250 bp) were sequenced by using a MiSeq platform (Illumina, San Diego, CA, USA) with v2 chemistry, according to the manufacturer's instructions.

Metataxonomic Analyses
QIIME2 software was used to import and analyze raw reads (FASTQ format) [41]. Barcode adapters and primers were first removed, and sequences were quality filtered applying DADA2 algorithm to obtain amplicon sequence variants (ASVs) [42]. The taxonomic assignment of the bacterial ASVs was executed using the Greengenes 16S rRNA database, whereas the manually build database described by Mota-Gutierrez et al. [43] was used for the taxonomic assignment of the fungal ASVs. The resulting ASVs with a relative frequency > 1.0% in at least two samples were retained. To confirm the taxonomic assignment of both bacterial and fungal ASVs, a double check on the Basic Local Alignment Search Tool (BLAST) was performed. The results of metataxonomic analysis were deposited in the Sequence Read Archive of the National Center for Biotechnology Information (NCBI) under the bioproject accession number PRJNA843917.

Statistical Analysis
The Tukey-Kramer's Honest Significant Difference (HSD) test (α = 0.05) was carried out to evaluate differences within thistle coagulants (Cc, CE_st, CE_ct), and cheese samples (CF_C vs CF_Ot_st and CF_C vs CF_Ot_ct) by one-way analysis of variance (ANOVA) using the software JMP ® Version 11.0.0 (SAS Institute Inc., Cary, NC, USA). The QIIME2 diversity plugin was used to calculate Alpha and Beta diversity indices for bacterial biota and mycobiota of thistle coagulants, raw milk, and Caciofiore cheese samples. Significant statistical differences (p-value < 0.05) were evaluated by the permutational multivariate analysis of variance (PERMANOVA). Principal component analysis (PCA) was performed on ASV relative frequencies using ade4 package in R v4.2.0. Differences in bacterial biota and mycobiota of Caciofiore cheese samples were evaluated by using the Kruskal-Wallis's test (p-value < 0.05); the ggplot2 and reshape2 packages in the R environment were applied to obtain boxplot showing the relative frequencies of differential ASVs.

pH, TTA, and Viable Counting
Raw ewe's milk used in the two cheesemaking rounds, R1 and R2, showed pH and TTA values of 6.74 ± 0.02 and 6.56 ± 0.02, respectively, for the first parameter, and of 0.15 ± 0.01 and 0.16 ± 0.02%, respectively, for the latter. The trends of pH and TTA in the cheeses during ripening are shown in Figure 3. In both rounds, a rapid drop in pH was observed during the first 5 days in all the samples; the minimum value was measured after 15 days, followed by a slight increase in the next ripening phase (30-60 days). An opposite trend was observed for TTA, but with a slower increase (maximum values reached between 30 and 45 days). No significant differences were observed between control and experimental cheese samples collected from the same round and the same sampling time (CF_C vs CF_Ot_st and CF_C vs CF_Ot_ct, respectively). Table 1 shows the viable counts of milk coagulants. The commercial coagulant (Cc) showed counts < 1 Log cfu g −1 for all the microbial groups, except for yeast (6.36 ± 0.02 Log cfu g −1 ). The load of presumptive lactobacilli and Escherichia coli was <1 Log cfu g −1 for all coagulants. The two experimental coagulants from flowers of spontaneous (CE_st) and cultivated (CE_ct) O. tauricum were characterized by similar counts of molds, whereas significant differences of about 1-1.5 order of magnitude were observed in the loads of total mesophilic aerobes, presumptive lactococci, and presumptive thermophilic cocci, with higher values in CE_st compared to CE_ct. Even greater differences emerged in the viable counts of the remaining microbial groups, with CE_st showing counts of 3.04 ± 0.06, 3.51 ± 0.16, 3.64 ± 0.11, and 2.90 ± 0.00 Log cfu g −1 for Enterobacteriaceae, coliforms, Pseudomonadaceae, and yeasts respectively, and CE_ct being characterized by counts always <1 Log cfu g −1 .
± 0.01 and 0.16 ± 0.02 %, respectively, for the latter. The trends of pH and TTA in the cheeses during ripening are shown in Figure 3. In both rounds, a rapid drop in pH was observed during the first 5 days in all the samples; the minimum value was measured after 15 days, followed by a slight increase in the next ripening phase (30-60 days). An opposite trend was observed for TTA, but with a slower increase (maximum values reached between 30 and 45 days). No significant differences were observed between control and experimental cheese samples collected from the same round and the same sampling time (CF_C vs CF_Ot_st and CF_C vs CF_Ot_ct, respectively).  Table 1 shows the viable counts of milk coagulants. The commercial coagulant (Cc) showed counts < 1 Log cfu g −1 for all the microbial groups, except for yeast (6.36 ± 0.02 Log cfu g −1 ). The load of presumptive lactobacilli and Escherichia coli was <1 Log cfu g −1 for all coagulants. The two experimental coagulants from flowers of spontaneous (CE_st) and cultivated (CE_ct) O. tauricum were characterized by similar counts of molds, whereas significant differences of about 1-1.5 order of magnitude were observed in the loads of total mesophilic aerobes, presumptive lactococci, and presumptive thermophilic cocci, with higher values in CE_st compared to CE_ct. Even greater differences emerged in the viable counts of the remaining microbial groups, with CE_st showing counts of 3.04 ± 0.06, 3.51 ± 0.16, 3.64 ± 0.11, and 2.90 ± 0.00 Log cfu g −1 for Enterobacteriaceae, coliforms, Pseudomonadaceae, and yeasts respectively, and CE_ct being characterized by counts always <1 Log cfu g −1 .   Values are expressed as Log cfu g −1 ± standard deviation. Within each row, overall means with different superscript letters are significantly different (p < 0.05).
The microbial viable counts of raw milk and cheese samples collected during ripening are shown in Tables 2 and 3, for R1 and R2 samples, respectively. As a general trend, a high similarity was seen in the viable counts of control and experimental cheeses (CF_C and CF_Ot_st in R1, CF_C and CF_Ot_ct in R2, respectively). In both rounds, a similar evolution over time was observed for total mesophilic aerobes, presumptive lactococci, and presumptive thermophilic cocci, with a significant increase, of about 3 Log, in the first two days, and a further increase to reach the maximum load between 5 and 15 days, followed by a slight reduction in the last stage of maturation. The same behavior was also observed for lactobacilli in the second round (R2), while in the first round (R1), the growth of this microbial group was slower, and the maximum load was reached after 30 days of ripening in both control (CF_C) and experimental (CF_Ot_st) cheeses. As for Enterobacteriaceae, no significant differences were found between control and experimental cheeses: the maximum loads, reached after 5 days of ripening (counts in the range of 7.35-7.99 Log cfu g −1 ), were followed by a reduction, up to 4 and 2 Logs in R1 and R2, respectively; a similar trend was also observed for coliforms. Escherichia coli was detected in both milk batches (2.18 ± 0.00 vs. 2.53 ± 0.83 Log cfu g −1 in R1_M and R2_M, respectively). Regarding R1, the counts of E. coli in cheese samples remained almost stable during maturation, with mean values of 3.82 and 4.07 Log cfu g −1 in CF_C and CF_Ot_st, respectively; in R2 cheeses, maximum counts of 6.75 ± 0.22 (CF_C) and 7.16 ± 0.21 (CF_Ot_ct) Log cfu g −1 were found after 5 days, followed by a slow decrease. For R2, after 30 days of ripening, a lower Pseudomonadaceae count was detected in cheeses coagulated with CE_st than those coagulated with commercial rennet (5.09 ± 0.69 vs. 6.51 ± 0.20 Log cfu g −1 , respectively), although at the end of ripening (60 days), no more differences could be evidenced. For viable counts of molds, the highest loads were found between 5 and 15 days of ripening in both cheesemaking rounds; moreover, no significant differences were observed between control and experimental cheeses, except for CF_C and CF_Ot_st analyzed after 45 days of ripening. Finally, as for yeasts, in both rounds the maximum loads were found at the fifth day of maturation, and significant differences between control and experimental cheeses were observed only in cheeses at 15 days of maturation in R2 (6.16 ± 0.63 vs. 5.09 ± 0.19 Log cfu g −1 in CF_C and CF_Ot_ct, respectively) and 60 days in R1 (3.54 ± 0.10 vs. 4.80 ± 0.79 Log cfu g −1 in CF_C and CF_Ot_st, respectively).
In the first round (R1), the comparison between control and experimental cheeses (CF_C vs CF_Ot_st) showed no significant differences, as depicted in the PCA plot (Figure 5, panel (a)). In fact, all cheeses were characterized by a substantial increase of Lactococcus lactis over time, with initial relative frequencies of 3.71 (R1_CF_C_t0) and 3.41% (R1_CF_Ot_st_t0), and relative frequencies at the end of the ripening period of 62.54 (R1_CF_C_t60) and 48.15% (R1_CF_Ot_st_t60), respectively. After five days of ripening, Lactobacillus spp. and Latilactobacillus sakei displayed a progressive increase, reaching up to 7.77 and 4.37% of the relative frequency in R1_CF_C_t60, and 6.66 and 7.61% of the relative frequency in R1_CF_Ot_st_t60, respectively. As expected, Staphylococcus aureus and Acinetobacter johnsonii were detected only at the beginning of the cheesemaking processes. Citrobacter, Serratia, and Enterobacteriaceae were stably detected at all sampling times in both production rounds, with overall mean relative frequencies of 24.08, 4.65 and 2.32%, respectively.   In the first round (R1), the comparison between control and experimental cheeses (CF_C vs CF_Ot_st) showed no significant differences, as depicted in the PCA plot ( Figure  5, panel (a)). In fact, all cheeses were characterized by a substantial increase of Lactococcus lactis over time, with initial relative frequencies of 3.71 (R1_CF_C_t0) and 3.41% (R1_CF_Ot_st_t0), and relative frequencies at the end of the ripening period of 62.54 Even in the second round (R2), no significant differences emerged by comparing control and experimental cheeses (CF_C vs CF_Ot_ct), as shown in the PCA plot ( Figure 5,  panel (a)). R2_CF_C_t0 was characterized by a low bacterial diversity, with Lactococcus lactis being the sole species detected among bacteria, whereas in R2_CF_Ot_ct_t0 Rummeliibacillus ssp. and Macrococcus caseolyticus were detected, with 38.68 and 13.79% of the relative frequency, respectively. In both cheese manufactures, Lactococcus lactis drove the fermentation processes over time, reaching final relative frequencies of 54.38 and 56.57%, in R2_CF_C_t60 and R2_CF_Ot_ct_t60, respectively. Since day 5, a gradual increase was also seen for Lactobacillus spp., which reached final relative frequencies of 20.85 and 20.03% in R2_CF_C_t60 and R2_CF_Ot_ct_t60, respectively. Staphylococcus spp. was among the preva-lent taxa in both cheese manufactures during the early maturation, with relative frequencies of 8.06 and 7.28% in R2_CF_C_t5 and R2_CF_Ot_ct_t5, respectively; as expected, this genus decreased below 1% at the end of ripening. Enterobacteriaceae, Citrobacter, Enterococcus, and Weissella were stably detected at all sampling times irrespectively of the milk coagulant used, with overall mean relative frequencies of 14.75, 3.45, 2.02, and 1.87%, respectively.
Microorganisms 2023, 11, x FOR PEER REVIEW 12 of 25 (R1_CF_C_t60) and 48.15% (R1_CF_Ot_st_t60), respectively. After five days of ripening, Lactobacillus spp. and Latilactobacillus sakei displayed a progressive increase, reaching up to 7.77 and 4.37% of the relative frequency in R1_CF_C_t60, and 6.66 and 7.61% of the relative frequency in R1_CF_Ot_st_t60, respectively. As expected, Staphylococcus aureus and Acinetobacter johnsonii were detected only at the beginning of the cheesemaking processes. Citrobacter, Serratia, and Enterobacteriaceae were stably detected at all sampling times in both production rounds, with overall mean relative frequencies of 24.08, 4.65 and 2.32%, respectively.  Globally, the Kruskal-Wallis test evidenced thirteen statistically different bacterial groups among the cheese manufactures, and their boxplots are showed in Figure 6. Citrobacter, Lactococcus, Latilactobacillus sakei, Leuconostoc mesenteroides, Serratia, and Staphylococcus aureus were mainly associated with cheeses produced in round 1. Conversely, Corynebacterium, Enterobacteriaceae, Enterococcus, Lactobacillus delbrueckii, Staphylococcus, Staphylococcus sciuri, and Weissella were mainly associated with cheeses produced in round 2.  Concerning the metataxonomic analysis of the mycobiota, 1,871,974 high quality reads were used, with an average of 53,485 sequence per sample and a coverage > 99%. The Alpha diversity analysis displayed different levels of complexity depending on the type of sample. In more detail, the fungal composition of raw milk showed a higher complexity regarding cheeses, whereas the thistle-based coagulants showed no significant differences if compared to raw milk and cheeses sampled at different maturation stages (p-value < 0.05) (Supplementary Figure S3). The results from Illumina DNA sequencing, expressed as relative frequency of fungal ASVs, are shown in Figure 7. The mycobiota of the commercial coagulant (Cc) was dominated by Debaryomyces hansenii, attesting at 99.51% of the relative frequency, whereas CE_st and CE_ct shared the presence of Saccharomyces cerevisiae (7.22 and 55.30% of the relative frequency, respectively), Aerobasidium proteae (15.24 and 2.96% of the relative frequency, respectively), Aspergillus (5.82 and 3.36%, respectively), Alternaria (4.85 and 0.96%, respectively), and Debaryomyces hansenii (2.36 and 0.46%, respectively). In addition, CE_st included Eremothecium coryli, Kluyveromyces marxianus, Metschnikowia, and Penicillium, with relative frequency of 13.73, 4.64, 3.01, and 1.84%, respectively. As for raw milk batches R1_M and R2_M, the common presence of Geotrichum silvicola (24.82 and 45.22% of the relative frequency, respectively), Cladosporium (10.21 and 21.93%, respectively), Yarrowia lipolytica (20.26 and 4.24%, respectively), Candida parapsilosis (3.61 and 10.01%, respectively), and Wickerhamiella pararugosa (4.42 and 2.44%, respectively) was highlighted. Saccharomyces cerevisiae, Kurtzmaniella zeylanoides, and Didymella were also detected in R1_M, with 12.35, 7.01, and 3.49% of the relative frequency, respectively, while Kluyveromyces marxianus and Pichia cactophila were identified in R2_M, with 2.87 and 2.64% of the relative frequency, respectively.
Cheeses manufactured in round 1 showed no significant differences related to the type of coagulant, as depicted in the PCA plot ( Figure 5, panel (b)). Debaryomyces hansenii was predominant in both cheese manufactures, reaching up to 85.96 (R1_CF_C_t60) and 73.14% (R1_CF_Ot_st_t60) of the relative frequency at the end of ripening. Penicillium spp. was detected at variable levels in almost all the cheese samples, with an overall mean relative frequency of 19.94%. A decreasing trend was seen for Geotrichum silvicola in both cheese manufactures; Yarrowia lipolytica was among the prevalent taxa exclusively at the beginning of the monitoring process, with relative frequencies of 8.33 and 5.42% in R1_CF_C_t0 and R1_CF_Ot_st_t0, respectively. Similar results were found for the cheeses manufactured in the second round (R2); in both control and experimental cheese manufactures, Debaryomyces hansenii increased with time up to 49.58 (R2_CF_C_t60) and 73.97% (R2_CF_Ot_ct_t60) of the relative frequency. An analogous trend was observed for Penicillium spp., which after 60 days of ripening reached 22.18 (R2_CF_C_t60) and 8.72% (R2_CF_Ot_ct_t60) of the relative frequency. On the contrary, Geotrichum silvicola and Candida parapsilosis decreased over time. Like Y. lipolytica, even Cladosporium spp. was identified among the main taxa just after molding, attesting at 13.24 (R2_CF_C_t0) and 2.53% (R2_CF_Ot_ct_t0) of the relative frequency. Pichia cactophila was stably detected in all cheese samples, with an overall mean relative frequency of 1.64%. Again, as shown in the PCA plot ( Figure 5, panel  (b)), no significant differences emerged by comparing control (R2_CF_C) and experimental (R2_CF_Ot_ct) cheeses.
Globally, the Kruskal-Wallis test evidenced four statistically different fungal groups among cheese manufactures, and their boxplots are showed in Figure 8. Aerobasidium proteae was mainly associated with Caciofiore cheeses produced in round 1, whereas Kluyveromyces marxianus, Pichia cactophila, and Pichia kudriavzevii with Caciofiore cheeses produced in round 2.  the PCA plot ( Figure 5, panel (b)), no significant differences emerged by comparing control (R2_CF_C) and experimental (R2_CF_Ot_ct) cheeses. Globally, the Kruskal-Wallis test evidenced four statistically different fungal groups among cheese manufactures, and their boxplots are showed in Figure 8. Aerobasidium proteae was mainly associated with Caciofiore cheeses produced in round 1, whereas Kluyveromyces marxianus, Pichia cactophila, and Pichia kudriavzevii with Caciofiore cheeses produced in round 2.

Discussion
To the author's knowledge, this is the very first report describing the exploitation of O. tauricum Willd. for cheese manufacturing. The present investigation is also among the very few ones till date carried out to evaluate the impact of thistle rennet on cheese microbial dynamics during fermentation and ripening. From a microbiological point of

Discussion
To the author's knowledge, this is the very first report describing the exploitation of O. tauricum Willd. for cheese manufacturing. The present investigation is also among the very few ones till date carried out to evaluate the impact of thistle rennet on cheese microbial dynamics during fermentation and ripening. From a microbiological point of view, as a main ingredient of thistle curdled cheeses, thistle rennet must not represent a source of undesirable microorganisms, with detrimental effects on cheeses (e.g., off-odors or offflavors, undesired changes in texture, appearance, or palatability, etc.) or consumers' health (e.g., foodborne disease, food poisoning). By contrast, the occurrence of pro-technological or functional microbes in thistle rennet might be advantageous for obtaining high-quality cheeses, especially when primary or secondary starters are not used, as in the cheeses herein manufactured.

Viable Counting
Except for presumptive lactobacilli, viable counts of the coagulant obtained from spontaneous O. tauricum populations (CE_st) were consistent with those previously reported for crude extracts of other thistle species with a documented milk-coagulating potential [33][34][35]44]. As a general trend, counts of CE_ct were significantly lower than those of the extract from spontaneous O. tauricum (CE_st). For Enterobacteriaceae, coliforms, and Pseudomonadaceae, viable counts of CE_ct < 1 Log cfu g −1 suggest a higher hygienic quality of these extracts compared to those prepared from spontaneous O. tauricum (CE_st); this finding might be tentatively ascribed to the differences in the management of cultivated fields than pastures or fallow lands, more often accessed or grazed by wild animals or even livestock, whose feces are known to harbor the aforementioned microorganisms [45,46]. In this regard, recent comparative metagenomic analyses also highlighted the strong influence of farming practices on taxonomic and functional diversity of phyllosphere microbes [47].
However, both control and experimental cheeses herein manufactured were dominated by lactic acid bacteria, which exhibited an intense growth in the early stage of maturation (0-15 days), associated with a significant pH drop, followed by a stabilization phase (30-60 days). Viable counts of presumptive lactococci and thermophilic cocci in raw milk and cheeses during their early maturation were comparable to those overall collected by Cardinali et al. [34] in the same specialty cheese manufactured with a crude extract of Carlina acanthifolia All. subsp. acanthifolia.
Enterobacteriaceae, coliforms, and E. coli were enumerated in raw milk used in R1 and R2, respectively. In cheeses, after five days of maturation, they decreased, but without disappearing completely up until the end of maturation (60 days), when higher counts were seen in R2 than R1. Members of the Enterobacteriaceae family have previously been detected in raw milk and raw milk cheeses, especially in the early cheese maturation [48][49][50]. The presence of these microorganisms in cheese is controversial for safety and technological issues [51]: in fact, though some species are known to produce aromatic volatile compounds [52], some others include strains with pathogenic traits [53].
Regarding eumycetes, the differences found in the viable counts of yeasts and molds in the three assayed thistle coagulants did not affect the fungal dynamics during cheese ripening. Similar evidence emerged in the study carried out by Cardinali et al. [35].

Illumina Sequencing
Over the past decades, to overcome the limitations of conventional culture-dependent analyses, DNA-based techniques have gained increasing acceptance, allowing for better investigation of the microbiota of foods [54][55][56]. Among these techniques, next-generation sequencing has previously been exploited for the investigation of complex microbial communities of thistle-curdled cheeses [13][14][15]. In the present study, the application of Illumina sequencing on nucleic acids extracted from thistle-based coagulants, raw milk, and cheese samples provided a comprehensive insight on the dynamics of bacteria and fungi, from the raw ingredients to mature cheeses.
Cc, CE_st, and CE_ct were mainly dominated by spoilage and environmental microorganisms, whereas lactic acid bacteria were only found at low relative frequencies, with Lactococcus lactis being detected in Cc, CE_st, and CE_ct, Latilactobacillus sakei in Cc and CE_st, and Lactobacillus spp. in Cc and CE_ct, respectively.
Acinetobacter johnsonii, Staphylococcus spp., Geotrichum silvicola, Cladosporium spp., Yarrowia lipolytica, Candida parapsilosis, and Wickerhamiella pararugosa were among the prevalent taxa in both batches of Sopravissana ewe's raw milk. The occurrence of Acinetobacter johnsonii, Staphylococcus spp., and Candida parapsilosis has previously been documented in raw milk of the same sheep breed [34,35]. Geotrichum species, mainly G. candidum, are also commonly associated with milk and cheese [57][58][59]; although recognized as synonym of Galactomyces candidus [60] due to phylogenetic affinity, to the best of the authors' knowledge this is the first study reporting the identification of Geotrichum silvicola in raw ewe's milk and cheese. Great interest is addressed to the yeast Yarrowia lipolytica for its proteolytic and lipolytic activities, and the consequent aroma development during cheese ripening [61][62][63]. To date, Y. lipolytica has frequently been identified in cheeses [64], whereas its occurrence in milk has more poorly been documented [65][66][67][68]. However, Corbo et al. [69] remarked the predominance of Y. lipolytica in ewe's milk, which is characterized by a higher fat content than cow's and goat's milk. As for Wickerhamiella pararugosa (synonym Candida pararugosa), this yeast has recently been found among the minority species in various dairy products [13,70,71] and milk [72]. Air is considered the main contamination source of Cladosporium spp. [73,74], a mold commonly identified in raw milk [75][76][77].
Several minority bacterial and fungal taxa were also identified in the two raw milk batches by Illumina DNA sequencing; though Latilactobacillus sakei, Enterococcus spp., Macrococcus caseolyticus, Lactobacillus spp., and Kluyveromyces marxianus were detected in both R1_M and R2_M, these batches differed significantly in number and relative abundance of the detected taxa, with R1_M harboring Kurtzmaniella zeylanoides, and Saccharomyces cerevisiae and R2_M being characterized by the occurrence of Lactobacillus delbrueckii, Weissella spp., and Pichia spp. As previously suggested by other authors, the differences emerged between the two raw milk samples, collected in spring (R1_M) and late summer (R2_M) 2021, respectively, might be likely ascribed to the seasonality [78][79][80].
However, in both production rounds, Lactococcus lactis among bacteria and Debaryomyces hansenii among eumycetes dominated the microbial dynamics during fermentation and maturation of both control and experimental cheeses. The latter finding agrees with those of previous studies carried out in typical Mediterranean sheep and goat cheeses [81][82][83][84][85][86][87]. L. lactis is a homofermentative starter microorganism known for its ability to lower the pH and drive cheese fermentation [88]. D. hansenii is a yeast species that plays a key role during cheese ripening through the production of proteolytic and lipolytic enzymes acting on milk proteins and fat [89].
Furthermore, concerning bacterial dynamics, in both control and experimental cheeses, an increase in the relative frequency of Lactobacillus spp. was detected starting from t 15 in R1, and t 5 in R2, respectively, until the last sampling time (t 60 ), thus showing the ability of this genus to adapt to the cheese habitat, and suggesting its important role during cheese fermentation and ripening [90]. A similar behavior has been also observed for Latilactobacillus sakei in R1, and Weissella spp. in R2. All the cited lactic acid bacteria have previously been found as part of the bacterial biota of raw milk cheeses [13,[91][92][93]. As for fungal evolution, in R1, most of the predominant taxa in milk and early ripened cheeses were found as minority components at the end of maturation. In R2, Candida parapsilosis and Geotrichum silvicola were among the dominant species until the end of maturation despite a significant reduction over time. In the same production round, Kluyveromyces marxianus was also detected; the occurrence of this yeast species in ewe's milk cheeses is well documented, where it contributes to maturation and aroma formation owing to its ability to metabolize lactose, proteins, and fat [81,83,84,94].
For the stable detection of the genus Penicillium in all cheeses herein analyzed, independent of the production round and the rennet used, it might be likely explained by a contamination from the air or the ripening environment, as previously been suggested for other cheeses [74].
Viable counts also showed a stable occurrence of presumptive Pseudomonadaceae in samples collected from both production rounds; however, neither Pseudomonas spp. nor Pseudomonas fragi were detected in any cheese sample by Illumina sequencing after 30 days of maturation, thus once more underlying the importance of using a combined analytical approach, based on both culture-dependent and -independent analyses.

Conclusions
The combination of conventional microbiological analyses and Illumina DNA sequencing applied on experimental and commercial thistle-based coagulants, raw milk, and cheeses sampled at different ripening times allowed an in-depth investigation of the microbial dynamics occurring during fermentation and ripening of Caciofiore cheese manufactured with raw ewe's milk curdled with thistle rennet. Despite the differences that emerged between the three milk coagulants, the results herein collected revealed that the microbiological quality of raw milk is among the main factors affecting the evolution of the bacterial and fungal biota during cheese fermentation and ripening. By contrast, in accordance with the results previously collected by Aquilanti et al. [33] and Cardinali et al. [34,35], a neglectable impact of thistle rennet on cheese microbial dynamics emerged, with the rennet microbiota being mainly dominated by adventitious microorganisms with neither dairy attitude nor detrimental traits. In addition, from a microbiological point of view, at the end of ripening, mature cheeses manufactured with the experimental coagulants were equivalent to those produced with a commercially available thistle rennet obtained from C. cardunculus, thus supporting the potential use of the newly explored rennet from O. tauricum for manufacturing of ewe's milk cheeses. In this regard, the evidence in present study supports the large-scale cultivation of O. tauricum for production of a new alternative thistle rennet to be proposed to dairy farmers for manufacturing of Mediterranean ovine cheeses.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/microorganisms11010219/s1. Figure S1: Flow chart of the cheesemaking process of Caciofiore. Figure S2: Alpha diversity plots of Pielou's evenness and Shannon indices for bacterial biota based on the sample matrices (Caciofiore cheeses, raw milks, and thistles coagulants). Figure S3: Alpha diversity plots of Pielou's evenness and Shannon indices for mycobiota based on the sample matrices (Caciofiore cheeses, raw milks, and thistle coagulants).