A Study of Greek Graviera Cheese by NMR-Based Metabolomics

1. Introduction

"Graviera" cheese is one of the most popular types of cheese consumed in Greece, along with "feta" cheese, and is produced by dairy farms located throughout the country. It is traditionally made from 80% ewe and 20% goat milk, however graviera cheeses can be prepared also either solely from ewe, goat and cow milk or from a mixture of different milks. Although graviera is produced all over the country, only three regions are certified to produce graviera cheese with the Protected Designation of Origin (PDO) label, as dictated by Greek Law and European Union legislation [1], namely "Graviera Agraphon", "Graviera Kritis" and "Graviera Naxou". Graviera "Agraphon" and "Kritis" are produced using ewe or a mixture of ewe/goat milk, whereas Graviera "Naxou" is produced using cow or cow and ewe milk. According to PDO label rules, PDO Gravieras must contain only milk from animal farms within the geographical region indicated.

During cheese production three important biochemical processes take place, glucolysis, proteolysis and lipolysis, liberating in the cheese lactate, amino acids and glycerides/free fatty acids respectively. These chemical compounds play a crucial role in determining the product’s quality since they contribute to and determine the flavor and aroma characteristics of cheese. Lactate can be metabolized by bacteria present in the cheese curd to products that may affect cheese texture by promoting eye formation (CO₂) or by increasing pH and eventually soften the interior of camembert-type cheese (CO₂, O₂). Lactate metabolism products can also cause defects in cheese texture and flavor, such as "late gas blowing", which refers to cracks formed in the cheese mass during ripening by butyrate and H₂, a procedure accompanied by the development of off-flavors[2]. Proteolysis refers to cheese protein hydrolysis by proteinases, enzymes that can either originate from bacteria indigenous to milk that have survived pasteurization or added during cheese-making, such as starter or non starter lactic acid bacteria, rennet, etc. The process of proteolysis contributes to cheese flavor and texture by breaking down the protein network and releasing amino acids [3]. Lipolysis is the hydrolysis of triglycerides by hydrolases producing free fatty acids and minor glycerides such as diglycerides or monoglycerides. Hydrolases are classified as lipases or esterases depending on the nature of substrate, the length of fatty acid chain and enzymatic kinetics. Free fatty acids, especially short and intermediate-chain FA, are responsible for the characteristic flavor and aroma of cheese, depending on their perception threshold, concentration and pH [4].

There are a number of factors that can have a significant effect on cheese composition and quality. Herd diet, which is usually based on local flora, can affect milk composition and microbiological content, therefore it is considered as one of the main factors that contribute to cheese amino acid and fatty acid content. For example, Tzora et al. demonstrated that dietary treatment of sheep enriched with omega-3 fatty acids can affect the omega-3 fatty acid content of cheese as well as the bacterial populations [5]. Alpine Asiago PDO cheeses produced with milk originating from pasture fed cows were characterized by higher amounts of lysine, choline and 2,3-butanediol, indicating the effect of animal feed on cheese composition [6].

The development of cheese quality characteristics can be also influenced by the pedoclimatic conditions in the geographical region the cheese is produced. A number of studies have been devoted to the analysis of greek graviera cheese obtained from different areas of production. Using ICP-MS, Danezis et al. determined the elemental profile of greek graviera cheeses and classified them into nine geographical categories/regions [7]. Notably, cretan gravieras were found to contain higher levels of praseodymium (Pr) and neodymium (Nd), a finding attributed to soil and vegetation in the island of Crete being rich in these elements. Another study of the microbiological and physicochemical analysis of greek graviera cheese showed that "Graviera Kritis" was characterized by lower pH values and amino acid content than "Graviera Naxou", possibly due to different NSLAB acidifying and proteolytic activity [8]. Analysis by means of SPME-GC-MS in the same study showed that although the two graviera cheeses had similar VOCs profiles, there were some compounds uniquely identified in each cheese label. Vatavali et al. analysed the physicochemical properties, mineral content, fatty acid composition and VOC profile of graviera samples produced in six different regions of Greece [9] and then later expanded the work to include another five Graviera-producing regions.[10] The statistical analysis of the combined analytical data set showed that gravieras from Naxos was the most clearly differentiated group of samples, presumably due to milk composition and geographical differences [10].

