基于稳定同位素和矿物元素的牛奶产地溯源技术

王默涵，高 铭，李雅丽，柳 嘉，富宏坤，孙晓辰，段盛林,

（1.中国食品发酵工业研究院有限公司，北京 100015；2.功能主食创制与慢病营养干预北京市重点实验室，北京 100015；3.中国合格评定国家认可中心，北京 100062）

乳业是畜牧产业发展的重要组成，我国牛乳产业虽然起步较晚，但日益兴盛。牛奶是优质钙和蛋白质的良好来源，随着生活水平的提高和牛奶销售量的大幅增长，保证牛奶质量是迫切所需。自2008年三聚氰胺事件后，乳品质量安全备受关注，具有特定产地、草饲散养等声称的高端乳制品层出不穷，但也使其更容易被混淆甚至伪造，乳源控制因而开始成为生产商和消费者的关注点，牛奶产地属性也对确定产品价值具有重要意义[1]。

关于产地属性问题，世界各国通过出台法律法规或建立追溯制度来推进食品质量安全管理[2]。欧盟、美国和日本是较早开展食品追溯标准化工作的国家和地区[3]，欧盟早在1992年就已提出食品产地的概念，并在2000年要求食品的生产、加工和流通等环节建立溯源制度来明确产品来源[4]。我国在《农产品质量安全》以及《食品安全法》中也明确指出了食品农产品的追踪与溯源的重要性，对依法建立食品安全追溯体系、保证食品的可追溯性进行了明确说明[4- 5]。

稳定同位素和矿物元素是生物天然携带的一种专属标记，稳定同位素比值和矿物元素指纹信息与其生存环境和生长方式密切相关[5]，它已成功应用于葡萄酒[6]、谷物[7]、茶叶[8]、果蔬[9-10]、畜禽肉品[11-12]及水产品[13]中，是产地属性的有效保护工具，能为食品溯源提供一种科学独立、稳定流动的身份信息[14]，且由于受人为因素和检测环境的影响较小而具有准确度高、灵敏度高、无污染等特点[15]。

现有的牛奶产地溯源研究主要为仅测定牛奶中的稳定同位素比值或结合同位素和十余种甚至几十种矿物元素进行建模[5]，但若仅测定同位素，可能影响训练模型的鉴别准确度，若纳入的元素种类过多，成本和技术层面在推广时具有一定难度。综合考虑模型鉴别精度和成本，本研究采集黑龙江省、云南省、河北省、甘肃省、内蒙古自治区五地的牛奶样本，以牛奶中碳、氮、氢、氧四种稳定同位素比值（δ13C、δ15N、δ2H和δ18O）和钙、铁、锌三种重要矿物元素的质量浓度为指标对检测结果进行多元统计分析，基于样品产地建立数据模型，探究其对牛奶产地的判别能力，开发一种便于应用推广的牛奶产地溯源技术，助力奶源控制，提振国民对国产乳品的消费信心。

1 材料与方法

1.1 材料与仪器

9 mm × 5 mm锡舟、5 mm × 3.5 mm银囊 德国赛默飞世尔科技有限公司；钙、铁、锌单元素标准液（1 mg/mL） 国家标准物质认证中心；浓硝酸（优级纯） 江苏晶锐化学股份有限公司；正己烷 天津市大茂化学试剂厂；NaN3青岛雪洁助剂有限公司；实验用水均为GB/T 6682规定的一级水。

MARS 6微波消解系统 美国CEM公司；iCAP6000电感耦合等离子体发射光谱仪、Flash EA 2000 元素分析仪、高温裂解元素分析仪、 Delta V Advantage稳定同位素比值质谱仪 美国 Thermo Fisher公司；SL2002N电子天平 瑞士 Mettler-Toledo公司； DHG-9075A电热干燥箱 上海一恒科技有限公司；LGJ-18真空冷冻干燥机 北京松源华兴科技发展有限公司；GN-RO-40型超纯水仪 北京双峰众邦科技发展有限公司。

1.2 实验方法

1.2.1 样品采集 本研究采集的牛奶样本均为荷斯坦乳牛常乳，于2020年7月～8月分别取自我国黑龙江省、内蒙古自治区、河北省、甘肃省、云南省牧场，各省份样本均为同一牧场不同牛区或不同产奶时间的奶罐样，即不同胎次、泌乳期的多头泌乳牛的混合乳，以减少其他因素干扰，每个牧场采集5个独立样本。各样品均添加0.03% NaN3用于防腐，编号、登记后立即运回，直至检测前始终储藏于-20 ℃。采样信息统计如表1所示。

表1 采样信息统计Table 1 Sampling information statistics

1.2.2 稳定同位素分析

