李 梁 李春霖 王贞红 刘振东 聂 晶 张永志 袁玉伟,*

(1 西藏农牧学院食品科学学院/西藏特色农牧资源研发省部共建协同创新中心，西藏 林芝 860000；2 浙江省农业科学院农产品质量安全与营养研究所/农业农村部农产品信息溯源重点实验室，浙江 杭州 310021；3 西藏农牧学院资源与环境学院，西藏 林芝 860000)

我国西藏人民饮茶历史悠久，茶叶对于藏区人民的日常饮食起着重要作用[1]。过去西藏茶叶产量较低，主要依赖从四川等地运输。近年来，西藏高度重视茶产业发展，把茶产业作为促进经济发展和助推脱贫攻坚的重要产业之一[2]。西藏茶叶主产于林芝市，2019年林芝市茶园种植面积为4.5万亩，2020年春季干毛茶产量约1.55万公斤，带动1 400余人人均增收3 445元，茶产业已成为西藏人民致富的“金叶子”[3]。西藏主要生产红茶、绿茶和黑茶类产品，由于高海拔和独特的地理环境，西藏茶叶品质较优，香气丰富，滋味醇厚，具有高原特色。然而，西藏茶叶的产地特征尚不清楚，市场中可能出现的假冒标识问题，对西藏茶叶的贸易公平和品牌保护产生不利影响。因此，挖掘西藏茶叶的产地特征，开展溯源机理研究对西藏茶叶质量监管和产业发展具有重要意义。

植物稳定同位素值(δ13C、δ15N、δ2H和δ18O)与栽培地的气候环境和栽培方式等因素密切相关[4-5]。已有研究表明，茶叶中稳定同位素存在一定的地域差异，具有特征指示性。Liu等[6]研究发现，浙江和山东茶叶中稳定同位素比率有明显差异，结合多元素和化学计量学方法，能够对产地进行有效判别。Jin等[7]分析了福建安溪和武夷山茶叶中稳定同位素比率，发现两地样品差异显著。植物稳定同位素在不同地域表现的差异主要是由环境和气候因子的变化导致的，这些因子通过调控植物生理过程从而引起同位素特征发生改变[8]。目前已有相关研究探索农产品中稳定同位素与环境因子的关联性。如Camin等[9]测定了意大利葡萄酒样本的稳定同位素比率和环境气候因子，发现δ18O值与气候和地理位置信息相关性最强；Akamatsu等[10]统计分析了日本水稻中δ13C和δ18O值与环境因子之间的关系，结果表明灌浆期的日最低气温是影响大米同位素的主要因素；Wang等[11]研究了气候、地形和土壤等生长环境对香蕉同位素组成的影响，发现降水和气温对δ18O值影响较大。已有研究发现，茶叶中稳定同位素受到环境温度和光照强度的影响[12]，δ2H值与气温、光照和降水量呈现显著相关性[13]。上述研究对探索茶叶同位素特征与环境因子关系具有一定的借鉴意义。

西藏茶叶作为我国特色农产品，其产地同位素特征和环境影响因素尚不明确。本研究通过分析西藏不同海拔地区茶叶的同位素特征，研究西藏茶叶与我国四大主产区(江北茶区、江南茶区、华南茶区、西南茶区)茶叶的同位素差异，挖掘环境因子对茶叶中同位素的影响，旨在为西藏茶叶品牌保护提供技术支撑和科学依据。

1 材料与方法

1.1 材料与试剂

本试验材料为93份茶叶样品，收集于2017年4月，来自全国12个省，包括27个地市及自治州，覆盖我国江北(山东、陕西)、江南(安徽、浙江)、华南(福建、广东、海南)和西南(贵州、四川、云南、重庆)四大茶区，其中江北茶区样品9个、江南茶区样品21个、华南茶区样品15个、西南茶区样品30个。西藏自治区样品18个，取自林芝市不同海拔的3个取样点，其中低海拔(750 m)样品5个、中海拔(1 250 m)样品8个、高海拔(2 120 m)样品5个。所有样品均于2017年春季按照绿茶工艺加工制成。样品分布情况见图1。

图1 样品分布及气象信息收集图Fig.1 Collection map of tea samples and climate information