Furthermore, different cheese types require different milk treatment, specific production processes and precise maturation conditions for the final product to occur. The lipid fraction of a large number of different greek PDO cheeses were studied using GC-MS, and the collected data was used to study whether physicochemical properties and fatty acid profiles could act as markers of PDO label, milk and cheese type discrimination [11]. In a study of the free fatty acid profile of traditional greek cheese varieties, Georgala et al. reported that Cretan graviera cheese contained more propionic acid while Kefalotyri had a higher acetic acid content [12].

¹H NMR spectroscopy has proven to be an extremely useful tool in cheese analysis, contributing in the determination of the lipid fraction as well as the water-soluble metabolite content of different types of cheese [13,14,15]. The metabolite profile obtained by NMR spectroscopy was used successfully to study the ripening stage of Grana Padano [16] and Fiore Sardo [17] cheeses from Italy. Samples of Parmigiano Reggiano cheese were successfully differentiated from other "grana type" cheeses produced in Eastern Europe by means of NMR spectroscopy combined with multivariate analysis, despite the fact that the cheeses were at a different ripening stage [18]. Using data obtained from a variety of analytical techniques, including ¹H NMR spectroscopy, Brescia et al. were able to discriminate PDO from PGI samples of mozzarella cheese utilizing untargeted multivariate statistical analysis [19]. Likewise, the production site of Asiago d' Allevo cheese samples was identified by analyzing the ¹H NMR spectra of their organic extracts using untargeted multivariate analysis [20]. It has also been demonstrated that European Emmental cheeses can be discriminated according to their geographical origin by HR-MAS NMR spectroscopy [21]. Samples of "Mozzarella di bufala Campana" produced in different sites, yet included in the PDO geographical region of the cheese, were discriminated by the same technique [22].

To our knowledge, the polar and lipid metabolite profile of greek graviera cheese has not been studied using the analytical NMR spectroscopy methodological approaches described above for a variety of other cheeses. The aim of this work is thus firstly to characterize the full polar and apolar metabolite profile of greek graviera cheese, produced in the mainland, Crete, and the Aegean islands, and to examine the ability of liquid phase ¹H NMR spectroscopy in combination with statistical analysis to classify graviera cheese according to area of production, with emphasis on the authentication of Cretan graviera.

2. Results and Discussion

2.1. Analysis of ¹H NMR spectra

2.1.1. Polar fraction

A representative ¹H NMR spectrum of the aqueous extract of greek graviera cheese is depicted in Figure 1. Thirty two compounds were identified in the ¹H NMR spectra of the aqueous extracts and their detailed peak assignments are listed in Table 1. Through the use of a suitable internal standard, these compounds were quantified in all the studied cheese samples by integrating their NMR spectra. Table S1 provides detailed production details and meta data for all samples studied, while Table S2 summarizes the mean values of these polar metabolites as a function of graviera geographic origin.

In general, the composition of greek gravieras was found to be similar to those of other hard and semi-hard yellow cheeses [13,16,19,21,23]. The graviera polar metabolite profile contained the organic acids lactate (8-22 g/Kg), citrate (0-1572 mg/Kg), succinate (0-1340 mg/Kg) and acetate (0-1238 mg/Kg), in order of decreasing average value. The acetate content of the gravieras determined in this study was in the range reported by Georgala et al. [24] for Graviera Kritis (727-1074 mg/Kg). In a previous study, Dudley and Steele[25] reported that some nonstarter Lactobacillus plantarum strains are able to produce succinate by citrate catabolism in Cheddar cheese. Succinate was present in all but one graviera cheese product in this study, and at levels similar to those reported for Emmental cheese [26], while citrate was present in 43 out of 74 graviera cheese samples. Based on literature, the presence of citric acid in cheese is highly affected by the presence of certain lactic acid bacteria species and their ability to ferment it to acetate [27]. Benzoic acid was the only aromatic acid identified in the graviera samples, ranging between 3-34 mg/kg of cheese. In a recent study Yerlikaya et al. examine the formation of benzoic acid and its relationship to the microbial properties of traditional Turkish cheeses [28].

