Fatty Acid and Stable Carbon Isotope Composition of Slovenian Milk: Year, Season, and Regional Variability

This study examined the percentage and stable isotope ratios of fatty acids in milk to study seasonal, year, and regional variability. A total of 231 raw cow milk samples were analyzed. Samples were taken twice per year in 2012, 2013, and 2014, in winter and summer, covering four distinct geographical regions in Slovenia: Mediterranean, Alpine, Dinaric, and Pannonian. A discriminant analysis model based on fatty acid composition was effective in discriminating milk according to the year/season of production (86.9%), while geographical origin discrimination was less successful (64.1%). The stable isotope composition of fatty acids also proved to be a better biomarker of metabolic transformation processes in ruminants than discriminating against the origin of milk. Further, it was observed that milk from Alpine and Mediterranean regions was healthier due to its higher percentage of ω-3 polyunsaturated fatty acid and conjugated linoleic acid.


Introduction
There is an increasing interest in the quality of foodstuffs, and agricultural research continues to focus on improving the nutrient profile of food products. Milk is a valuable commodity that is consumed directly, as well as being a minor and major ingredient in many food products [1]. It is also rich in nutrients, a resource of proteins (2.6-4.2%), vitamins (0.1%), minerals (0.7%), and fatty acids (95-96% of total milk fat), and is considered by many as an essential part of a healthy diet [2,3]. In Europe, milk provides 9% of the dietary energy supply, 19% of the dietary protein supply, and between 12-14% of dietary fat supply [4].
According to the Agriculture and Rural Development, Analysis and Outlook Unit European's farms produced 172.2 million tonnes of raw milk in 2018, a number that they expect to increase to 179 million tonnes by 2030, but at a slower pace (+0.6%) than in 2008-2019, and represents a year-on-year increase of 1.6 million tonnes (+1% per year on average) [5]. In Slovenia, dairy farmers produced 564,000 tons of cows' milk in 2019 (1.2% less than in 2018), and 164,000 tons of drinking milk, which represents a decrease of more than 2% compared to 2018. Despite the decrease, imports increased by almost 13% [6]. Consumers prefer milk with favorable sensory qualities, such as flavor, which depends on animal management factors, including breed, lactation stage, processing techniques, and diet regime. However, diet is considered the primary and the most sensitive factor responsible for creating a specific milk flavor [7] and also a significant factor affecting its fatty acid (FA) composition. These differences also reflect differences in regional dairy production systems, and linking a milk product to a specific region is of high importance as a recognition of quality [1]. Quality linked to geographical origin has created significant interest in building local and regional food systems across Europe, including Slovenia.
Lipid content and FA profile are considered to have an important role in the human diet [7,8]. For example, certain milk FAs, such as linoleic and α-linolenic FAs (C18:2ω-6 and C18:3ω-3), are essential and cannot be synthesized by the human body [2]. Further, a lower proportion of saturated fatty acids (SFA) and a higher proportion of unsaturated fatty acids (UFA), especially polyunsaturated fatty acids (PUFA), is desirable [9]. The consumption of SFA (especially lauric (C12:0), myristic (C14:0), and palmitic acids (C16:0)) [9] is associated with cardiovascular disease, while UFAs are considered as beneficial for human health. Other functional food components are omega-3 polyunsaturated FA (ω-3 PUFAs), which improve neurological functions and are essential in early human development years, and conjugated linoleic acid (CLA), which has anti-carcinogenic, anti-diabetic, anti-atherogenic effects [10,11].
Farmers can significantly alter the FA composition of milk by managing the nutrition of their dairy cows [12]. Fresh forages (especially grass) in the diet enhance the concentration of unsaturated FA, such as α-linoleic acid and CLA, and decrease the levels of SFA, while adding grass conservation products (hay or silage) increases the proportion of SFA in milk [13]. Besides, the total SFA concentration of milk fat is usually lower in summer compared to winter due to the higher consumption of summer forage crops that have a higher content of PUFAs. The higher content of monounsaturated fatty acids (MUFA) and PUFA in milk from cows grazing in highland regions can be attributed to differences in the botanical composition since higher botanical diversity and energy shortage in cows grazing in mountain pastures modifies the bacterial population and lipid accumulation of UFAs in the rumen [14]. Fatty acid composition of milk fat has also been used as a criterion to detect adulteration of milk, such as the substitution of part of the fat or proteins of plant or animal origin and admixtures of milk of different species. The method to characterize the matrixes according to the animal species can be greatly improved when combined with the content of triacylglycerols [15].
The effects of feeding regime can also be tracked using stable isotope analysis using isotope ratio mass spectrometry (IRMS). In dairy products, the isotope composition is also primarily determined by the cows' diet, as well as their metabolisms [16]. Differences in stable C isotope abundance is mainly the result of the different ingested proportion of feeds derived from C 3 and C 4 plants. The metabolism of C 4 plants, such as corn or sugar, uses a different biosynthetic pathway to fix atmospheric CO 2 compared to C 3 plants (grass and cereals), resulting in the accumulation of the heavier 13 C isotope compared to C 3 plants, which leads to higher δ 13 C values. These differences also show up in the isotopic composition of FAs [17,18]. Since the spectrum of plants used for feeding purposes depends mainly on local soil and climatic condition, stable isotope composition can also encompass information about the geographical origin of dairy products useful for authenticity control.
Up to now, studies have focused on stable isotope ratio measurements of δ 13 C, δ 15 N, and δ 34 S values in casein, δ 2 H, δ 18 O in milk water [19][20][21][22] and Sr in dry milk samples [21]. One of the most recent significant advances has been the development of compound-specific stable isotope analysis (CSIA) that enables the identification of relevant biomarkers, such as fatty acids, in food authenticity studies. Some studies have even looked at the correlation between long-chain FAs (α-linolenic FA-C18:3ω-3, eicosapentaenoic FA-C20:5ω-3) and stable isotope measurements to distinguish between organic and conventionally produced milk [17,23]. In another study, the authors investigated the relationship between the carbon isotope compositions of CLA in milk and fatty acids in different dietary sources (C 3 versus C 4 plants) to understand the transformation of FA from diet to milk [17]. The results suggested that the δ 13 C-change in CLA did not originate from the microbial biohydrogenation in the rumen, but most probably from endogenous desaturation of t11-18:1.
Only a few studies, however, have focused on δ 13 C and δ 2 H values in milk fatty acids. Ehtesham and co-workers [24] demonstrated that δ 2 H values in precipitation were highly correlated with the δ 2 H values in bulk milk and specific fatty acids (butyric, myristic, palmitic, and oleic acid). Separation of milk powder origin to geographic sub-regions within New Zealand has been achieved, with clear discrimination between South Island and North Island and with significant clustering of regions. Further, a linear discriminant analysis (LDA) model based on δ 2 H and δ 13 C values of butyric, myristic, palmitic, and oleic fatty acids has provided the best separation of samples from the Northern and Southern parts of New Zealand [1]. Ehtesham et al. [25] also focused on the relationship between δ 2 H and the FAs in milk and the water component of milk (feed and farm water). The authors found that the δ 2 H values of FAs in feed and δ 2 H values of farm water were related to the δ 2 H values of milk FAs, milk solids, and milk water and could be used as a biogeochemical marker surrogate to environmental conditions. This approach can also be applied to determine the dietary regime and geographical origin of milk products where the water has been removed (milk powder) or is altered in the case of processed milk products (cheese and butter).
The main objective of this study was to evaluate both the need and options for origin tracing, focusing on Slovenian milk. In 2016, the Slovenian Ministry of Agriculture, Forestry, and Food began a three-year project to promote and raise awareness of the importance and characteristics of locally produced and processed foods. During this time, the protective mark "Selected Quality" national scheme and the "Selected Quality-Slovenia" was established. Nowadays, most milk and dairy products, produced and processed in Slovenia, use the "Selected Quality -Slovenia" protective mark. This mark indicates that the products have the following distinguishing properties: (i) the origin of raw milk-the milk and other dairy products are produced entirely in the same country, Slovenia; (ii) milk quality-excellent microbiological quality of milk; (iii) freshness of raw milk-it's 15 h at most from the takeover of milk at the collection site to its acceptance into the processing plant. Therefore, to protect the obtain mark and prevent possible fraud, a robust screening method to determine the authenticity and regional traceability of milk is needed. In this study, two methods (fatty acid composition and isotope ratios of carbon within fatty acids) were applied. Additionally, the effect of the production year and season on the FA composition of milk was also investigated. Finally, the FA profiles pertaining to the milk's essential and beneficial components were also studied, and a comparative study was made to examine the differences in the quality of milk from different geographic regions in Slovenia.
A lower C16:0 and HFA percentage, as well as a higher amount of CLA and C18:1ω-9, most likely results from increased grazing [9]. However, the amount of each FA depends on lipid metabolism in the rumen and the mammary gland and the release of FAs from body reserves during a period of negative energy balance during early lactation since diet (especially type of forage, lipid supplements, forage:concentrate ratio, starch level, and their interactions) can influence metabolic processes in the rumen and the composition of rumen biome [9]. Variations in FA composition can also be explained, in part, by the micro-environmental conditions influencing the quality of feed. The composition of forage and plant lipids is influenced by factors, such as species, growth stage, temperature, and light intensity. Besides, environmental temperature and humidity can affect both milk yield and FA composition since dairy cows are sensitive to climatic variations. The production of the short-chain FAs (C6:0-C8:0) and C18:0 decrease at the higher temperatures and humidity. In this study, 2013 and 2012 were exceptionally hot and wet years, respectively, and there were no differences in the levels C6:0-C10:0 FAs ( Table 1).
The results showed no significant seasonal effects (p > 0.05) on the percentage of C15:0, C16:1, C20:1ω-9, C20:3ω-6, C21:0, C20:3ω-3, C20:4ω-6, C20:5ω-5, C22:0, and C24:0, although the percentage of total SFAs was lower in summer (68.04 ± 3.20% of total FA) than in winter (70.40 ± 2.06% of total FA, p < 0.05). The percentage of MUFAs of total FA was significantly higher (p < 0.05) in summer (26.30 ± 3.23%) than in winter (24.20 ± 1.87%) due to increased amounts of C18:1 in summer compared to winter (p < 0.05). The percentage of PUFAs in summer (4.11 ± 0.87% of total FA) was higher than in winter (3.43 ± 0.33% of total FA, p < 0.05). On average, amounts of SFA decreased in the summer and increased in the winter, whereas the amount of PUFAs increased in summer and decreased in winter. The percentage of ω-3 PUFA was also higher in summer than in winter. Differences in FA percentage related to the season could be explained by differences in the cow' dietary regimes since during the summer, the cows typically graze on pasture or are fed fresh forage. Replacing grass silages with corn silage increases the concentration of medium-chain saturated FAs (C6:0 to C14:0) and C18:2ω-6 and decreases the concentration of C16:0, C18:0, C18:3ω-3, and C20:5ω-3 in milk fat. Replacing grass silage with corn silage will decrease the amount of total ω-3 PUFA and increase the amount of ω-6 PUFA in milk, leading to an elevated ω-6:ω-3 ratio [12,21,23]. The effect of replacing grass with corn silage on the CLA content of milk is inconclusive. Studies have shown that corn silage-based diets increase the CLA content of milk fat compared with grass-based diets [24,25], while others have found no differences in CLA content of milk fat from cows offered grass or corn silage-based diets [21,26,27]. The present study indicated that the percentage of C18:2ω-6 and CLA was higher during the summer. However, higher environmental temperatures during the summer also affected the FA composition of milk, and milk fat in summer was lower in C16:0 relative to C18:0 than in milk fat from the same cows during winter. Although these changes were related to changes in blood plasma lipids, such observations were confounded by dietary changes.
As part of this study, the data were analyzed using a clustering approach based on the percentage of FA, production year, season, and corn silage and cereals in the diet. The δ 13 C values of casein fraction were also used to characterize the presence of corn silage at a threshold limit of −23.5% , i.e., samples with δ 13 C values higher than −23.5% indicated the presence of corn in the diet [16]. Before clustering, the Silhouette curve method was applied to estimate the optimum number of clusters [29]. In this study, the best split was achieved using 14 clusters. The partition around medoids (PAM) method was then applied, which divided the milk samples into 14 clusters (Figure 1). The results showed that the production year and season and the presence of corn silage and cereals in the animal's diet played a significant role in the clustering. Figure 1 also shows the combined effects of year and season. Production year had a more significant effect on FA composition than the season. Specific FAs showed no global pattern of differentiation and were probably influenced by micro-factors, such as temperature, humidity, feeding regime, lactation stage, breed, and processing techniques. year and season. Production year had a more significant effect on FA composition than the season. Specific FAs showed no global pattern of differentiation and were probably influenced by microfactors, such as temperature, humidity, feeding regime, lactation stage, breed, and processing techniques. Multivariate analysis was used to identify the key factors responsible for discriminating between production year and season. In this case, 231 samples consisting of authentic Slovenian milk samples collected in summer and winter 2012 (12 s, 12 w), 2013 (13 s, 13 w), and 2014 (14 s, 14 w) were statistically evaluated. Figure 2 presents the results of the discriminant analysis (DA) as a discriminant function score plot (a) and a discriminant loadings plot (b). In the functional score plot, each group (centroid) was represented by a scatter plot, while in the loadings plot, they appeared as a set of vectors, indicating the degree of association of the corresponding initial variables with the first two discriminant functions (F1 and F2).
In latter, the degree of distribution of each parameter in the classes was revealed. Except for the 2013 and 2014 winter samples (13 w and 14 w, respectively), good separation between year and season of production was achieved. The most influential parameter for the summer/winter 2012 (12 s/12 w) classes (see top left of biplot a) was C20:4-6 (direction of its vector was towards the right part of the biplot), while CLA separated 12 s from 12 w. Winter samples 2013 (13 w) and 2014 (14 w) overlapped, while the separation between summer classes 14 s and 13 s was likely a result of the higher percentage of C17:0iso, C20:1-9, C18:3-6, C20:3-3 (bottom left (a)). A leave-one-out crossvalidation (LOOCV) classified 86.9% of the samples correctly. Multivariate analysis was used to identify the key factors responsible for discriminating between production year and season. In this case, 231 samples consisting of authentic Slovenian milk samples collected in summer and winter 2012 (12 s, 12 w), 2013 (13 s, 13 w), and 2014 (14 s, 14 w) were statistically evaluated. Figure 2 presents the results of the discriminant analysis (DA) as a discriminant function score plot (a) and a discriminant loadings plot (b). In the functional score plot, each group (centroid) was represented by a scatter plot, while in the loadings plot, they appeared as a set of vectors, indicating the degree of association of the corresponding initial variables with the first two discriminant functions (F1 and F2).
In latter, the degree of distribution of each parameter in the classes was revealed. Except for the 2013 and 2014 winter samples (13 w and 14 w, respectively), good separation between year and season of production was achieved. The most influential parameter for the summer/winter 2012 (12 s/12 w) classes (see top left of biplot a) was C20:4ω-6 (direction of its vector was towards the right part of the biplot), while CLA separated 12 s from 12 w. Winter samples 2013 (13 w) and 2014 (14 w) overlapped, while the separation between summer classes 14 s and 13 s was likely a result of the higher percentage of C17:0iso, C20:1ω-9, C18:3ω-6, C20:3ω-3 (bottom left (a)). A leave-one-out cross-validation (LOOCV) classified 86.9% of the samples correctly.