稳定同位素标准物质：IAEA-CH-6(蔗糖，δ13CV-PDB= -10.45‰±0.03‰)、IAEA-N-2(硫酸铵，δ15Nair=20.30‰±0.20‰)，购于国际原子能机构(international atomic energy agency, IAEA，奥地利)；B2155(酪蛋白，δ13CV-PDB=-26.98‰±0.13‰，δ15Nair=5.94‰±0.08‰)，购于英国Elemental Microanalysis公司；USGS64(甘氨酸，δ13CV-PDB=-40.81‰±0.04‰)、USGS40(δ15Nair=-4.52‰±0.06‰)、USGS54(加拿大松，δ2HV-SMOW=-150.40‰±1.10‰，δ18OV-SMOW=17.79‰±0.15‰)、USGS56(南非红象牙木，δ2HV-SMOW=-44.0‰±1.8‰，δ18OV-SMOW=27.23‰±0.03‰)，购于美国地质勘探局Reston同位素实验室。

1.2 仪器与设备

XP6型天平，瑞士Mettler-Toledo公司；Vario PYRO cube、Vario Isotope cube型元素分析仪，Isoprime 100型、Biovision稳定同位素比率质谱仪，德国Elementar公司；HR2864粉碎机，中国飞利浦电子公司。

1.3 试验方法

1.3.1 样品预处理 茶叶样品烘干至含水量不变，经粉碎机粉碎处理，过80目网筛，避光常温保存，待测。

1.3.2 碳、氮稳定同位素检测 称取茶叶样品3.0～4.0 mg，放入锡箔舟中包样，按顺序放入元素分析仪固体样品自动进样盘。样品中的碳、氮元素经燃烧还原后分别转化为纯净的CO2和N2气体，进入同位素比率质谱仪进行检测。元素分析仪的燃烧管主要填料为WO3,温度保持1 150℃，还原管填料为Cu，温度保持850℃，载气为高纯氦气(He, 99.999%)，流量为250 mL·min-1，参考气体为CO2和N2。元素分析仪和热导检测器(thermal conductivity detector, TCD)检测流速分别为230 mL·min-1和40～50 mL·min-1。氧气流量为40 mL·min-1，注氧时间70 s。C模式质谱参数：加速电压为3 964 V，捕集电流为200 μA，磁场强度为4 000 mA；N模式质谱参数：加速电压为4 162 V，捕集电流为400 μA，磁场强度为3 000 mA。样品检测前进行气密性检测、离子源真空度和稳定性测试，CO2参考气压力为12 psi，N2参考气压力为5 psi。检测方法参照文献[6]和[13]。

1.3.3 氢、氧稳定同位素检测 称取茶叶样品0.5 mg，放入银舟中包样，按顺序放入元素分析仪固体样品自动进样盘。样品中的氢、氧元素经高温裂解后分别转化为H2和CO气体，进入同位素比率质谱仪进行检测。元素分析仪裂解炉温度保持1 450℃，氦气吹扫流量为150 mL·min-1，参考气体为CO和H2。H模式质谱参数：加速电压为4 497 V，捕集电流为600 μA，磁场强度为890 mA；O模式质谱参数：加速电压为 4 191 V， 捕集电流为200 μA，磁场强度为3 000 mA。样品检测前进行气密性检测、离子源真空度和稳定性测试，H2参考气压力为20 psi，CO参考气压力为25 psi。检测方法参照文献[6]和[13]。

1.3.4 稳定同位素比率计算

(1)

式中，R样品为检测样品中重同位素与轻同位素丰度之比，即13C/12C、15N/14N、2H/1H、18O/16O；R标准为国际标准品中重同位素与轻同位素丰度比[14]。测试数据均采用两点法进行校正。

1.4 气象信息采集

在试验样品采集点中，选取16个采样点(图1)进行气象信息采集，涵盖了西藏和我国茶叶主产区，收集当地环境气象信息，包括平均相对湿度、最小相对湿度、日降雨量、日照时间、平均气温、最高气温、最低气温，每个环境因子均包括采样当月(2017年4月)和采样前3个月(2017年1月至4月)的2组平均数据，共计14项环境数据。

1.5 数据处理

采用正态检验和方差齐性检验判断数据是否符合正态分布和方差齐性，采用单因素方差分析中的Duncan法对不同海拔样品的同位素比率进行差异显著性分析，采用Tamhane’s T2方法对不同产区样品的同位素比率进行差异显著性分析，在IBM SPSS Statistics 21软件中进行。

通过偏最小二乘法(partial least squares，PLS)分别建立4种同位素比率与环境因子之间的关联性模型，以样品14项环境因子数据为自变量，稳定同位素比率为因变量，建立定量预测模型。模型共包含16个样本，随机选取其中14个样本为训练集，2个样本为预测集。模型结果根据所有样品的决定系数(R2)、均方根误差(root mean square error，RMSE)和相对预测偏差(relative prediction deviation，RPD)表示，计算方法见公式(2)～(4)[15-16]。结合连续变量投影法进行变量重要性分析(variable importance for projection，VIP)，并计算自变量的标准化相关系数。PLS和VIP在XLSTAT 2019软件中进行。