The graviera samples were found to contain 18 out of 20 proteinogenic amino acids, products of proteolysis, namely Val, Ile, Leu, Ala, Lys, Pro, Met, Glu, Gln, Phe, Thr, Tyr, Gly, Trp, Asn, Asp, Arg, Ser. Almost the same amino acids were reported for Swiss cheeses [29], but with His, instead of Trp found in the greek graviera. Sixteen of the eighteen amino acids identified in greek graviera were reported in Grana Padano cheese using ¹H NMR spectroscopy, missing Trp and Arg.[16]. Using GS/TOF-MS, approximately the same amino acids have been identified in various types of hard cheese, but with His instead of Arg [30]. The amount of total free amino acids (TFAA) in the graviera samples varied significantly, ranging between 0.8 and 54 g/kg, with a mean value of 13 g/kg.

Some non proteinogenic amino acids, namely ornithine and γ-aminobutyric acid (GABA) were also identified in the polar extracts of greek gravieras. GABA is an L-glutamate metabolite that has been quantified at high concentrations (up to 7g/Kg) in Cheddar cheese [31], and at lower concentrations (<400 mg/Kg) in 22 different Italian cheeses. In the present study it was identified in 27 of 74 gravieras at levels up to 3000 mg/Kg. Ornithine, a product of arginase activity on arginine, was identified in 65 of 74 samples, with the highest ornithine content determined in a cretan graviera sample (2.2 g/kg). These two amino acids, γ-aminobutyric acid and ornithine have recently been the focus of studies regarding their bioactive functions and impact on health [32].

Tyramine, a decarboxylation product of tyrosine, was also identified in a large number of graviera samples, evidently due to the presence of decarboxylase-positive microflora in the cheese. It has been reported that the tyramine content of Dutch type cheese was highly affected by storage time and storage temperature [33], so this biogenic amine may represent a helpful marker of the aging stage of cheeses. Cretan graviera was not found to contain any sugars, in agreement with data reported in literature for similar hard type yellow cheeses.

2.1.2. Lipid fraction

The ¹H NMR spectrum of the lipids extract of greek graviera is depicted in Figure 2, with the spectral region between 5-6.5 ppm in expanded fashion, while the peak assignment of the apolar graviera components is listed in Table 2. The ¹H NMR spectrum of graviera is very similar to those of other cheeses [14,15], and contains all the expected proton peaks of long and medium chain fatty acids, as well as peaks attributed to butyric acid, caproleic acid and the conjugated linoleic acid, rumenic acid (CLA), which are fatty acids typically identified in animal-derived lipids. CLAs have been reported to be significantly increased in milk fat from cows grazing pasture compared to typical dairy diets, thus CLAs in cheese could be related to animal feeding type [34]. It is worth noting that although caproleic acid and CLAs have been identified by NMR in greek bovine milk [35], only CLAs were reported in a recent GC-MS study of the FA profile of Greek protected designation of origin (PDO) cheeses [11], while neither caproleic acid nor CLAs were reported in recent studies of the FA profile of greek gravieras [10,36]. The CLA content of greek cheeses and its evolution during aging has also been the focus of earlier studies [37,38].

Finally, NMR spectroscopy also provides access to the glyceride profile of cheese, analytical data which are not available from GC-MS analyses of fatty acids. The characteristic peaks of 1,2-diglycerides were identified in the graviera cheese apolar extract spectra, indicating the activity of lipase enzymes that hydrolyze triglycerides producing free fatty acids, diglycerides and monoglycerides (MG). Additionally, both 2-MG and 1-MG monoglyceride peaks were identified at very low levels in several (18 of 74) cheese spectra, possibly as a result of factors irrelevant to cheese making [39], since no repetitive pattern concerning milk type, area of production or maturation was observed.