Influence of the Geographical Origin on FA Composition
Fatty acid composition combined with various chemometric tools has been successfully applied in different food commodities to discriminate the geographical origin, including fruits [30], oils [31,32], and also milk and dairy products [26,33,34]. In this study, discriminant analysis (DA) was used to discriminate between milk from Alpine, Mediterranean, Pannonia, and Dinaric regions of Slovenia using the whole dataset, including summer and winter samples from 2012, 2013, and 2014. The results are presented in Figure 3a

Influence of the Geographical Origin on FA Composition
Fatty acid composition combined with various chemometric tools has been successfully applied in different food commodities to discriminate the geographical origin, including fruits [30], oils [31,32], and also milk and dairy products [26,33,34]. In this study, discriminant analysis (DA) was used to discriminate between milk from Alpine, Mediterranean, Pannonia, and Dinaric regions of Slovenia using the whole dataset, including summer and winter samples from 2012, 2013, and 2014. The results are presented in Figure 3a,b. A comparison of the scatter of data points and vector plots indicated that the Mediterranean group (Figure 2a) was a result of high amounts of C15:0, C17:0, C17:0iso, and C17:0aiso FA. For the Alpine region and Dinaric region groupings (Figure 3c: lower left) C18:3ω-3, CLA, C18:3ω-6, and C20:5ω-3 had the most separating power, while the C12:0, C11:0, C21:0, C18:0, and C18:2ω-6 differentiated the Pannonian group. A cross-validation test showed 64.1% correct classification. The classification rate was the highest (80.8%) for Pannonian milk samples and the lowest for Dinaric samples (20.5%). The more precise and efficient separation between the four regions was obtained when taking into account the season (summer-winter). In this case, the classification rate was >90%. Figure 3c shows the DA plot for 39 milk samples obtained in winter 2013. The analysis generated four well-separated groups. Separation of the regions was also achieved based on FA composition when the data were treated altogether. For F1 (89.0%), C15:0 and its isomers were the main discriminatory parameters, while for F2, the most discriminatory parameters (10.5%) were C18:3ω-6 and CLA. However, cross-validation was low, and only 35.9% of the samples were classified correctly. The highest rate of classification was for the Mediterranean samples (66.7%) and the lowest for the Alpine samples (8.3%). The low classification rate was related to the low number of samples.  5%). The more precise and efficient separation between the four regions was obtained when taking into account the season (summer-winter). In this case, the classification rate was > 90%. Figure 3c shows the DA plot for 39 milk samples obtained in winter 2013. The analysis generated four well-separated groups. Separation of the regions was also achieved based on FA composition when the data were treated altogether. For F1 (89.0%), C15:0 and its isomers were the main discriminatory parameters, while for F2, the most discriminatory parameters (10.5%) were C18:3-6 and CLA. However, cross-validation was low, and only 35.9% of the samples were classified correctly. The highest rate of classification was for the Mediterranean samples (66.7%) and the lowest for the Alpine samples (8.3%). The low classification rate was related to the low number of samples.

