All procedures were conducted according to the European guidelines for the care and use of animals in research (Directive 2010/63/EU). The experimental protocol was approved by the Ethical Committee of the Department of Veterinary Science of the University of Messina (code 041/2020).

The experimental period lasted from February to June 2021 and was performed on a commercial dairy farm located in the Sicilian region (Ragusa, Italy) at an altitude of 520 m above sea level. One hundred and twenty healthy multiparous dairy Friesian cows were divided into two groups (60 animals each), called Control (Ctr) and Experimental (Exp), homogenous for the Body Condition Score (3 ± 0.5), distance from calving (90–120 days), and milk production (25 ± 3 kg/day), fed a diet (20 kg of Dray Matter (DM)/head per day as Total Mixed Ration) with concentrate (Crude Protein: 170 g/kg of DM; Crude Cellulose 69 g/kg; Crude oils and fats: 4.1 g/kg of DM; Ashes 76 g/kg of DM) and meadow hay (Crude Protein: 110.9 g/kg of DM; Ether Extract: 25.0 g/kg of DM; Neutral Detergent Fiber: 521.9 g/kg of DM). The Control group (Ctr) received a concentrate without any olive cake supplements, whereas the Experimental group (Exp) received a concentrate supplemented with the olive cake as 8% on a DM as reported in Table 1. A 3-week period of adaptation to diets was assured before collecting samples. The by-product used was represented by the residue of the mechanical pressing of olives with the two-stage process, to produce extra virgin olive oil. In addition, the by-product contained about 5% of leaf residue with which the olives normally arrive at the mill. This olive cake was subsequently pitted by centrifugation and dried in the open air. Its moisture content went from 60% to a value of 10%. Its nutritional characteristics are reported in Table 2.

Milk samples, obtained from control and experimental groups, were monthly collected in distinct tanks after mechanical milking performed in the evening (3 p.m.) and in the morning (7 a.m.). Samples were transferred, under refrigerated conditions (4 ± 2°C), to the Natura & Qualità (Soc. Agr. ARL) dairy company (Ragusa, Italy) and subjected to both physicochemical and microbiological analyses before cheese-making. Milk samples were analyzed for the fat, protein, casein, and lactose content by using a Fourier transform infrared (Milkoscan FT2, FOSS, Hillerød, Denmark) and the somatic cell count was determined using a FossomaticTM FC (FOSS, Hillerød, Denmark). Microbiological count was performed to assess the presence of Enterobacteriaceae, Escherichia coli and coliforms, total mesophilic bacteria, and Listeria monocytogenes, using the following agar media and growth conditions: Violet Red Bile Agar (VRBA) anaerobically incubated at 37°C and 45°C for 24–48 h, for the count of Enterobacteria and fecal coliforms, respectively; Plate Count Agar (PCA), aerobically incubated at 30°C for 48–72 h, for total mesophilic aerobic bacteria count; and Chromatic ECX gluc agar aerobically incubated at 37°C for 18–24 h for E. coli. In addition, according to the International Organization for Standardization (ISO), a two-stage enrichment method for L. monocytogenes (ISO (International Organization for Standardization) (ISO 11290-1:2017), 2017a) and Salmonella spp. detection (ISO (International Organization for Standardization) (ISO 6579-1:2017), 2017b) was used. All media were purchased from Liofilchem (Roseto degli Abruzzi, Italy).

Cheese samples (weight 500 g) were obtained following the flowchart as reported in Figure 1. Provola cheese samples, from each group, were monthly collected, from March to July, vacuum-packed, and transported under refrigerated conditions (4 ± 2°C) to the Department of Veterinary Sciences of the University of Messina for nutritional analysis, and to the Department of Agriculture, Food and Environment, University of Catania for microbiological and 16S rRNA gene meta-taxonomic analysis.