The ¹H NMR spectra of the lipid extracts were integrated and the data were used to calculate the fatty acid profile of the graviera samples as % moL of total fatty acids according to a previously published methodology [40], with the relevant equations modified to include fatty acids specific to cheese lipids. Rumenic acid (CLA) and caproleic acid provide independent proton signals (G and H respectively, letter codes in Table 2) that can be used to quantify them directly, while a combination of signals is needed for the quantification of the rest of the fatty acids. For example, the percentage of ω-3 polyunsaturated acids, mainly linolenic acid (LN), in the samples was calculated using the following equation:

$$\left[\mathrm{L}\mathrm{N}\right]=\frac{\mathrm{E}}{\mathrm{E}+\mathrm{F}+\mathrm{I}}\times 100$$

(1)

Accordingly, butyric acid (BA), linoleic acid (LO), monounsaturated acids (MUFA) and saturated fatty acids (SFA) were calculated by the following adjusted equations:

$$\left[\mathrm{B}\mathrm{A}\right]=\frac{\mathrm{I}}{\mathrm{E}+\mathrm{F}+\mathrm{I}}\times 100$$

(2)

$$\left[\mathrm{L}\mathrm{O}\right]=\frac{\left(\mathrm{A}-\left(4\times \frac{\mathrm{E}}{3}\right)\right)}{B}\times 100$$

(3)

$$\left[\mathrm{M}\mathrm{U}\mathrm{F}\mathrm{A}\right]=\frac{\mathrm{0,25}\times \left(\mathrm{C}-\left(4\times \mathrm{G}\right)-\left(2\times \mathrm{H}\right)-\left(4\times \frac{\mathrm{E}}{3}\right)-\left(4\times \frac{\mathrm{A}-\left(4\times \frac{\mathrm{E}}{3}\right)}{2}\right)\right)}{\frac{\mathrm{B}}{2}}\times 100$$

(4)

$$\left[\mathrm{S}\mathrm{F}\mathrm{A}\right]=[\frac{F}{E+F+I}-\frac{\left[CLA\right]}{100}-\frac{\left[LO\right]}{100}-\frac{\left[MUFA\right]}{100}]\times 100$$

(5)

The quantitative analytical data obtained for the lipid fraction of the graviera samples studied are also reported in Table S2.

2.2. Analysis of Variance (Anova) of graviera composition

In order to examine compositional differences between gravieras produced in different geographical locations, the polar and apolar compositional data obtained by NMR spectroscopy were analysed by Anova. Focusing on gravieras produced in the island of Crete, Table S3 lists the relevant values of the factor F and the level of significance p for each metabolite. Several metabolites show statistically significant concentration differences at p level <0,05, indicating a large variability in chemical composition is present already within the island of Crete, and some representative cases are depicted in Figure 3. The distribution of γ-aminobutyric acid according to geographic origin, for example, appears to be very interesting, as GABA is present in 11 of 19 gravieras from Rethymnon, while only 5/36 gravieras from the rest of the island were found to contain any. Recently, Tsafrakidou et al. [41] studied the phenotypic and genotypic variability of lactobacilli isolated from mature Graviera Kritis produced at two traditional dairies, concluding that lactobacilli are associated with their particular production dairy ecosystem. Since the production of increased GABA has been attributed to presence of selected NSLAB [31], this difference in Rethymnon gravieras may be related to locally present LAB strains.

It is worth noting that gravieras from eastern Crete (Heraklion, Lasithi) contain almost double the amount of free amino acids compared to western Crete, and relatively lower amounts of saturated fatty acids. Focusing on fatty acids that have not been reported in the literature so far in detail, we observe that gravieras produced in Crete contain similar amounts of CLA regardless of origin within the island, while gravieras produced in Heraklion have the lowest amounts of caproleic acid and the highest amounts of linoleic acid.