Stable Carbon Isotope Composition of FA in Milk
This study used the stable carbon isotope composition of individual FAs to check the variability in production year and season in 2013 and 2014. The differences measured in the isotopic composition between individual milk FA reflected diet regime, the isotopic fractionation during biosynthesis, metabolic rates, and isotope turnover time. Table 2 presents the stable carbon isotope values of individual FAs in cow milk samples and stable carbon isotope values in bulk freeze-dry milk. The given values were the mean and standard deviation (SD). Results were shown as box plots of δ 13 C values of milk fatty acids and bulk milk samples (Figure 4). Data points, exceeding two units of SD, were regarded as outliers, and in the cases where measurements were below the 5% confidence level, the data were discarded. Interestingly, FAs formed during biosynthesis had the lowest δ 13 C values due to isotope fractionation. compared to C16:0, which originates in part from dietary lipids and to a significant extent from de novo synthesis from acetate, derived mainly from dietary carbohydrates [36]. The  13 C value (δ 13 C-18:0−δ 13 C-16:0) in the samples were < −3.3‰, which is typical for ruminant fat.
Stable carbon isotope composition of individual FA indicated less variability between season and year than fatty acid composition, which suggested a more distant relationship. There was also no data supporting the correlation between a time-significant percentage of FA and δ 13 C values within the same FA. The results of a multivariate analysis of the stable isotope content of FA are presented in Figure 4. The temporal separation between selected groups described 93.8% of the variability. The most important stable isotope parameters that contributed to this separation were: δ 13 C-16:0, δ 13 C-18:1-9c, δ 13 C-bulk, δ 13 C-14:1, and δ 13 C-18:0, and a leave-one-out cross-validation (LOOCV) classified 90.9% of the samples correctly.   Table S1. The δ 13 C in FA in milk collected from the Mediterranean region differed significantly (p < 0.05) from the other three regions. The differences were more apparent in summer than winter. The Mediterranean region also had more negative δ 13 C values, which indicated that cows had a higher In all the samples, C15:0 had significantly lower average δ 13 C values (−30.3 ± 3.7% ) compared to other FAs. The highest average δ 13 C value (−22.6 ± 2.9% ) was observed in C14:0. Higher δ 13 C values were found in winter, which reflected the higher amounts of corn in the diet. Although the reported δ 13 C values from specific regions, e.g., the Mediterranean region, were influenced by the year of production, significant differences in reported δ 13 C values between the summer and the winter season were observed. The fatty acid C18:1ω-9 was an exception with consistently lower values during the summer. It was assumed that this variability was a result of different production practices and the amount of grass and corn forage in a cow's diet. The C18:1ω-9 in milk from cows fed a C 3 diet is enriched in 13 C relative to C18:0, whereas no enrichment occurs when cows' are fed on a C 4 diet [17]. In this study, an average enrichment of 1.6% was observed in milk sampled during the summer and supported the assumption that a C 3 -based diet could result in the transfer of higher amounts of 18:1 ω-9 directly from the feed. As part of a dietary source, C18:0 can be produced by microbial biohydrogenation of C18:0 as the final product of ruminal biohydrogenation [35]. Microbial biohydrogenation can significantly affect the isotopic values of C18:0, leading to a lower δ 13 C values compared to C16:0, which originates in part from dietary lipids and to a significant extent from de novo synthesis from acetate, derived mainly from dietary carbohydrates [36]. The ∆ 13 C value (δ 13 C-18:0−δ 13 C-16:0) in the samples were <−3.3% , which is typical for ruminant fat.
Stable carbon isotope composition of individual FA indicated less variability between season and year than fatty acid composition, which suggested a more distant relationship. There was also no data supporting the correlation between a time-significant percentage of FA and δ 13 C values within the same FA. The results of a multivariate analysis of the stable isotope content of FA are presented in Figure 4. The temporal separation between selected groups described 93.8% of the variability. The most important stable isotope parameters that contributed to this separation were: δ 13 C-16:0, δ 13 C-18:1ω-9c, δ 13 C-bulk, δ 13 C-14:1, and δ 13 C-18:0, and a leave-one-out cross-validation (LOOCV) classified 90.9% of the samples correctly. Figure 4 shows the summer and winter δ 13 C values of different FAs in milk for the four studied regions. The δ 13 C values of individual fatty acids (season, year, and production region) are given in Table S1. The δ 13 C in FA in milk collected from the Mediterranean region differed significantly (p < 0.05) from the other three regions. The differences were more apparent in summer than winter. The Mediterranean region also had more negative δ 13 C values, which indicated that cows had a higher amount of grass-based (C 3 plant) forage in their diet, although the higher temperatures in the Mediterranean region could be a factor. The highest δ 13 C values were in FA in milk from the Pannonian region. Besides, the δ 13 C-18:1ω-9c determined in the Pannonia region differed significantly from Alpine and Dinaric regions.