1.2.2.1 样品前处理 参考WIJENAYAKE等[16]的方法提取牛奶中的酪蛋白进行同位素分析。将牛奶样品冷冻干燥24 h后捣碎，准确称取12 g样品于50 mL离心管中，加入10 mL超纯水稀释，再加入5 mL正己烷，涡旋振荡5 min，静置30 min后弃去上层溶液，用0.1 mol/L稀盐酸将离心管下层溶液调节至pH4.6，以8000 r/min的条件离心10 min，离心后取下层沉淀，加入5 mL超纯水，重复用5 mL超纯水洗涤2～3次，水洗后的沉淀冷冻干燥24 h后待测。

1.2.2.2δ13C和δ15N值的测定 使用元素分析仪联用稳定同位素比值质谱仪测定样品的δ13C和δ15N值。元素分析仪左反应管用氧化铬填充，测样前升温至650 ℃；右反应管用铜丝填充，测样前升温至980 ℃，色谱柱温度为50 ℃；氦气流速为100 mL/min。分别测定参考物质IAEA-CH-6（δ13C=10.4‰±0.2‰）和USGS-34（δ15N=-1.8‰±0.2‰）的δ13C值和δ15N值，确定仪器处于稳定状态后，称取0.15～0.25 mg样品粉末于锡舟中检测。

1.2.2.3δ2H值和δ18O值的测定 使用高温裂解元素分析仪联用稳定同位素比值质谱仪测定样品的δ2H值和δ18O值。高温裂解元素分析仪反应管用玻璃碳填充，测定δ18O值前反应管阶段升温至1250 ℃，色谱柱温度为60 ℃，氦气流速为120 mL/min；测定δ2H值前反应管阶段升温至1400 ℃，色谱柱温度为85 ℃，氦气流速为120 mL/min。分别测定参考物质NBS-22（δ2H=-120.0‰±1.0‰）和USGS-34（δ18O=-27.9‰±0.6‰）的δ2H值和δ18O值，确定仪器处于稳定状态后，称取0.25～0.35 mg样品粉末于银纸杯中检测。

1.2.2.4 同位素比值的计算 采用国际通用的相对测量法表示物质的同位素比值，碳、氢、氧、氮元素的同位素特征均以国际标准物质为参照，结果以δ（‰）表示：

式中：X为重同位素的丰度，即13C或18O 或2H或15N；R为重同位素与轻同位素的丰度比，即13C/12C 或18O/16O或2H/1H或15N/14N；R样品为被测样品的同位素比值；R标准为国际标准物质的同位素比值。

1.2.3 矿物元素分析 分析方法及仪器操作条件均参考GB 5009.268-2016《食品中多元素的测定》[17]。准确称取各牛奶样品0.50 g于微波消解内罐中，加入5 mL浓硝酸，旋紧罐盖放置于通风橱中酸解过夜，次日放入微波消解仪中消解，冷却后取出，缓慢打开罐盖排气并用少量水冲洗内盖，将消解罐放在控温电热板上于100 ℃加热30 min赶酸排气，取出恢复至室温后用超纯水定容至50 mL，混匀备用，同时做空白试验。采用电感耦合等离子体发射光谱仪，通过外标法定量测定样品中的钙、铁、锌三种矿物元素的质量浓度。

1.3 数据处理

本研究中的所有指标均进行了至少3次的平行测定，结果表示为平均值±标准差，数据处理和作图采用IBM SPSS Statistics 25。各组间通过单因素方差分析和邓肯多重比较确定显著性差异（P＜0.05）；所有数据标准化后，通过聚类分析（系统聚类，Ward法，欧氏平方距离）、主成分分析和判别分析（Fisher 线性判别分析，“留一法”交叉验证）区分每个牛奶样品的地理来源。

2 结果与分析

2.1 稳定同位素分析

稳定同位素是理想的产地溯源参数，因为其在动物组织中的分布水平与动物的饮食摄入和地理环境密切相关[18]。其中，δ13C主要与奶牛的食物成分有关，δ15N主要反映饲料种植及奶牛生存环境的土壤情况，δ2H和δ18O则与当地的水源、空气等环境因素及气候条件密切相关[5]。结果显示，来自黑龙江省、内蒙古自治区、河北省、甘肃省、云南省等5个产地的牛奶样品，其酪蛋白的δ13C、δ15N、δ2H和δ18O值地域差异性显著，结果如表2所示。

表2 五地牛奶样品稳定同位素比值测定结果（‰）Table 2 The stable isotope ratio in milk samples from five region(‰)

