Stable Isotope and Element Profiling for Determining the Agroclimatic Origin of Cow Milk within a Tropical Country

Information on the geographic origin of milk is important in determining quality attributes and for economic gain through building brand value associated with origin. Stable isotope signatures and trace element concentrations have been increasingly used in authentication of milk, though information on the power of such technology in verifying agroclimatic origin in small continents with diverse climatic, environmental conditions, and animal management practice is scarce. Therefore, the main objective of this study was to investigate the possibility of using a stable isotope composition of C, N, O, and H and element fingerprints to determine the agroclimatic origin of milk produced in different agroclimatic zones of Sri Lanka. Stable isotopes ratios of C, N, H, and O, and elemental fingerprints of milk samples were determined by IRMS and ICP-MS, respectively. Significant variations were observed in stable isotope ratios, especially δ18O and the mean content of Li, Al, Cr, Mn, and Sr in the bulk milk samples obtained from different agroclimatic zones. A linear discriminant analysis differentiated cow milk produced from four agroclimatic zones based on stable isotope ratios, and the inclusion of elemental ratios enhanced the discriminating ability.


Introduction
Authentication of geographical origin is a rapidly emerging topic in food safety, quality, traceability, and consumer protection. Tools to validate the geographic origin of food including milk are extremely valuable, as the origin is often associated with quality attributes, for instance, free from disease and pollution. Recent scares, such as melamine adulteration of milk powder [1], have focused wider public attention on the importance of knowing the origin of milk.
Various analytical techniques, including liquid and gas chromatography, isotope ratio and elemental analyses, spectroscopy, DNA, and sensor techniques, have been widely used to authenticate the geographical origin of various foodstuffs [2]. Among these techniques, trace element concentrations and stable isotope signatures have been increasingly applied in authentication studies, since their chemical patterns have a direct causal relationship to geographical factors, such as the environment in which the products were grown and/or processed. Stable isotope signatures mainly depend on climatic or geographical conditions, while the elemental composition is primarily affected by geological and pedological Table 1. Milk production zones in Sri Lanka (after [23]

Sample Collection and Preparation
Milk samples were collected from four different agroclimatic zones in Sri Lanka, including dry zone, coconut triangle, mid country, and up country ( Figure 1). Six large to medium-scale dairy farms were selected from each of the agroclimatic zones for sampling except the dry zone. Seven farms/locations were selected in the dry zone as a contingency plan since one farm indicated that they were planning to convert their dairy farm to a buffalo farm. Two liquid milk samples of 250 mL were collected on three consecutive days from a selected farm for trace element and stable isotope analysis. The same farms were sampled during wet and dry periods. Samples were collected from the bulk collection tank immediately after milking into pre-cleaned (acid washed) high-density polyethylene bottles. Samples were transported in cooling containers with ice and then stored in a −20 • C freezer until analyses. Table 1. Milk production zones in Sri Lanka (after [23]

Sample Collection and Preparation
Milk samples were collected from four different agroclimatic zones in Sri Lanka, including dry zone, coconut triangle, mid country, and up country ( Figure 1). Six large to medium-scale dairy farms were selected from each of the agroclimatic zones for sampling except the dry zone. Seven farms/locations were selected in the dry zone as a contingency plan since one farm indicated that they were planning to convert their dairy farm to a buffalo farm. Two liquid milk samples of 250 mL were collected on three consecutive days from a selected farm for trace element and stable isotope analysis. The same farms were sampled during wet and dry periods. Samples were collected from the bulk collection tank immediately after milking into pre-cleaned (acid washed) high-density polyethylene bottles. Samples were transported in cooling containers with ice and then stored in a −20 °C freezer until analyses.  Milk samples were prepared based on the method proposed by Feng et al. [24] and validated by Luna et al. [25]. An aliquot of 1.5 mL of milk was centrifuged at 20,000× g for 30 min at 4 • C. With the aid of a pipette, without disturbing the fat layer the skim milk fraction in the bottom of the tube was transferred into a micro-centrifuge tube. The initial tube with milk fat was saved to analyze stable isotope ratios of milk fat. Then 6.0 M HCl (about 60 µL) was added to the skim milk sample to facilitate precipitation of milk casein at its isoelectric point (pH 4.6). Next, the micro-centrifuge tube was vortexed for 5 min and allowed to settle for 10 min at room temperature to facilitate casein precipitation. Finally, the tube was centrifuged at 20,000× g for 30 min at 4 • C to precipitate milk casein. The supernatant containing the whey fraction was transferred to another tube to separate the milk whey fraction. All these fractions (fat, casein, and milk whey) of the samples were stored at −20 • C until lyophilization.

Stable Isotope Analysis
For δ 13 C and δ 15 N analysis, whole milk, milk fat, milk casein, and whey were measured as 1.0 ± 0.05 mg, 0.5 ± 0.05 mg, 1.0 ± 0.05 mg, and 2.0 ± 0.05 mg into tin capsules, respectively, and encapsulated. The encapsulated samples were placed in a desiccator until the isotope ratio mass spectrometer (IRMS) analysis was performed. Bulk δ 13 C and δ 15 N were measured by combustion to CO 2 and N 2 , respectively, using Euro-EA Elemental Analyzer (EuroVector). The isotope ratio of the gases was measured on a Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK) operating in continuous flow mode. The isotopic compositions of the sample gases were normalized and reported against the international scales for carbon and nitrogen, VPDB, and ambient air. Normalization was made by three-point calibration with international standards, USGS 40 (L-glutamic acid, δ 13 C = 26.24‰ and δ 15 N = −4.52‰), USGS 41 (L-glutamic acid, δ 13 C = 37.76‰ and δ 15 N = 47.56‰), and EDTA-OAS laboratory standard (Elemental Microanalysis Ltd., Cornwall, UK, δ 13 C = −38.93‰ and δ 15 N = −0.73‰). Time-based drift correction was calculated by placing laboratory standards (IAEA 153 milk powder) with EDTA-OAS after every 12 samples. Precision was assessed by placing sequential duplicates of the samples after every 10th sample.
For δ 18 O and δ 2 H analysis, 0.6 ± 0.05 mg of bulk milk samples were weighed into silver capsules which were closed loosely to allow the exchange of water vapor for correction of the effect of exchangeable hydrogen. Capsule bottoms were slightly flattened to keep them standing upright. Samples were left for six days in ambient conditions to achieve equilibrium of exchangeable hydrogen with the lab atmospheric water. Three replicates were prepared for each sample.
Samples were placed into the autosampler carousel. Spaces were left for polymer standard materials (IAEA-CH7, Kilwell, Nylex, Shah Alam, Malaysia). The carousel was placed in a desiccator and evacuated to remove air. Samples were vacuum dried on a hot plate set to 40 • C for four days. IAEA Milk Powder 153 (δ 2 H = −40.78 ± 5.2‰, δ 18 O = 16.81 ± 0.59‰), USGS-53 (Lake Shala Distilled Water, δ 2 H = +40.2 ± 0.4‰, δ 18 O = +5.47 ± 0.03‰), USGS 47 (Lake Louise Drinking Water, δ 2 H = −150.2 ± 0.5‰, δ 18 O = −19.80 ± 0.02‰), CH7 (δ 2 H = −99.2‰), and IAEA601 (Benzoic Acid, δ 18 O = 23.14‰) were used as the reference samples. CH7 and IAEA601 were repeated as control samples at each 15th position throughout the batch for drift correction by linear regression against time. Control samples and USGS-53 were prepared on the same day in a separate tray. Hydrogen and oxygen isotope ratios of the samples were determined by pyrolysis to H 2 and CO gases in a Thermal Conversion Elemental Analyzer (TC/EA, Thermo Scientific, Bremen, Germany), followed by isotope analysis of the gas using Thermo Delta-V IRMS in continuous flow mode. Raw results were corrected to the international VSMOW (Vienna Standard Mean Ocean Water) isotope scale using two-point calibration. The measurement uncertainty was calculated by analyzing six to eight replicates of the quality control standard interspersed regular intervals (after every 12 samples) in every batch of samples analyzed throughout the sample analysis period. The measurement uncertainty reported herein is the highest uncertainty value obtained for a single batch of samples throughout the entire sample analysis. The uncertainty values were ±0.18‰ for δ 13 C, ±0.25‰ for δ 15 N, ±0.32‰ for δ 18 O, and ±1.46‰ for δ 2 H. The isotopic difference between the sample and an international standard reported in delta notation (δ) units in parts per thousand (‰) is determined by the formula [26]: where i E denotes the higher and j E the lower atomic mass number of element E. The subscript P denotes the samples used to determine the respective values. R indicates the isotope number ratios and Ref indicates the reference material; Vienna Pee Dee Belemnite (V-PDB) for carbon (δ 13 C), atmospheric N 2 (Air-N 2 ) for nitrogen (δ 15 N), and VSMOW for δ 18 O and δ 2 H.

Trace Element Analysis
An aliquot of 1 mL of each milk sample was digested with 5 mL of supra-pure HNO 3 (>69%, Fluka, Switzerland) and 1 mL of H 2 O 2 (35%, Sigma-Aldrich, Darmstadt, Germany) at 200 • C for 10 min in a MARS-6 Microwave digester (CEM Mathews, Charlotte, NC, USA) equipped with an EasyPrep Plus high-pressure vessel system. The resulting solutions were then diluted to 25 mL with deionized water (<0.055 µScm −1 ) and stored in polypropylene bottles at 4 • C until analysis.
Digested samples were filtered through 0.45 µm cellulose acetate membrane filters, and concentrations of trace elements were measured using Thermo ICapQ high-resolution ICP-MS (Thermo Fisher, Bremen, Germany). Before the analysis, the instrument was optimized with a solution containing Ba, Bi, Ce, Co, In, Li, and U to cover the entire mass range and obtain high sensitivity. A series of multi-element standards prepared by diluting 10 mg L −1 AccuStandard ® (New Haven, CT, USA) stock solution was used to calibrate the instrument. Internal calibration was carried out using 100 µg L −1 of 103 Rh and 185 Re standards. All measurements were conducted at least in triplicate, and analytical results were expressed as mean values. Analytical quality controlling was assessed using Certified Reference Material BCR-063R (skim milk powder from the Institute for Reference Materials and Measurements, European Commission). The given reference values and test values showed that the recoveries of available elements ranged from 85-115%, and relative standard deviations by sextuplicate analyses of elements were less than 10%. Further, synthetic external standards, duplicates, and reagent blanks were analyzed at every 10th sample and measurement uncertainties were always less than 2%. Deionized water blanks were also analyzed recurrently to check for sample carry-over.

Statistical Analysis
Results of isotope ratios and elemental concentrations were statistically evaluated using R-Studio version 3.3.0 and SAS (version 9.1; SAS Institute Inc., Cary, NC, USA). Data were preliminarily explored by preparing box plots and scatter diagrams. Descriptive statistics were taken for each data set. One-way analysis of variance (ANOVA) was used to assess stable isotope data to evaluate whether the variables used differ between the agroclimatic zones, the sampling period's effect, and their interaction. The means were separated using Duncan's multiple range tests at a confidence level of 95%. Principle component analysis (PCA) was performed to check if the combined isotopic parameters could be used to distinguish different agroclimatic zones. The linear discriminant analysis (LDA) test was performed with a p-value of ≤0.05 to discriminate samples according to agro-climatic zones based on the elemental composition and isotope ratios as well as isotopic data alone. The predictive ability of the LDA model was evaluated by leave-one-out cross-validation.

Stable Isotopic Composition in the Milk and Its Constituents
Stable isotope ratios of 142 milk samples collected from four different agroclimatic zones of Sri Lanka varied from −15.5 to −28.0‰ for δ 13 C, 11.6 to 3.3‰ for δ 15 N, −48 to −155‰ for δ 2 H, and 26.0 to 16.7‰ for δ 18 O. The mean δ 18 O of whole milk samples showed a significant difference among the four agroclimatic zones ( Table 2). Milk from the dry zone region showed relatively higher δ 18 O values than the other zones, while up country milk showed the lowest values ( Figure 2). The present study showed a significant correlation of δ 18 O values with latitude, a moderate correlation with altitude (Figure 3), and a lower correlation with distance to the coast. It has been revealed that the isotopic enrichment of exchangeable oxygen and hydrogen in plants and animals is usually shifted in relation to the local precipitation, as the enrichment of all sources of absorbed water is altered by plant or animal metabolism. Temperature is an important factor governing the enrichment of δ 18 O and δ 2 H in precipitation [27]. The milk samples collected from Sri Lanka showed higher δ 18 O values at lower elevation and lower values at high elevation. This observation could have resulted from the influence of climatic factors, such as temperature and rainfall, as relative humidity affects the fractionation of oxygen isotopes. Each agroclimatic zone receives rainfall from different rainfall regimes by the two main monsoon rains and two convectional (inter-monsoon) rains known to have characteristic isotopic signatures [28]. The isotopic composition of monsoon rains in Sri Lanka is affected by terrestrial and oceanic moisture sources, depression, and cyclonic conditions of the Bay of Bengal, Arabian Sea, and the country's topography [29]. The present study showed a significant correlation of δ 18 O values with latitude, a moderate correlation with altitude (Figure 3), and a lower correlation with distance to the coast. It has been revealed that the isotopic enrichment of exchangeable oxygen and hydrogen in plants and animals is usually shifted in relation to the local precipitation, as the enrichment of all sources of absorbed water is altered by plant or animal metabolism. Temperature is an important factor governing the enrichment of δ 18 O and δ 2 H in precipitation [27]. The milk samples collected from Sri Lanka showed higher δ 18 O values at lower elevation and lower values at high elevation. This observation could have resulted from the influence of climatic factors, such as temperature and rainfall, as relative humidity affects the fractionation of oxygen isotopes. Each agroclimatic zone receives rainfall from different rainfall regimes by the two main monsoon rains and two convectional (intermonsoon) rains known to have characteristic isotopic signatures [28]. The isotopic composition of monsoon rains in Sri Lanka is affected by terrestrial and oceanic moisture sources, depression, and cyclonic conditions of the Bay of Bengal, Arabian Sea, and the country's topography [29]. Higher evapotranspiration rates in the dry zone may also lead to enriched δ 18 O values in plants and animals. In Sri Lanka, the temperature decreases at a steady rate of about 6.5 • C for each 1000-m rise. Higher temperatures are experienced generally in the dry zone areas of northern, northcentral, and eastern regions of the island that range between 33.3-34.7 • C [30]. The dry zone receives a mean annual rainfall of less than 1000 mm with a distinct dry season, while the intermediate zone receives a mean annual rainfall of about 1750 mm. The wet zone receives over 2500 mm of annual rainfall [30].
In the present study, the δ 15 N values of whole milk, milk casein, and whey revealed a moderate potential for discriminating cattle farming zones. However, this study found that δ 15 N as well as δ 13 C values of dry zone samples were significantly different from other zones ( Table 2). The highest δ 15 N values were reported in the dry zone samples, and the lowest δ 15 N values were reported in the up country. The whole milk δ 15 N values were between their casein and whey δ 15 N values for all the agroclimatic zones except the upcountry. In up country, the whole milk δ 15 N values were higher than their respective casein and whey δ 15 N values. Further, up country recorded the lowest temperature range (18-23 • C) and the other zones showed a higher temperature range (27-34 • C) during sampling. According to literature, ambient temperature could influence casein δ 15 N values [17]. The isotopic fractionation of nitrogen in dairy cows was estimated to be high in metabolic origin and affected by temperature as the metabolic rates of animals are indirectly affected by environmental temperatures. In up country, animals are subjected to minimal heat stress compared to the other zones.
Complete milk protein is composed of 20% whey protein and 80% casein. δ 15 N values of milk whey were more depleted than that of milk casein and whole milk in all four agroclimatic regions. A similar observation has been reported by Kornexl et al. [31]. This finding may be attributed to the different amino acid compositions and relative abundance of amino acids originating from body tissues or diet [31,32]. Further, the δ 15 N value in agricultural products is influenced by the differences in soil conditions, amount of rainfall, soil age, nitrogen turnover rate in soils, and other complex parameters [17,31]. Regional differences in the δ 15 N value of milk may reflect the differences in feeding regimes depending on the regional dairy farm [11]. There is a trend for increasing foliar δ 15 N with decreasing rainfall [33].  Higher evapotranspiration rates in the dry zone may also lead to enriched δ 18 O values in plants and animals. In Sri Lanka, the temperature decreases at a steady rate of about 6.5 °C for each 1000-m rise. Higher temperatures are experienced generally in the dry zone areas of northern, northcentral, and eastern regions of the island that range between 33.3 °C-34.7 °C [30]. The dry zone receives a mean annual rainfall of less than 1000 mm with a The difference in δ 13 C reflects the relative amounts of C3 and C4 plants used by the regional feeding regimes. According to the Calvin cycles in photosynthetic CO 2 fixation, C3 plants show δ 13 C values between −24 and −32‰. Plants following the C4-dicarboxylic acid pathway show δ 13 C values between −10 and −19‰ [34]. Except in the dry zone, all other cattle management zones used particularly the grasses of the genera Brachiaria and Panicum to feed cattle, which fall into the C4 plant category. Gliricidia sepium, which is in the C3 category, was fed to cattle of all the zones except the up country. In up country, Ryegrass (Lolium multiflorum) was given to the animals which is also a C3 plant. However, a considerable amount of concentrate feed, such as rice bran, rice polish, coconut poonac were given as cattle feed in most of the farms. Cattle reared in dry zone eat various types of plants, while freely grazing on communal grazing lands. Therefore, the δ 13 C values showed a value in between C3 and C4 plants in cattle feed for most of the samples. However, in this study, the δ 2 H value displayed weak potential to differentiate milk samples of different agroclimatic zones. The variations in δ 18 O, δ 13 C, and δ 15 N values for each of the agroclimatic zone investigated are presented in Figure 4. This representation helps visualize the potential of stable isotope values as a 'screening' parameter to differentiate milk produced in different agroclimatic zones. Samples from the dry zone could be visualized separately from the mid country and up country.

Influence of Sampling Period on Stable Isotope Values of Milk
Although δ 13 C, δ 15 N, δ 18 O, and δ 2 H in whole milk and milk components were significantly influenced by agro-climatic zone, these variables were also moderately influenced by the sampling period (Table 3)

Influence of Sampling Period on Stable Isotope Values of Milk
Although δ 13 C, δ 15 N, δ 18 O, and δ 2 H in whole milk and milk components were significantly influenced by agro-climatic zone, these variables were also moderately influenced by the sampling period (Table 3) and less influenced by interactions of the sampling period and the region. Mean δ 18 O values in whole milk displayed a significant difference in wet and dry sampling periods for each agroclimatic zone except for the dry zone (p < 0.001). The only variable that displayed a significant difference in wet and dry periods for the dry zone was δ 15 N of milk casein (p < 0.001). Significant variation between sampling periods in the δ 15 N values of milk casein and whole milk samples collected from coconut triangle and up country was observed. Mean δ 15 N of whole milk and milk fractions and δ 18 O showed relatively enriched values for the wet period compared with the dry period. Forage production and availability in Sri Lanka mainly depend on the rainfall. The availability of lush pastures is high during rainy season and forage is scanty during dry season. Due to the scarcity of forage during the dry season, cattle are often fed with different feeding materials, such as silage and concentrate feeds. The difference of δ 15 N during wet and dry seasons may have resulted from variations in ambient temperature and feed material [17,31].
The topographical features and regional scale monsoons wind regimes dominate the climate dynamics in each agroclimatic zone. Monsoonal, convectional, and depressional rain contribute to the seasonal rainfall in Sri Lanka. The topographical features of the various agroclimatic zones mean they receive rain in different months of the year. The two main monsoonal rains and two inter-monsoon rains have characteristic isotopic signatures [28]. The up country agroclimatic zone falls in the most affected area of the southwest monsoon. North-east monsoon and second inter-monsoon rains are the significant rainy seasons in mid country. The dry zone and coconut triangle receive abundant rainfall from the northeast monsoon [30]. The source of precipitation, temperature, isotopic depletion of moisture in clouds moving inland, elevation effect, and cooling of the vapor masses as they rise over the landscape may change the isotopic signatures of δ 18 O and δ 2 H in different months of the year in agroclimatic zones.

Determination of Geographical Origin by Statistical Analysis
The combined isotope parameters obtained from milk samples were evaluated to distinguish agroclimatic zones using multivariate statistics. The first two principal components explained 69% of the total variance (Table 4). However, clear clustering was not observed in the score plot of PC1 and PC2 ( Figure 5). As the PCA could not clearly group the samples based on the origin, LDA was employed with backward stepwise refinement of stable isotope ratios to provide a mathematical model to classify samples simultaneously according to their agroclimatic origin. The selected variables in the model were whole milk δ 15 N, milk casein δ 15 N, milk whey δ 15 N, whole milk δ 13 C, milk fat δ 13 C, milk casein δ 13 C, whole milk δ 2 H, and whole milk δ 18 O. The efficacy of the model in discriminating milk samples originated in coconut triangle, dry zone, mid country, and up country from the rest agroclimatic regions were 83%, 92%, 89%, and 97%, respectively. After cross-validation, the above percentages were 75%, 92%, 81%, and 91%, respectively. Table 5 summarized original and cross-validated classification by discriminant analysis for the samples collected from four different agroclimatic zones. The discriminant analysis also showed a degree of overlap. As illustrated in Figure 1, farms located closer to the border areas of these agroclimatic zones share similar climatic conditions. Samples collected during the dry period were evaluated for elemental composition and stable isotope ratios. Previous studies have shown that there is a correlation between the elemental profile of animal products and geographic origin. Further, the elemental profiles of milk have shown to be a useful tool for verification of geographic origin [35][36][37]. The elements present in the natural environment go through a complex metabolic process inside the feed (eg., grass) and animal, and the variability of trace elements in milk may be affected by both intrinsic and extrinsic factors related to feed and the animal. A significant variation in climate, soil, and cattle management practices can be observed between the agroclimatic zones of Sri Lanka, though geographic distances between these zones are small (Figure 1). Al, Cr, Mn, Fe, Co, Ni, Cd, and Pb may be related to intrinsic factors of the animal, such as breed, stage of lactation, parity, and health status together with feeding and dietary adaptations based on the availability of feeding material. Li, Cu, Zn, As, Se, Sr, Ba, and Pb could be attributed to dietary supplements and feed additives used to feed animals [38][39][40][41][42]. Se, Cu, and Zn are present in mineral mixtures fed to intensively manage high yielding animals [40]. However, there are very limited data available on correlations between elemental content of milk and geographic origin for the South Asian region. The value of such information is vital to study the applicability of elemental content analysis for geographic origin verification for the region, as cattle management practices vary widely in the region even within the same agro-climatic zone.  Samples collected during the dry period were evaluated for elemental composition and stable isotope ratios. Previous studies have shown that there is a correlation between the elemental profile of animal products and geographic origin. Further, the elemental profiles of milk have shown to be a useful tool for verification of geographic origin [35][36][37]. The elements present in the natural environment go through a complex metabolic process inside the feed (eg., grass) and animal, and the variability of trace elements in milk may be affected by both intrinsic and extrinsic factors related to feed and the animal. A significant variation in climate, soil, and cattle management practices can be observed between the agroclimatic zones of Sri Lanka, though geographic distances between these zones are small (Figure 1). Al, Cr, Mn, Fe, Co, Ni, Cd, and Pb may be related to intrinsic factors of the animal, such as breed, stage of lactation, parity, and health status together with feeding and dietary adaptations based on the availability of feeding material. Li, Cu, Zn, As, Se, Sr, Ba, and Pb could be attributed to dietary supplements and feed additives used to feed animals [38][39][40][41][42]. Se, Cu, and Zn are present in mineral mixtures fed to intensively manage high yielding animals [40]. However, there are very limited data available on correlations between elemental content of milk and geographic origin for the South Asian region. The value of such information is vital to study the applicability of elemental  Geologically, underlying over 90% of Sri Lanka is Precambrian high-grade metamorphic rocks which can be subdivided into three main litho-tectonic units, namely Highland Complex, Vijayan Complex, and Wanni Complex. The remaining 10% is underlain by Miocene sedimentary rocks. Rocks underlain in the Highland Complex are mainly granulite facies suites, including charnockites, garnet-biotite gneisses, quartzite, marbles, and calc-silicate gneisses. The Vijayan region consists of gneisses and granitoid varieties, while the Wanni Complex is dominated by supracrustal rocks and meta-igneous suits, such as granodioritic gneisses. All the sampling locations of up country and mid country, and D1 and D4 locations of the dry zone belong to the Highland Complex of the metamorphic terrain while other sampling locations belong to Wanni Complex. The major and trace element content in the natural environment is mainly determined by their content in rocks, which significantly varies within a very short distance in a metamorphic terrain like Sri Lanka [43]. However, it is difficult to directly correlate the local geology with the resulted elemental contents in milk, as there are other factors influencing the elemental content in milk such as management practices. Li, Al, Cr, Mn, and Sr contents in samples collected in the dry period indicated significant differences in their mean contents among four agroclimatic zones (Table 6). Fe, Ni, Cu, Zn, As, Se, Ba, Pb, and Bi displayed moderate potential in differentiating agroclimatic zones, and Co and Cd composition were similar among the different farming zones. Sr and Mn concentrations together with isotope ratios have been used as markers for geographical traceability of raw milk in previous studies [36,37]. The first two principal components explained 58% of the total variance and clear clustering was observed in the score plot of PC1 and PC2 ( Figure 6). The LDA carried out using both trace element and stable isotope data with backward stepwise variable selection (i.e., milk casein δ 15 N, whole milk δ 15 N, milk casein δ 13 C, whole milk δ 18 O, whole milk δ 2 H, Li, Al, Cr, Mn, Fe, Co, Ba, and Sr) showed excellent discrimination in samples from the dry zone, coconut triangle, mid country, and up country (100%) even after cross-validation (Figure 7).      One limitation of the study is that the authors could not select farms in each of the agro-climatic zones with a good geographic distribution, due to unavailability of suitable farms throughout the zones. This uneven distribution of farms resulted from suitable agricultural land, water sources and other infrastructure availability, and environmental considerations. In this study, only large to medium scale farms were selected. Even though some of these farms were very close to the agroclimatic boundaries, they could be clearly categorized into respective agroclimatic zones in the LDA carried out with trace elements and stable isotopes data. This could have been possible due to the differences in elemental contents in the environment. However, we could not confirm whether the LDA model can differentiate milk coming from locations that were not sampled, especially in the eastern and southern regions of the dry zone where suitable farms are not available.

Conclusions
This study has demonstrated the ability of stable carbon, nitrogen, oxygen, and hydrogen isotope ratios, and multi-element analysis of cow milk to verify the agroclimatic origin of milk in a small tropical island with a diverse climate, environmental conditions, and cattle management practices. This study provides preliminary information on the applicability of stable isotope and elemental content analysis in verifying agroclimatic origin of milk, even within a small geographic area with diverse conditions. The isotope signatures significantly varied in accordance with the agroclimatic origin of the milk samples and were moderately influenced by the sampling period, detailed data can be found in the Supplementary Materials. The δ 18 O was the most correlated stable isotope with the agroclimatic zone and the LDA analysis of stable isotopes (whole milk δ 15 N, milk casein δ 15 N, milk whey δ 15 N, whole milk δ 13 C, milk fat δ 13 C, milk casein δ 13 C, whole milk δ 2 H, and whole milk δ 18 O) correctly classified more than 75% of the milk samples to their respective agroclimatic zones after cross-validation. The inclusion of selected elemental concentrations with stable isotope data (milk casein δ 15 N, whole milk δ 15 N, milk casein δ 13 C, whole milk δ 18 O, whole milk δ 2 H, Li, Al, Cr, Mn, Fe, Co, Ba, and Sr) of the milk samples enhanced the discriminating power of the LDA model up to 100%. A combination of selected elemental concentrations and stable isotopic compositions analyzed by a chemometric model is a promising tool to verify agroclimatic origin of cow milk besides high variations in climatic conditions and cattle management practices. Further, the data generated in the present study are expected to be used as baseline data to develop a model for non-targeted screening for potential chemical contaminants, such as dicyandiamide (DCD) and adulterants.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.