(2)

(3)

(4)

式中，N为样品数量,ypred为同位素预测值,yref为同位素检测值,yref为同位素检测值的平均数，SD为同位素检测值的标准差。

2 结果与分析

2.1 西藏不同海拔地区茶叶稳定同位素特征

本研究共采集西藏茶叶样品18份，稳定同位素比率分析结果见表1。由于西藏海拔差异较大，样品取自3个不同海拔的取样点，不同海拔的样品稳定同位素比率差异整体显著。δ13C随海拔变化最为明显，茶叶样品中δ13C值随海拔高度的升高而减小。低海拔地区(750 m)样品δ13C值介于-28.8‰ ～ -27.8‰ 之间，平均值为-28.3‰；中海拔地区(1 250 m) 样品δ13C值介于-30.1‰ ～ -28.9‰之间，平均值为-29.6‰；高海拔地区(2 120 m)样品δ13C值介于-31.9‰ ～ -30.9‰之间，平均值为-31.2‰。δ15N 值在低海拔地区样品中最高，平均值为3.0‰，显著高于中海拔(0.9‰)和高海拔地区(0.0‰)的样品。δ2H和δ18O随海拔变化趋势一致，均为低海拔样品最高，分别为-48.8‰和27.3‰，中海拔样品最低，分别为-100.9‰和23.4‰。

2.2 西藏茶叶与四大主产区茶叶稳定同位素分布差异

西藏由于其独特的地理位置，与我国茶叶主产区在环境气候条件上差异明显，稳定同位素比率也呈现一定的特征差异。将西藏样品分别与江北、江南、华南、西南(除西藏)四大主产区代表性样品进行差异显著性分析，结果见图2。西藏样品中δ13C值特征最为明显，介于-31.9‰ ～ -27.8‰之间，平均值为-29.7‰，显著低于其他主产区样品。δ15N值介于-2.2‰ ～ 5.0‰之间，平均值为1.2‰，与其他产区样品无显著差异。δ2H值介于-111.5‰ ～ -40.5‰ 之间，平均值为-83.3‰，显著低于西南、江南和江北产区样品，与华南产区样品无显著差异。δ18O值介于22.2‰ ～ 27.6‰之间，平均值为25.3‰，显著高于华南产区样品，与其他产区样品均无显著差异。

注：**代表P<0.01，*代表P<0.05，ns代表P>0.05。Note: **represent P<0.01, * represent P<0.05, ns represent P>0.05.图2 西藏与四大主产区茶叶中稳定同位素箱线图Fig.2 Boxplot of stable isotopes in tea samples from Tibet and four main producing area

2.3 茶叶稳定同位素与环境因子关系模型

为探索西藏茶叶同位素特征形成的气候环境影响，本研究选取西藏和我国茶叶主产区范围内16个取样点进行环境因子数据采集和茶叶同位素比率检测，分别将δ13C、δ15N、δ2H、δ18O与14项环境因子数据进行PLS分析，定量预测模型结果见表2。其中R2和RPD越高，RMSE越低，模型效果越优[15-16]。4个模型的R2均高于0.55，RMSE均低于5.09，RPD介于1.51～3.08之间。δ2H模型预测效果最优，R2和RPD均最高，分别为0.90和3.08，表明δ2H与14项环境因子关联性更强。由于δ2H绝对值较其他同位素比率偏大，而RMSE受自变量绝对值影响，导致δ2H模型RMSE最高，为5.09。其次为δ13C模型，R2为0.84，RPD为2.63，RMSE最低，为0.48。δ18O模型R2为0.75，RMSE为0.64，RPD为2.13。δ15N模型预测效果最差，R2和RPD均最低，分别为0.55和1.51，RMSE为1.20，表明δ15N与14项气候环境因子关联性较弱。

表1 西藏不同海拔茶叶样品稳定同位素分布表Table 1 Distribution of stable isotopes in tea samples from Tibet with different altitudes /‰