One-day-old Provola cheese samples (100 g) were analyzed for chemical composition by FoodScan™ Dairy Analyser (FOSS, Italy). The fatty acid methyl esters (FAME) of cheese were analyzed by a GC–FID. For the chromatographic analyses, on 15 mg of lipids of cheese, the fatty acids methyl esters were prepared by using a solution of sulfuric acid/methanol (1:9, v/v) and submitted to a HRGC analysis. A gas chromatograph with a FID detector (Agilent Technologies 6890N, Palo Alto, CA, United States) equipped with a SP-2560 fused silica capillary column (100 m × 0.25 mm i.d. × 0.2 μm film thickness, Supelco, Inc., Bellefonte, PA, United States) was used. Helium was the carrier gas, 1 mL/min. The volume of the injection was 1 μL and the split ratio 1:50. The column temperature was programmed: an initial isotherm of 70°C (2 min), an increment of 15°C/min to 155°C (held for 25 min), then an increment of 3°C/min and a final isotherm of 215°C (8 min), following the procedure described by Tudisco et al. (2015). The fatty acids were identified by comparing the relative retention times of FAME peaks from samples containing standards from Supelco (mix 37 FAMEs, Supelco, Bellefonte, PA, United States). Chromatogram peak areas were acquired and calculated using a ChemStation software (Agilent, Santa Clara, CA, United States). The concentration of each fatty acid was expressed as g/100 g, considering 100 g as the summation of the areas of all the FAME identified. For each sample, the chromatographic analysis was repeated three times. The total polyphenol content was determined spectrophotometrically using the Folin–Ciocalteu’s reagent according to the method of Shetty et al. (1995). Approximately 50 mg of homogenized Provola was weighed, and 2.5 mL of 95% ethanol was added and kept at 0°C for 48 h. Around 1 mL of the supernatant was transferred into a tube and mixed with ethanol (95%) and ultrapure water. Folin–Ciocalteu’s reagent (50%) and Na₂CO₃ (5%) were added. Absorbances were read at 760 nm with 95% ethanol as analytical blank, using a UV–visible spectrophotometer (UV-2401PC Shimadzu). A calibration curve was constructed using appropriate dilutions of a gallic acid standard solution (95% in ethanol). Determinations were performed in triplicate together with blank solutions.

Microbiological analysis of 1-day old control and experimental cheese samples, monthly collected from March to July, was performed according to Commission Regulation (EC) No 2073/2005 of 15 November 2005 On microbiological criteria for foodstuffs (n.d.) by sampling plans with five random collected samples. In detail, from each cheese sample, both core and surface sections (25 g) were sampled and treated as previously described (Pino et al., 2018). Microbiological count was performed according to Pino et al. (2021) and Randazzo et al. (2021) using the following agar media and conditions: Kanamycin Aesculin Azide (KAA) Agar, aerobically incubated at 37°C for 24–48 h, for enterococci count; Mannitol Salt Agar (MSA), incubated at 32°C for 48 h, for staphylococci count; M17, supplemented with 5 g/L of lactose, incubated at 30°C for 24–48 h, for lactococci; de Man Rogosa and Sharpe agar (MRS), adjusted to pH 5.4, incubated at 32°C for 72 h under microaerophilic condition, for lactobacilli; and Sabouraud Dextrose Agar (SDA), supplemented with chloramphenicol (0.05 g/L) incubated at 25°C for 3–5 days, for yeasts and molds count. VRBA and PCA were used as previously described for Enterobacteria, fecal coliforms, and for total mesophilic aerobic bacteria count. Listeria monocytogenes (ISO (International Organization for Standardization) (ISO 11290-1:2017), 2017a) was enumerated as previously reported. All media were purchased from Oxoid (Basingstoke, United Kingdom). Moreover, according to the International Organization for Standardization (ISO), a two-stage enrichment method for Salmonella spp. detection (ISO (International Organization for Standardization) (ISO 6579-1:2017), 2017b) was used.

Control and experimental cheeses were subjected to total bacterial DNA isolation using the DNeasy^® mericon^® Food kit (Qiagen, Germany) following the manufacturer’s instructions. DNA concentration was determined using the fluorimeter Qubit 4.0 (Invitrogen, Carlsbad, CA, United States) before storing at −20°C until use.

Results of microbiological analysis were calculated as mean values of three determinations and standard deviations and were expressed as log cfu/g.

The DNA isolated from cheese samples was subjected to 16S rRNA gene meta-taxonomic analysis. In detail, the V3 region of the 16S rRNA gene sequence was amplified using the Probio_Uni e/Probio_Rev primer pair (Milani et al., 2013) and subjected to MiSeq (Illumina) by GenProbio srl (Parma, Italy). Bioinformatics analysis was performed as reported by Vaccalluzzo et al. (2022). 16S rRNA raw data from the study were deposited at NCBI Sequence Read Archive (SRA¹) under accession code PRJNA 903126.