Quality of Milk
Some comparisons can be drawn about the ratios of SFA, MUFA, and PUFA in milk. A healthy diet should contain more PUFA and less of SFA since SFAs are a risk factor for cardiovascular diseases [37]. Alternatively, MUFA and PUFA, except for trans unsaturated fatty acids, reduce the level of low-density lipoprotein (LDL) cholesterol and increase the level of high-density lipoprotein (HDL) cholesterol [38].
The lowest percentage of SFA (avg. 66.97 ± 2.54 % of total FA) was observed in milk produced in summer in the Alpine region, for all three sampling years, while in winter, there were no significant differences between the regions in the SFA percentage of milk sampled in the summer ( Table 3). The feeding regime in the Alpine region is based on the managed grazing of pastures during the summer, while in the winter, the cows' diet consists of a higher amount of cereal concentrates and conserved forages. The feeding of grass conservation products (hay or silage) increases the SFA content of milk fat.   The consensus is that PUFAs, especially ω-3, ω-6, and CLA, have a beneficial effect on human health. Besides, the ω-6:ω-3 ratio, expressed as the concentration of linoleic acid (C18:2ω-6) to α-linolenic (C18:3ω-3) acid, i.e., the most abundant ω-6 and ω-3 FA, as well as atherogenic (AI) and thrombogenic (TI) indices are used to evaluate the nutritional value of milk fat since the milk fat with low AI and TI values could provide protection against coronary heart disease. Further, ω-3 PUFA reduces the risk of cardiovascular disease, hypertension, cancer and has been linked to improved neurological function [5,6,10,13]. The most abundant PUFA in fresh forages is ω-3 (especially α-linolenic acid). In this study, the highest percentage of α-linolenic acid was in the Alpine region in the summer (avg. 0.87 ± 0.38% of total FA) and in the Mediterranean region in winter (avg. 0.71 ± 0.2% of total FA). While the milk from the Pannonian region had the lowest percentage of α-linolenic acid in both summer and winter (avg. 0.62 ± 0.22% of total FA in summer and 0.46 ± 0.07% of total FA in winter). The lower levels of α-linolenic in winter are the result of oxidative losses of PUFA during the wilting of grasses [13]. A correlation also existed between α-linolenic acid and ω-3 FA. The highest percentage of ω-3 FA was present in milk from the Alpine region in summer (avg. 0.97 ± 0.38% of total FA), while in winter, the highest levels of ω-3 FA were observed in the Mediterranean region (avg. 0.83 ± 2.1% of total FA). The lowest percentage of ω-3 FAs was found in milk collected from the Pannonian region in all three production years (avg. 0.76 ± 0.29% of total FA in summer and avg. 0.57 ± 0.09% of total FA in winter). Changes in the ω-3 FA percentage were smaller in the Pannonian region, during both summer and winter compared to the other studied regions. The reason was that the cows' were fed a diet of corn silage (intensive milk production) and conserve forage throughout the year. Low values of ω-3 FA in low altitude regions of Slovenia (Pannonian and Dinaric regions) were consistent with the literature data for lowlands [39,40]. Since a high ω-6:ω-3 ratio can be detrimental to human health and a factor in cardiovascular disease, cancer, inflammatory and autoimmune disease, an optimum ratio in the human diet of 1:1 ratio is considered, while a ratio <5:1 is recommended [39]. However, a higher ω-6:ω-3 ratio is indicative of higher amounts of cereals (corn, soybean, barley, and oats) in the cows' diet. The ω-6:ω-3 ratio in Slovenian milk is between 4.2:1 and 2.5:1 and is lower in summer compared to winter. The lowest ω-6:ω-3 ratio was in milk produced in the Mediterranean region through the whole season and did not exceed 2.6:1, which was indicative of a fresh forage-based diet. In the Dinaric and Pannonian regions, the effect of feeding cows corn silage and cereals throughout the year increased the ω-6:ω-3 FA ratio from 3.2:1 to 4.2:1. The data of AI and HI indices were comparable to the study performed by Ferlay et al. [41] determined in the milk of Tarentaise and Montbéliarde cows (3.14 and 3.43, respectively), but they were lower than the data reported by Nantapo et al. [42] in the milk of Friesian, Jersey, and Friesian x Jersey cows with values ranged from 4.08 to 5.13. The highest AI indices in our study were determined in the Mediterranean region (2.89-4.13) and the lowest in the Alpine region (2.31-3.79).
Environmental factors related to year and season together with diet are the most important factors influencing the composition of CLA in milk. In this study, the percentage of CLA ranged from 0.52 ± 0.18% of total FA to 0.73 ± 0.31% of total FA in milk from all four studied regions and increased with dietary intake of 18-carbon PUFAs. Increasing the amount of PUFAs in milk could be achieved by feeding cows' a diet based on fresh silage. Milk collected in summer in the Alpine region (avg. 0.99 ± 0.28% of total FA) had the highest percentage of CLA, most likely from grazing on pasture. In the winter, there was no significant difference between the Alpine and other regions since the animal's diet is composed of both corn and grass silage. Because Alpine and Mediterranean milk in summer contained a higher percentage of ω-3 PUFA and CLA, compounds known to have positive effects on human health [40,43], this milk is likely "healthier" than milk produced in other regions in Slovenia. However, despite the benefits of feeding cows fresh grass silage, the overall effects on the FA composition of milk would not improve human health over the long-term [44].