由于C3植物和C4植物光合作用的途径不同，可能导致碳同位素比值存在差异[5]。牧草、小麦和大豆等C3植物的δ13C值范围通常在-24‰至-34‰之间，而玉米和高粱等C4植物的δ13C范围则通常介于-10‰至-23‰之间[19]。经奶牛摄入和消化后，这种特征可被浓缩至奶牛组织中，因此乳汁通常带有饲料类型或比例特征[20]。MOLKENTIN等[21]的研究表明，无玉米喂养的奶牛乳汁δ13C的极限值为-23.5‰，而无任何C4植物喂养的奶牛乳汁中的δ13C值更负，草饲奶牛乳汁的δ13C值显著低于精饲奶牛（如以玉米、高粱、青贮等为原料经精加工而成的饲料）。本研究中的δ13C值均介于-18‰至-20‰之间，表明5个采样牧场均表现为缺乏新鲜牧草、高谷物水平的精饲料喂养模式，显著高于同类研究中爱尔兰、新西兰、澳大利亚等具有丰富牧草资源国家的酪蛋白δ13C值，这些国家草饲奶牛乳汁的δ13C值通常低于-24‰[22]。具体δ13C值表现为黑龙江省和河北省＞甘肃省＞云南省和内蒙古自治区（P＜0.05），呈现东部、低海拔地区（黑龙江和河北）显著高于西部、高海拔地区（内蒙古、甘肃、云南）的趋势（P＜0.05）。由于采样牧场均为以玉米、棉籽、青贮等为主要原料的精饲喂养，因此在饲料类型大致相同的情况下的差异可能与不同地域饲料作物的光合作用有关，但尚无δ13C值与海拔之间相关趋势的直接研究。

氮同位素比值可用于区分饲料作物的施肥方式和固氮作物种类在奶牛饮食中的比例[19]。合成肥料的δ15N值通常介于-4‰至4‰之间，而有机氮肥（包括粪便和堆肥）的δ15N值范围可宽至0～30‰[23]，部分固氮作物（如豆科植物）的δ15N值则可能接近0‰[24]。本研究中牛奶酪蛋白具有相似的δ15N 值，均在2‰至4‰之间，且内蒙古自治区和河北省＞黑龙江省和云南省＞甘肃省（P＜0.05），均显著低于爱尔兰、新西兰、澳大利亚等国家[22]，这可能是由于奶牛饮食中合成氮肥的使用比例较大而有机氮肥比例较小，且饲料中豆粕等固氮作物的比例较高，差异主要反映在耕作及喂养模式上，与经纬度、海拔的相关性较小。

δ2H值和δ18O值的结果趋势表现出高度的一致性：甘肃省＞河北省和内蒙古自治区＞黑龙江省和云南省，该趋势也与δ15N值表现出一定的一致性，这可能与气候条件、距海远近、海拔、空气、水源等养殖环境有关[18]，但与爱尔兰、新西兰、澳大利亚等国家及我国的其他同类数据有较大差距[1,22]，这可能与不同的采样时间和采样地点有关。

可见，稳定同位素比值是鉴别养殖模式（草饲、精饲等）的有效参数，也为地理溯源提供可靠依据。

2.2 矿物元素分析

矿物元素常具有丰度效应，不同地域间土壤、水源及当地奶牛饲料中的矿物元素含量可能存在差异，故地域可能影响乳中矿物元素含量[25]。ZHOU等[26]评估了我国家畜乳中元素含量在地域间的差异，结果表明，不同地域家畜乳中元素含量不同，饮水与饲料是可能的影响元素。本研究通过单因素方差分析和邓肯多重比较（见表3）显示出牛奶中钙和锌的质量浓度在5个产地的分布差异显著（P＜0.05），且分布规律一致，均在甘肃省最高、内蒙古自治区最低；而牛奶样品中铁元素的含量较低，在5个产地的分布差异不显著（P＞0.05）。鉴于采样时已排除个体差异的影响，此结果反映了不同养殖环境下差异化的营养状态，表明其对区分五地区的样品具有一定的可行性。

表3 五地牛奶样品矿物元素质量浓度测定结果Table 3 The concentration of mineral elements in milk samples from five region

2.3 聚类分析

聚类分析是归类划分数据中相似性数据的统计方法[27]。使用以上四项同位素比值及三项矿物元素质量浓度的标准化得分进行聚类分析，产生的聚类树状图如图1所示。将聚类距离为5时作为分割线，五产地的牛奶样品可被分为4类。其中甘肃省的样品与其他省距离较远，可被完全分离；河北省的样品虽同属一类，聚类效果较好，但其他产地的个别样品也分类至这类中；云南省和黑龙江省的大部分样品被归为一类，表明两产地的聚类距离相近，无法完全分离；内蒙古自治区的样品聚类效果不好，分散在各类别之中。可见，聚类分析为产地鉴别提供了可能，但仍存在交叉现象。