Further analysis by Anova showed that metabolite profile differences were also present for samples produced in different Aegean islands (Naxos, Tinos, Lesvos, Ios) but the small number of samples (n=2) available for each island only allows for qualitative analysis at present. Nevertheless, it should be pointed that the CLA content of gravieras produced in Naxos (0.51±0.1 %) and Tinos (0.54±0.09 %) from cow milk was found to be significantly lower than those produced from sheep/goat milk in Crete (0.92±0.2 %). The increased levels of CLA when sheep/goat milk is used for cheese-making are well documented in the literature and for a variety of cheeses [42]. To allow for a better comparison between wider geographical regions, gravieras produced in mainland Greece were grouped together, and compared with the Cretan and Aegean groups as a whole. The Anova analysis results, presented in Table S4, demonstrate that statistically significant differences between mainland, cretan and aegean gravieras are present in the NMR metabolite profile, with the most important differences observed for CLA, GABA, 1,2-DG, sterols, and LN, as depicted in Figure 4. Cretan gravieras contain larger amounts of CLA and ω3-polyunsaturated fatty acids, but lower amounts of caproleic acid and 1,2-diglycerides, compared to aegean and mainland gravieras.

An attempt was also made to examine possible correlations between the metabolite profile of the gravieras with their maturation time, which varied between 3 months (the minimum allowed by law) and 22 months, with most samples aged less than 12 months. Aged gravieras possess superior organoleptic properties and are thus considered of increased commercial value. A Partial Least Squares (PLS) model of maturation time was developed using the metabolite profile obtained through NMR, and Figure 5 depicts a plot of predicted vs actual maturation time based on this PLS model, along with the Variable Importance Parameter (VIP) values of the metabolites. The metabolites mostly correlating with increased maturation time where (VIP>1) were glycerol and several amino acids, indicating the importance of lipolysis and proteolysis in developing the special sensory attributes of aged graviera cheese. Although the correlation of the graviera NMR metabolome with maturation time through the PLS model of Figure 5 is not entirely satisfactory, it is worth noting the PLS model appears to cluster samples regardless of geographical origin.

2.3. Discriminant Analysis

As demonstrated by Anova in the previous section, there is a large within-region metabolite variability for gravieras produced within Crete and the Cyclades, thus in order to facilitate the multivariate analysis of the data set as a whole, samples were initially categorized in three broad geographical areas (Crete, Cyclades, mainland). First, a targeted approach was pursued, using OPLS-DA unsupervised multivariate analysis of the polar and apolar metabolite data of the graviera samples. The score plot of this OPLS-DA analysis is depicted in Figure 6, and it can be observed that the metabolite profiles of the gravieras are able to classify them according to broad geographical origin successfully. Gravieras from mainland Greece are more widely distributed than the two island origins, as expected due to the significantly larger geographical area covered under this group. Gravieras originating from Crete, which represent the larger sample group, show good clustering. The loadings plot of the OPLS-DA model in Figure 6, in which the Variable Importance parameter (VIP) is reflected in the size of the metabolite data points, shows that the most important metabolites for the observed classification were 1,2-DG and sterols for the Cyclades group, and GABA for the mainland group, while cretan gravieras were differentiated mainly due their fatty acid composition differences (CLAs and butyric acid).