对模型进一步进行变量重要性分析和标准化相关系数分析，VIP值大于1的环境因子对同位素比率影响较大[17]，结果见图3。δ13C主要受日照和气温影响，重要性最高的3个因子为4月(采样当月)日照时数、平均气温和最低气温，3个变量因子的标准化相关系数均为正值，表明这3个因子与δ13C成正相关。在δ15N预测模型中，4月日照时数、1-4月(采样前3个月)日照时数和1-4月平均相对湿度影响较大。δ2H模型显示，日照、气温和相对湿度均会对δ2H影响较大，VIP大于1的因子以相对湿度为主，VIP最高的3个因子分别为4月日照时数、1-4月最高气温和1-4月最小相对湿度。δ18O则主要受气温影响，VIP大于1的因子主要包括1-4月和4月的最高气温、平均气温和最低气温，6个因子的标准化相关系数均为负值，表明这6个因子与δ18O成负相关。

表2 偏最小二乘模型结果Table 2 Results of partial least squares models

图3 基于偏最小二乘的变量重要性分析与标准化相关系数Fig.3 Variable importance for projection and standardized coefficients based on partial least squares

3 讨论

3.1 环境因子对碳、氮稳定同位素的影响

本研究结果发现，相较于其他同位素比率，西藏茶叶中δ13C值随海拔变化规律最明显，茶叶中δ13C值随海拔升高而降低。经PLS-VIP分析表明，δ13C值主要受采样当月日照和气温影响，且δ13C值与日照时数和平均气温成正相关。而气温随海拔的升高而降低，因此茶叶中δ13C值随海拔的升高而降低。在西藏与四大主产区茶叶同位素比率比较研究中发现，西藏茶叶δ13C值显著低于四大主产区茶叶。而西藏林芝市4月平均气温仅10℃，较其他主产区偏低。该结果与PLS分析中δ13C值与气温成正相关的结果一致。已有研究分析了庐山不同海拔茶树光合响应的差异，结果表明茶树在低海拔处的光合作用和光系统Ⅱ活性显著低于高海拔处[18]。此外，气温和CO2浓度升高能够明显提高茶树叶片的光合效率[19]，而光合作用与植物13C稳定同位素分馏密切相关[20-21]，进一步解释了茶叶δ13C值随海拔升高而降低的原因。

经单因素方差分析，西藏不同海拔样品中δ15N值差异不如δ13C值明显，而西藏与不同产区茶叶样品中δ15N值差异不显著，并且在PLS分析中δ15N值与环境气象因子关联性较弱。这说明相较其他同位素比率，δ15N值受气候影响较小。已有研究表明，植物中δ15N值与栽培和土壤关系更为密切，主要受到栽培方式[22]、土壤类型[23]、根系情况[24]等多方面影响。本研究有待进一步收集西藏茶叶种植环境样品，探究土壤背景、栽培方式以及肥料种类对茶叶同位素的影响。

3.2 环境因子对氢、氧稳定同位素的影响

不同海拔茶叶样品中δ2H和δ18O值分布规律较为一致，但西藏茶叶与主产区茶叶同位素比较中，δ2H和δ18O值分布规律有较大差异，推测与环境影响因素有关。已有研究表明，δ2H值与降雨、气孔开合等密切相关，具有显著的纬度效应、海拔效应和内陆效应[25]。本研究PLS结果表明，δ2H与环境气象因子关联性最强，利用14项环境因子能对δ2H进行定量预测。VIP结果显示，光照、气温、湿度多种环境因素均会影响茶叶中δ2H值，VIP大于1的变量中以相对湿度为主。该结果与已有研究结果较为一致，如Deng等[13]利用同位素和多元素技术对我国绿茶进行产地溯源，探究了不同产地绿茶的同位素特征及其与气候环境的相关性，发现茶叶中δ2H值与气温、光照和降水量相关性显著；Liu等[26]研究了不同产地绿茶的稳定同位素比率，发现湖北的绿茶样品中δ2H和δ18O值显著高于贵州和广东省样品，推测与环境相对湿度和蒸腾速率等因素有关。本研究PLS-VIP分析结果表明，气温是茶叶中δ18O值的主要影响因素，影响δ18O值的气温包括了采样当月和采样前3个月的平均气温、最低气温和最高气温。此外，除了气温对δ18O值的影响外，也已有研究表明植物纤维素中δ18O值与羰基和水之间的平衡同位素效应等因素相关[27]。

4 结论

本研究结果表明，西藏茶叶中稳定同位素比率具有一定的海拔规律和区域特征，其中δ13C值特征最为明显，随海拔升高而降低，显著低于其他主产区样品。茶叶稳定同位素比率均受到环境因子影响，其中δ2H值与气象环境关联性最强，主要与相对湿度等因素相关。在同位素特征分析的基础上，后续研究可结合代谢组学、光谱学、质谱学等技术，进一步挖掘西藏茶叶品质特征，助力西藏茶叶品牌建设与推广，推动西藏茶产业快速发展。