图1 不同产地牛奶的系统聚类树状图Fig.1 Dendrogram of cluster analysis from different regions

2.4 主成分分析

利用稳定同位素比值和矿物元素组合参数的主成分分析（PCA）可对不同产地的牛奶进行分类[5]。将所有数据标准化后，KMO检验统计量大于0.6，Bartlett球形检验显著性P＜0.05，表明本实验模型适合做主成分提取[28]。PCA共提取了7个主成分，总方差解释如表4所示，其中初始特征值大于1的主成分有3个，累积方差百分比为81.95%，可代表大部分数据特征。以前三个主成分为代表的模型为：

表4 总方差解释Table 4 Total variance explained

利用此三个主成分的标准化得分值做散点图，结果如图2所示，甘肃省的样品与其他省区的差异明显，可进行区分，其他地区样品仍有交叉性。

图2 不同产地牛奶的主成分得分散点图Fig.2 Scatter plot of principal component score from different regions

2.5 判别分析

判别分析是一种有监督的数据分类方式[29]。其可利用各类样品之间的相关性建立判别模型而对未知样品进行分类，还可利用交叉验证把除某一个样品之外的其余样品作为训练集，建立判别函数对该样品进行判别，从而极大地检验了判别模型的稳定性[30]。基于牛奶产地进行元素差异分析、聚类分析和主成分分析后可以看出，本研究的组合参数对牛奶产地溯源具有一定可行性。为进一步了解各参数对牛奶产地的判别效果，利用Fisher函数建立判别模型[31]。

本次判别分析建立了4个典型判别式函数，其中函数1的方差贡献率为96.1%，函数2的方差贡献率为2.8%，利用这2个判别函数分别作为横纵坐标建立坐标系，得到5个组的散点图，如图3所示，可以看出各产地样品相对集中，根据产地分为五组后各组质心也可较好地分离。

图3 不同产地牛奶的判别分析散点图Fig.3 Scatter plot of discriminant analysis from different regions

采用交叉验证检验模型的可靠性，结果如表5所示。原始分类结果显示：7种元素分析的原始整体判别率为100%，五产地的样品均可被100%正确识别；“留一法”交叉验证的分类结果中，有96%的牛奶样本被正确分类，其中黑龙江省的样本存在20%的误判，其余省区的样本均可被100%正确判别。说明本次判别分析建立的判别函数对来自5个地区的牛奶样品具有很好的分类效果。所建立的判别函数模型为：

表5 不同产地牛奶的判别分析分类结果（%）Table 5 Classification result of discriminant analysis from different regions (%)

同类的牛奶产地溯源研究中，ZHAO等[1]针对δ13C、δ15N、δ2H、δ18O四种稳定性同位素比值建立数据模型，得到河北省、宁夏回族自治区、陕西省、内蒙古四个区域的交叉验证准确率71.55%，河北省、天津市、江苏省、内蒙古自治区四个区域的交叉验证准确率92.11%；侯建波等[5]除稳定性同位素外，纳入了40种矿物元素，选出As、Se、Rb、Sr、Mo、Cs、Ba、Tl、δ13C、δ15N和δ2H组合指标，对德国、西班牙、意大利、澳大利亚、新西兰、中国和奥地利七个区域的牛奶整体判别率为 92.4%。本研究的分析方法与仅纳入稳定同位素的数据模型相比判别准确率提高，同时又控制了成本投入，无需检测过多矿物元素即可实现国内多省市区牛奶样本的判别，是一种更适合应用推广的方法。

3 结论

本研究分别测定黑龙江省、云南省、河北省、甘肃省、内蒙古自治区五个产地的牛奶样品的稳定同位素及矿物元素指纹信息，建立的数据模型可对奶源地喂养模式和牛奶产地实现有效判别。结果表明：δ13C、δ15N、δ2H、δ18O在产地间的差异显著（P＜ 0.05），这可能与地理属性和养殖模式有关，可区分如草饲、精饲等的喂养模式，也可为地理溯源提供可靠依据；钙、锌在产地间的差异显著（P＜0.05），而铁在产地间的差异不显著（P＞0.05），表现出不同养殖环境差异化的营养状态，是地理溯源的有力指标。基于产地的系统聚类和主成分分析对牛奶产地溯源具有一定可行性，但仍存在部分交叉；Fisher线性判别分析则可有效判别牛奶产地，原始整体判别率为100%，交叉验证的正确判别率为96%，其中对内蒙古自治区的正确判别率为80%，其余均为100%。本研究结果证明利用稳定同位素碳、氮、氢、氧与钙、铁、锌相结合的溯源技术可对我国黑龙江省、云南省、河北省、甘肃省、内蒙古自治区五个产地的牛奶进行有效的产地判别，是产地溯源的有利技术。