Extraction and Esterification of Milk FAs
Frozen milk samples were thawed (37 • C) and homogenized. Milk fat consists predominantly of triacylglycerides, which were extracted using the in-situ-trans-esterification method without prior extraction of the fat from samples according to the method of Park and Goins (1994). The total lipids fraction was extracted from the milk samples (450 µL) with dichloromethane (300 µL) and 0.5 M sodium hydroxide in methanol (3 mL), purged with nitrogen and heated at 90 • C for 10 min, after which the samples were cooled rapidly. Fatty acid methyl esters (FAMEs) were formed by adding 14% BF 3 -methanol solution (3 mL) and heating the sample at 90 • C for a further 10 min. After cooling, the fatty acids methyl esters (FAMEs) were extracted with hexane (1.5 mL) and transferred to a GC vial and stored at −20 • C.

Gas Chromatography of Fatty Acids
The characterization of FAMEs was performed using Agilent 6890N (Network GC System, Agilent Technology, Santa Clara, CA, USA) gas chromatography with FID detector (GC-FID). The separation was achieved on an Omegawax 320, 30 m × 0.32 mm, × 0.25 µm, on a capillary column (Supelco, Bellefonte, PA, USA). The temperature program was as follows: 50 • C (held 2 min) to 220 • C at 4 • C/min (held 20 min). The carrier gas was helium in constant flow mode (1 mL/min). The injector was set to 260 • C and the detector to 280 • C. The individual fatty acids were identified and quantified by comparing their retention times with those of a standard Supelco 37 component FAME Mix (Supelco) and expressed in weight percent of total fatty acids. Procedural blanks were analyzed with each set of samples. Standard Supelco 37 component FAME mix was analyzed after every ten samples to verify the stability of the analytical system. Method precision based on measurements of replicate (n = 2) real samples was 5%. Results were expressed as the percentage of each fatty acid with respect to the total fatty acids [45]. This percentage was calculated using the peak area of the samples corrected with the respective correction factors, according to Christie and Han [46].