Four grams of grated cheese were added with 10 μL of internal standard solution (2-octanol, at 10 ppm), placed into 20 mL glass vials, and sealed with polytetrafluoroethylene-coated silicone rubber septa (20 mm diameter; Supelco, Bellefonte, PA, United States). To obtain the best extraction efficiency, the micro-extraction procedure was performed as described by Salum et al. (2017), with slight modifications. After sample equilibration (10 min at 54.75°C), a conditioned 50/30 μm DVB/CAR/PDMS fiber (Supelco, Bellefonte, PA, United States) was exposed for 60 min. The temperature was kept constant during analysis, and the vials were maintained on a heater plate (CTC Analytics, Zwingen, Switzerland) of a CombiPAL system injector Autosampler (CTC Analytics). The extracted VOC were desorbed in splitless mode (3 min at 220°C) and analyzed through a Clarus 680 (Perkin Elmer) gas-chromatography (GC) system equipped with a capillary Rtx-WAX column (30 m × 0.25 mm i.d., 0.25 μm film thickness; Restek, Bellefonte, PA, United States). The column temperature was set initially at 35°C for 8 min, then increased to 60°C at 4°C min⁻¹, to 160°C at 6°C min⁻¹, and finally to 200°C at 20°C min⁻¹ and held for 15 min. Helium was used as the carrier gas at flow rate of 1 mL min⁻¹. A single quadrupole mass spectrometer (MS) Clarus SQ 8C (Perkin Elmer) was coupled to the GC system. The source and transfer line temperatures were kept at 250 and 230°C, respectively. Electron ionization masses were recorded at 70 eV in the mass-to-charge ratio (m/z) interval 34–350 (Randazzo et al., 2021). Each chromatogram was analyzed for peak identification using the NIST (National Institute of Standard and Technology) 2008 library. A peak area threshold of 1,000,000 and at least 85% probability of match were used for dentification, followed by visual inspection of the fragment patterns when required. The concentrations of VOC (calculated on internal standard base) were expressed as mg kg⁻¹.

With the aim of searching for statistically significant changes in alpha and beta diversity estimates and starting from Bray-Curtis, Jaccard, Weighted UniFrac, Unweighted UniFrac distance matrices, QIIME II nested plugins were used. The PERMANOVA test was computed in same software environment.

Taxa abundance in the two compared groups was inspected by mean of a compositionality-aware test, i.e., the Analysis of Composition of Microbiomes test (ANCOM), but also with the White’s non-parametric test corrected for multiple tests. In the case of White’s test, results have been shown in terms of boxplot.

For VOC statistical analyses, we stratified samples based on fed group and sampled month by make use of Partial least square discriminant analysis (PLSDA) and Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA).

Physicochemical composition and microbiological parameters of milk samples are reported in Supplementary Tables 1, 2, respectively. Cheese composition in terms of physicochemical parameters in both control and experimental samples, collected from March to July, is reported in Supplementary Table 3. The experimental cheese showed a higher content of polyphenols, proteins, n-3, n-6, PUFA, C18:0, C18:1n-9, C18:2n-6 cis, and C18:2n-6 trans, than control ones, whereas the concentrations of C14:0, C16:0, and C20:3n-3 were higher in the control cheeses. Physicochemical cheese detected variables were inspected for their differential abundances in the two compared groups in a Wilcoxon rank-sum test. Table 3 reports only parameters that were significantly different (FDR < 0.05).

Results of microbiological analysis carried out on both control and experimental cheese samples are shown in Table 4. Overall, Salmonella spp., L. monocytogenes, E. coli, fecal coliforms, and molds were never detected in the analyzed samples, except for the control cheese collected in July, where coliforms reached the value of 2.00 log cfu/g (data not shown). Lactic acid bacteria (LAB) achieved an average value of about 6.84 log cfu/g in control samples and of 7.19 log cfu/g in experimental cheeses. In detail, control samples reached the highest values (7.99 log cfu/g) in April and the lowest (6.04 log cfu/g) in June, while the experimental ones achieved the highest (7.96 log cfu/g) in April and the lowest (6.55 log cfu/g) in July (Table 4). The lactococci group showed the same average value (about 8.00 log cfu/g) in both samples, while marked fluctuations were revealed during sampling period in the counts of enterococci, Enterobacteriaceae and coagulase-positive staphylococci groups (Table 4), with the lowest average value always found in the experimental samples. Significant differences (p ≤ 0.05) among control and experimental samples were revealed for total mesophilic bacteria in April and May, while the yeast count showed similar values in both samples, reaching average value of 5.17 log cfu/g (Table 4).

High-throughput sequencing output determined the relative abundance of bacterial species in Provola cheese samples esteemed based on denoised and taxonomically assigned reads.

A total of 496,430 reads from 10 cheese samples were retained after denoising steps. Sample metadata, denoising statistics, the taxonomy classification, and the relative abundances for each taxonomic level are reported in Supplementary File 1.

The mostly abundant genus was Streptococcus, which accounts for the great majority of relative abundance (99%) and the contribution of other detected taxa was very low. Noteworthy, among those genera that contributed with abundance percentages lower than 1%, Lactobacillus, Lactococcus, and Acinetobacter were found in all samples and can thus be considered as “core” taxa, whereas unclassified genera belonging to the family of Moraxellaceae and Chromohalobacter were not evenly distributed. Specifically, the Chromohalobacter genus was only detected in the experimental sample gathered in June. The great majority of taxa resulted not statistically significant in the performed Welch’s statistical test. The only one statistically significant difference found at the phylum level was between the control and treated groups sampled during the summer seasons. Specifically, although without a statistical significance after multiple test correction, Proteobacteria abundance decreased in experimental samples with respect to controls (Supplementary Figure S1). No other significant taxa were found when other taxonomic levels were inspected.