Since the majority of the graviera samples available in this study were produced in the island of Crete, it was interesting to explore whether an OPLS-DA model could be obtained from the NMR data that could classify the graviera samples based on a dual cretan/non-cretan origin. Such models would be very useful as a first step towards future efforts to construct chemical composition databases for the authentication of Cretan PDO gravieras. To pursue this goal two different OPLS-DA models were constructed, one based on metabolite data (targeted), and one based on spectral bucketing of the whole NMR spectra (untargeted, polar plus apolar extract spectra), and Figure 7 depicts the respective score plots of these two exploratory OPLS-DA models. Both modeling approaches correctly classify the gravieras based on Cretan/non-cretan origin hold promise in establishing authentication models for Cretan graviera based on NMR metabolomics. The untargeted approach based on spectral bucketing has a higher discrimination ability, as reflected in the CV-ANOVA analysis of this model, depicted in Table S5, with a CV-ANOVA p value of were 2.9x10^-6 and the respective permutation test, depicted in Figure S1. It understood that non-Cretan gravieras sampling is under-represented in these models, thus current work in our group is focusing on enlarging the mainand gravieras sampling space. The large variability in graviera chemical composition observed in this study emphasizes the importance of multiplatform approaches in geographic origin discrimination studies of foods in order to take advantage of not just the metabolome, but also the proteome and metalome of the product in study.

3. Materials and Methods

3.1. Chemicals

The standard reagents deuterium oxide, D₂O (99.9 atom % D, contains 0.05 wt % 3-trimethyl-silyl propionic-2.2.3.3-d₄ acid sodium salt TMSP) and chloroform, CDCl₃ (99.8% atom % D contains 0.03 % (v/v) tetramethylsilane TMS) were obtained from Deutero GmbH, Germany.

3.2. Samples

For this study, 74 hard cheese samples were used for NMR analysis, and were obtained either directly from local dairy farms or bought from local markets. 53 cheese samples were produced in the island of Crete, while the remaining 21 were produced in northern Greece in the islands of Mitilini, Naxos, Tinos and Ios, as well as in Peloponnesus, Macedonia, Epirus, Thrace, Thessaly and Central Greece during the cheesemaking periods of 2012-14 and 2017-19. The maturation time of the cheese samples varied from 3 months (minimum period of maturation) up to 22 months.

3.3. Sample Preparation-NMR analysis

The graviera samples were prepared according to a published solid food NMR analysis protocol [43], with small alterations. 1-D and 2-D ¹H and ¹³C NMR spectra were obtained on a Bruker Avance-III-500 spectrometer operating at 500.137 MHz for proton and 125.75 MHz for carbon. The water signal was suppressed by pre-saturation , when necessary. The following conditions were used for the acquisition of 1-D ¹H spectra: 64k data points, 256 scans, 8 dummy scans, spectral width of 12.016 ppm, line broadening of 0.3 Hz. For ¹³C dept135 spectra, the conditions were: 65000 data points, 7168 scans, 4 dummy scans, spectral width of 160 ppm, line broadening of 1.0 Hz. Phase and baseline correction were applied to all 1-D spectra and 2-D spectra where necessary. All spectra were processed by standard TopSpin Software (Bruker, v3.1). Signal assignment was achieved by standard 2D NMR gradient spectroscopy (gCOSY, gHSQC, gHMBC, gTOCSY, gHSQC-TOCSY) experiments, comparison with online databases and the Chenomx NMR Suite software (8.02).

3.4. Spectra integration and multivariate statistical analysis

For targeted analysis, the ¹H NMR spectra of aqueous and chloroform extracts were integrated manually using TopSpin software and the option for automated baseline correction of the integrals was selected for better accuracy. 3-(Trimethylsilylpropionic-2,2,3,3-d₄ acid sodium salt was used as internal standard for the quantification of the polar graviera extracts. Representative LOQ for signals with S/N=10 were at 0.05 μmoL/g of cheese, and the relative standard deviation was calculated via replicant analysis to be <15% for the metabolites quantified.

For untargeted analysis, the ¹H NMR spectra of aqueous and chloroform extracts were bucketed using Amix Software (Bruker, v3.9.14) and the data were imported to Simca Software (Umetrics, v13.02) to develop multivariate statistical analysis models (PCA, OPLS-DA). Internal cross validation of the OPLS-DA models was performed automatically by the software by repeatedly dividing the samples in 7 random groups, developing a DA model using the 6 groups (training set) and confirming its validity with the last group (test set) until a stable model was produced.