Elemental Analysis-Isotope Ratio Mass Spectrometry (EA-IRMS)
The stable carbon isotope ratio measurements were reported in δ notation expressed in per mile (% ) relative to the Vienna-Pee Dee Belemnite (V-PDB) standard following Equation (1) [47]: where superscripts i and j denote the highest and the lowest atomic mass number of element E, and R P and R Ref is the ratio between the heavier and the lighter isotopes (in case of carbon, 13 C/ 12 C) in the sample (P) and reference material (Ref ).
To determine the bulk carbon isotope composition, approximately 1 mg of freeze-dried milk was weighed into a tin capsule, closed with tweezers, and placed into the automatic sampler of the elemental analyzer. All analyses were performed separately on a Europa Scientific 20-20 continuous flow mass spectrometer with an ANCA-SL solid-liquid preparation module (Sercon, Crewe, UK). Analyses were normalized against international standards: IAEA-CH-7 and IAEA-CH-6 with δ 13 C values of −32.15 ± 0.04% and −10.45 ± 0.04% , respectively. For independent control, laboratory reference material IAEA-CRP 2013 with δ 13 C value −20.3 ± 0.09% and UreaC with δ 13 C value −28.3 ± 0.2% were used. The analytical precision, expressed as the standard deviation of control material, was ± 0.2% . The 13 C/ 12 C ratio of synthetic fatty acid standard C19:0 (methyl nanodecanoate, RESTEK) was determined using a Europa Scientific 20-20 continuous flow mass spectrometer with an ANCA-SL solid-liquid preparation module (Sercon, Crewe, UK). Results were further normalized against two international reference standards: IAEA-CH-7 and IAEA-CH-6. For independent control, the international reference material USGS40 with δ 13 C values of −26.39 ± 0.04% was used. Analytical precision, expressed as the standard deviation of the control materials, was ± 0.2% .