No statistically significant differences were observed in alpha diversity index when control and experimental groups were compared, i.e., the two groups have comparable distributions both in Shannon and Faith’s PD metrics. Beta diversity distances including Bray-Curtis, Jaccard, and Weighted and Unweighted UniFrac distances did not showed divergent groups. The PCoA plot based on the unweighted UniFrac computed distances reveal no clustering. The pairwise PERMANOVA statistic based on Bray-Curtis distance values returned absence of statistical significance (Supplementary Figure S2).

Looking for a possible stratification of samples based on available metadata, the complete raw panels of detected VOCs and FoodScan physicochemical parameters were inspected by using a Principal component analysis (PCA). As explained by the PCA plot (Supplementary Figure S3), the two plotted principal components that roughly accounted for the 55% of the total distribution, seasonality does not allow for a good cluster separation based on sampling months, whereas a partial division of control versus experimental samples appeared.

The plotted vectors, reported in the PCA biplot graph, evidence a stronger contribution of VOCs in determining the sample placement onto the quadrants. Some fatty acids from physicochemical analysis also contributed (C18:2 n6 trans and C15:1 in the first quadrant).

To obtain a better detailing of VOC and FoodScan variable contribution as separate matrices, sample variables from the two analyses were singly analyzed in an Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) regression analysis. As a matter of fact, the two plots clearly show how the two groups belong to distinct plotted clouds. By separating predictive from non-predictive variation, the resulted VOC T-score and its orthogonal derived values totally accounts for the 40% of the distribution and precisely describe approximately the 29% and 12% on Y and X, respectively (Figure 2). To statistically support this evidence, a permutation analysis of one predicted (p1) and three orthogonal components (o1, o2, o3) that practically relies in the measurement of the observed cross-validated R2X, R2Y, and Q2 coefficients is reported in Supplementary Figure S3.

More in detail panel A of Supplementary Figure S4 reported the most contributing variables, based on the “scores of variable importance on projection” (VIP scores). This weighted sum of OPLS-DA loading squares considers the amount of explained Y-variable in each dimension. 2-Octene, p-cymene, o-cymene, trans-carane, and carane were the most contributing VOCs in terms of VIP score.

In parallel, FoodScan detected cheese variables were investigated by made use of the same multivariate model approach. In this case, the OPLSDA plot reveals some interesting physicochemical variables that much more than other contributed in separating the experimental and control groups (Figure 3). The most contributing FoodScan discriminating variables were evidenced at the top of OPLS-DA VIP score graph (Supplementary Figure S5) and precisely are erucic acid (C22:1 n9), octadecanoic acid (C18:0), polyphenols, tetradecenoic acid (C14:1), octadecatrienoic acid (C18:3 n3), eicosanoic acid (C20:0), and cetoleic acid (C20:1 n9).

To inspect the significance of VOC with a dedicated two group statistics, we made use of a non-parametric corrected (BH) White’s test. The difference in mean proportions was meaningful of different abundance in experimental versus control samples. Acetoin, hexanoic acid, butanoic acid, hexanal, and 2-ethyl-1-hexanol were detected to be higher in experimental group, whereas p-cymene, o-cymene, D-limonene, trans-carane, and beta pinene resulted to be increased in the control group. Figure 4 resumes the mean proportions together with their differences in control and treated groups (corrected p < 0.05). The VOC set from OPLS-DA analysis quite totally overlapped with the variables resulted after the White’s non-parametric test application.

Among FoodScan detected fatty acids, C18:0, C18:3 n3, C20:1 n9, and C22:1 n9 were significantly higher in experimental samples.

The two sets of significant variables (VOC and physicochemical parameters from FoodScan) were merged in a single matrix that was used as input for a subsequent correlation analysis (corrected Person’s correlation test). The analysis returns few interesting positive and negative statistically significant correlations.

More in details, as evidenced in Figure 5, the atherogenic index (relative to the monounsaturated fatty acids) positively correlated with D-Limonene and negatively with acetoin. In turn, acetoin positively correlated with polyunsaturated fatty acids of n3 series (∑ n3) and alpha-linolenic acid (C18:3n 3) and negatively with thrombogenic index, percentage of humidity, and sum ∑ SFA.

Both p-cymene and trans-carane positively correlated with C20:3 n3, whereas the first negatively correlated with polyphenols. Finally, X2 Ethyl 1 hexanol and butanoic acid were negatively correlated with Hexadecenoic acid (C16:1) and ∑ SFA, respectively.