Compound Specific Isotope Analysis of δ 13 C in Fatty Acids
The Isotopic compositions of individual FA were determined using an Agilent 6890N GC-C system coupled to an IsoPrime GV IRMS (GV Instruments, Manchester, UK)). The separation was achieved using a DB-1MS (60 m × 0.32 mm × 0.25 µm) capillary column (Agilent J&W, USA). The temperature program was as follows: 120 • C (held 1 min) to 300 • C (held 20 min) at 3 • C/min. The carrier gas was helium in the constant flow mode (1 mL/min). The injector temperature was 260 • C, and the oxidation reactor (Cu/O) in the 6890N GC/C system was 900 • C. Peak identification of individual fatty acids was performed by comparing their retention times with those in the standard Supelco 37 component FAME Mix, which was analyzed in the measurement sequence. Instrument stability was ≤0.08% . C19:0 FAME (methyl nanodecanoate) with δ 13 C values of −30 ± 0.1% (prior determined on EA-IRMS) was used as a laboratory reference material. For data normalization, a single-point normalization method was used [48]. The precision of measurements expressed as the standard deviation of the laboratory reference material was 0.3 to 0.5% .

Statistical Evaluation
For statistical evaluation, FAs were considered individually and grouped in classes as follows: saturated fatty acids (SFA), unsaturated fatty acids (UFA), MUFA, ω-3 PUFA, ω-6 PUFA, and CLA. Statistical calculations and multivariate analysis were carried out using the XLSTAT software package (2020.1.3, Addinsoft, New York, NY, USA). The statistical evaluation was performed using R, where non-parametric tests were used, i.e., Kruskal-Wallis followed by a post hoc test. In the case of pairwise comparison, a Mann-Whitney test was used. Basic statistics included mean values (median and arithmetic mean), standard deviation (SD), minimum and maximum, while for multivariate analysis, discriminant analysis (DA) was used.

Conclusions
Slovenian milk samples had been classified according to the year, season, and region of production based on the determination of the percentage and stable isotope composition of individual fatty acids. It was found that the FA composition differed significantly, depending more on year and season of production than geographical region. A leave-one-out cross-validation (LOOCV) classified 86.9% of the samples correctly. For geographical origin determination, it was found that the Mediterranean group was separated from other regions based on a high percentage of C15:0 and C17:0 and their isomers. The cross-validation test showed 64.1% correct classification with the highest rate of classification for the Pannonian (80.8%) samples. Although the classification rate according to the year/season of production improved when stable isotope composition of individual FA was included, describing 93.8% of variability, it was found that δ 13 C in FA could be more useful for studying metabolic transformation processes in ruminants. It appeared that the natural δ 13 C differences in FA depended on both-the diet type and the precursor FA. The "healthiest" milk for human consumption, based on FA composition, was produced in the Alpine and Mediterranean regions. The quality of milk was also affected by summer/winter season production. The results gave proof of the concept that FA composition could help to discriminate the geographical origin of milk and thus identify products of better quality. Finally, to confirm the validity of the markers, future studies should focus on how factors, such as breed and age, lactation period, production level, the farming system, influence the FA and stable isotope composition of milk.
Supplementary Materials: The following are available online, Table S1: δ 13 C of fatty acids in milk from a different geographical location in Slovenia during summer/winter in 2013 and 2014; Table S2: Geographical information of the sampling location.

Funding:
The work was performed within the project V4-1108 entitled "The use of specific methods for determination and prevention of adulteration of milk and dairy products", financially supported by Slovenian Research Agency and Ministry of Agriculture and the Environment, and IAEA project "The use of stable isotopes and elemental composition for determination of authenticity and geographical origin of milk and dairy products" (Contract No. 17897). This research represents a part of the ERA Chair ISO-FOOD-for isotope techniques in food quality, safety, and traceability (FP7, GA no. 621329) and MASSTWIN-Spreading excellence and widening participation in support of mass spectrometry and related techniques in health, the environment and food analysis (H2020, GA no. 692241).