用蒙特卡洛方法评定转基因玉米59122的测量不确定度

宋 君， 刘月悦， 王 东， 张富丽， 李 洁

(1.四川省农业科学院分析测试中心，四川成都 610066； 2.四川省农业科学院农村经济与农业信息研究所，四川成都 610066

转基因玉米59122品系是一种既能杀虫(鞘翅目)又能抗除草剂的“双抗”型转基因作物。该转基因品系由美国陶氏农业科学公司(DOW AgroSciences LLC)和先锋国际公司(Pioneer Hi-Bred International Inc.)开发，2001—2003年通过食用和环境安全评价。2004年转基因玉米59122被投放到美国和墨西哥的食品和饲料市场。2005年美国和墨西哥开始商业化种植该品系。2010年德国开始对59122玉米品系进行环境释放。截至目前，59122玉米品系商业化种植时间已有12年[1]。我国于2006年批准进口转基因59122玉米用作食品和饲料的加工原材料，但不允许商业化种植。截至目前，我国已经连续进口该转基因玉米品系11年[2]。转基因玉米59122品系包含2个杀虫基因Cry34Ab1、Cry35Ab1以及1个耐除草剂基因pat，共3个目的基因。抗虫基因分别由启动子Ubi Zm1、P-peroxidase和终止子PinII调控，耐除草剂基因由启动子CaMV35S和终止子NOS调控。为了区分转基因产品与非转基因产品，贯彻我国《农业转基因生物安全管理条例》《农业转基因生物标识管理办法》和《转基因食品卫生管理办法》等法律、法规，宋君等针对转基因玉米59122品系建立了转基因玉米59122品系的实时荧光定量PCR检测方法，并验证了该法的特异性和灵敏度[3]。

近年来，随着我国经济和技术的快速发展，“中国制造”的各领域产品出口到世界各国，尤其是食品和农产品贸易遍及全球。很多进口国的法律和法规对含有转基因成分的食品和农产品有超过限量标识的规定，因此，精准定量分析产品中转基因成分含量对我国进出口贸易有重要意义。测量不确定度评定是继误差理论后发展起来的对测量质量进行科学评估的重要手段。目前，测量不确定度已经被世界各国广泛认可和接受。产品进口国往往要求出口方提供产品的质量和安全检测数据及其不确定度以评估检测数据的可靠性。不确定度评定在我国的大量应用始于21世纪初，但应用范围主要局限于机械、工程、航空、物理和化学等领域[4-8]，在生物领域的应用相对较少[9]。由于转基因生物安全学是近年发展起来的新学科，因此，测量不确定度评定在转基因生物及产品安全检测中的应用更是非常有限[10]，尚未得到推广和普及。

测量不确定度评定的蒙特卡洛法(Monte Carlo Method，MCM)是2008年由国际计量学组织(JCGM)针对《不确定度评定表达指南》(Guide to the Expression of Uncertainty in Measurement，GUM)存在线性模型近似假设、复杂模型灵敏度的偏导数难求、输入量概率密度函数对称和输出量概率密度函数正态分布的假设等缺陷[11]，提出的基于数字计算仿真模拟实施“概率分布传播”的测量不确定度评定的另一种方法。与GUM法[12]相比，可能是因为须要涉及到计算机编程等问题，测量不确定度的蒙特卡洛评定方法在国际上应用相对较少[13-15]。截至目前，国际上尚未有采用MCM评定转基因成分测量不确定度的研究报道。为了明确影响转基因成分定量检测不确定度的因素及其贡献大小，提高转基因成分检测质量，维护我国农产品贸易利益，本研究首次采用MCM法评定混合样品中转基因玉米59122的测量不确定度，以期为转基因生物及产品的测量不确定度评定提供新方法探索和实践指导。

1 测量过程

1.1 标准物质及试样

采用含量为10%的转基因玉米59122有证标准物质(certified reference material，CRM)作为仪器的校准子(calibrator)。用含有0.5%转基因玉米59122的混合样品作为测试样品。

1.2 主要仪器设备

用美国热电公司生产的型号为Nanodrop-1000的分光光度计测量标准物质和试样的核酸浓度。转基因玉米59122片段(载体与受体植物基因组相结合的侧翼片段)的实时扩增和59122片段的拷贝数测量在美国ABI公司生产的7500型实时PCR系统上进行(全部试验工作于2014年在四川省农业科学院分析测试中心转基因生物安全研究室完成)。

1.3 关键化学试剂

本研究涉及到的主要试剂包括核酸提取和实时PCR化学试剂，购自天根生物技术有限公司(北京)。

1.4 引物、探针合成

本研究使用的引物/探针和试验方法按照文献[3]报道的引物/探针序列和方法开展试验。引物/探针由上海生工合成和标记。

1.5 校准曲线拟合和试样的测量

将分离得到的10%含量的转基因玉米59122(CRM)的DNA溶液分别稀释成100、20、4、0.8、1.6×10-4μg/μL 5个浓度。根据玉米基因组大小2 504 Mbp[16]计算，0.1 μg/μL DNA溶液相当于1 μL DNA溶液中含36 400 copies DNA分子。用3倍体积的上述5个浓度点的DNA(Q)和经PCR反应得到的仪器响应值(threshold of cycles，Ct)作为校准子拟合校准曲线。实时PCR反应体系含探针型主要混和液(dNTPs、Mg2+、聚合酶等成分)12.5 μL，浓度为10 μmol/L的上/下游引物各1 μL，浓度为10 μmol/L的探针0.5 μL，然后添加灭菌的二次蒸馏水至25 μL。5个DNA浓度点分别做3次平行试验；测试样品做16次重复测量。实时PCR程序为95 ℃，10 min； 95 ℃，15 s；60 ℃，1 min(40 cycles)。

2 转基因玉米59122分子片段的数学测量模型

59122片段绝对含量的对数值(lgQ)与仪器响应值(Ct)之间存在函数关系，

Ct=m×lgQ+k。

(1)

式中：Q为内源基因zSSIIb或59122片段的绝对含量，μg或copies；Ct为PCR仪对Q的响应值；k为校准曲线的截距；m为校准曲线的斜率。

将试样的zSSIIb和59122片段的Ct值分别代入各自的校准曲线，计算得到试样zSSIIb和59122片段的绝对含量Q，然后将Q代入式(2)，获得59122片段的相对含量C。

C=Qex/Qed×100%。

(2)

式中：Qex为zSSIIb基因的绝对含量；Qed为59122片段的绝对含量；C为59122片段的相对含量，%。

3 MCM的实施

3.1 m和k的联合分布及计算公式

按照GUIDE 98-3/SUPP.1概率分布传递原理[17]，使用MATLAB软件实施MCM评定转基因玉米59122片段的测量不确定度。由于斜率m和截距k线性相关，所以m和k属于联合分布(双元变量正态分布)，用矢量E及其协方差矩阵V来表示。

(3)

(4)

式中：u(k，m)为k和m的协方差；u(m)为m的不确定度；u(k)为k的不确定度。

u(m)和u(k)分别由式(5)、式(6)计算：

(5)

(6)

而s(Ct)由式(7)计算：

(7)

式中：n为校准测量次数；

Cti为PCR仪对第i个校准子的响应值。

协方差u(k，m)由式(8)计算：

u(k,m)=r(k,m)u(k)u(m)。

(8)

相关系数

(9)

3.2 MCM的输入

5个校准子的对数值(Qi和Qi#)与其对应的响应值Cti和Cti#(表1、表2)以均匀分布R(a，b)(a为均匀分布的下限，b为上限)作为输入量，拟合内源基因zSSIIb和59122片段的校准曲线。在zSSIIb基因扩增过程中，0.1 μg/μL浓度点的响应值Ct1为22.44)，标准差(SD)为0.04(表1)，因此，循环阈值Ct1的分布范围为Ct1～U[22.40，22.48]。类似地，计算获得各浓度点的Ct响应值分布范围。均匀分布的抽样方法为从标准均匀分布R(0，1)中抽取随机数r，构造ξ=a+(b-a)r。MATLAB软件在拟合内源基因zSSIIb和59122片段的校准曲线时，矢量E、协方差矩阵V、u(m)、u(k)、u(k，m)和r(k，m)分别按照式(3)至(6)、式(8)至(9)计算。

16次重复测量试样获得zSSIIb基因和59122片段的Ct值(表3)概率分布遵守正态分布N[y，u2(y)][y为最佳估计值，u(y)为标准不确定度]，抽样方法为从标准正态分布N(0，1)中抽取随机数r，构造ξ=y+u(y)r。

3.3 MCM的传播

根据表1、表2、表3中的数据、“3.2”节中设定的各输入量概率分布以及抽样方法，编写MATLAB程序代码模拟运行抽样，获得106个随机样本以及每个输入量输入数学模型[式(1)至(9)]计算得到106个均匀分布或正态分布数据组，获得输出量的最佳估计值、标准不确定度及95%包含概率(coverage probability)下的包含区间(coverage interval)。

3.4 MCM的输出

MATLAB程序模拟Cti、Cti#、CtSi、CtSi#、m、k、QzSSIIb、Q59122和C等9个输入量的概率分布、抽样及通过数学模型式(1)、式(3)至(9)获得的平均值、标准不确定度和95%包含概率的包含区间。

表1 zSSIIb内源基因校准曲线拟合

表2 59122片段校准曲线拟合

表3 试样中内源基因zSSIIb和59122片段的Ct值

4 结果

由表4可知，zSSIIb基因和59122片段的校准曲线的斜率分别为-3.60和-3.29，截距分别为31.50和30.41，即zSSIIb基因和59122片段的校准曲线分别为y=-3.60x+31.50和y=-3.29x+30.41。QzSSIIb和Q59122既作为输入量又作为输出量的平均值分别为178.47和1.02，不确定度分别为1.99和0.03，95%包含概率的包含区间分别为174.59～182.41和0.96～1.08。获得的59122片段相对含量C的平均值、标准不确定度和95%包含概率的包含区间分别为 0.57%、1.82×10-4和0.54%～0.61%。此外，MATLAB程序运行产生的混合样品中转基因玉米59122的相对含量的概率密度分布为典型的“钟形”正态分布(图1)。

表4 各输入/输出量的概率分布、均值、不确定度和95%置信度的包含区间

续表4

影响因子或输入/输出量概率分布平均值标准不确定度95%包含概率的包含区间CtS1均匀分布23.410.0723.29～23.53CtS2均匀分布23.450.0723.33～23.57CtS3均匀分布23.480.0723.36～23.60CtS4均匀分布23.450.0723.33～23.57CtS5均匀分布23.450.0723.33～23.57CtS6均匀分布23.490.0722.37～22.61CtS7均匀分布23.450.0723.33～23.57CtS8均匀分布23.320.0723.20～23.44CtS9均匀分布23.400.0723.28～23.52CtS10均匀分布23.450.0723.33～23.57CtS11均匀分布23.500.0723.38～23.62CtS12均匀分布23.580.0723.46～23.70CtS13均匀分布23.560.0723.44～23.68CtS14均匀分布23.490.0723.37～23.61CtS15均匀分布23.510.0723.39～23.63CtS16均匀分布23.360.0723.24～23.48CtS1#均匀分布30.250.1729.97～30.53CtS2#均匀分布30.370.1730.09～30.65CtS3#均匀分布30.380.1730.10～30.66CtS4#均匀分布29.940.1729.67～30.22CtS5#均匀分布30.520.1730.25～30.81CtS6#均匀分布30.370.1730.09～30.65CtS7#均匀分布30.350.1730.07～30.63CtS8#均匀分布30.390.1730.11～30.67CtS9#均匀分布30.230.1729.95～30.51CtS10#均匀分布30.630.1730.35～30.91CtS11#均匀分布30.340.1730.06～30.62CtS12#均匀分布30.510.1730.23～30.79CtS13#均匀分布30.570.1730.29～30.85CtS14#均匀分布30.630.1730.35～30.91CtS15#均匀分布30.320.1730.04～30.60CtS16#均匀分布30.290.1730.01～30.57QzSSIIb正态分布178.471.99174.59～182.41Q59122正态分布1.020.030.96～1.08C正态分布0.005 71.82×10～40.005 4～0.006 1

注：Ct11～Ct53(统称Cti)表示zSSIIb基因测量中标准物质各浓度点的响应值；Ct11#～Ct53#(统称Cti#)表示59122玉米品系测量中标准物质各浓度点的响应值；CtS1～CtS16(统称CtSi)表示试样16次重复测量zSSIIb基因获得的Ct值；CtS1#～CtS16#(统称CtSi#)表示试样16次重复测量59122片段获得的Ct值；QzSSIIb表示试样中zSSIIb的绝对含量；Q59122表示试样中59122片段的绝对含量；C表示试样中59122片段的相对含量；mzSSIIb和m59122分别为zSSIIb基因和59122片段的校准曲线的斜率；kzSSIIb和k59122分别为zSSIIb基因和59122片段的校准曲线的截距。

5 讨论

本研究首次报道了采用蒙特卡洛法实施“概率分布传播”评定转基因成分的测量不确定度。转基因玉米59122的相对含量为0.57%，与其“理论含量”0.5%的相对偏差(bias)为14%，在国际广泛认可的转基因成分定量检测偏差范围(-25%≤bias≤25%)内。本次试验中，在95%的包含概率下转基因玉米59122的相对含量检测值落在0.54%～0.61%范围内， 不确定度极小(1.82×10-4)， 包含区间范围较窄，显示检测质量较高，数据可靠。混合样品中转基因玉米59122的相对含量的概率密度分布呈正态分布，与用GUM法评定测量不确定的假设条件“输出量概率密度函数为正态分布”一致，因此，揭示了转基因成分测量过程中各个输入量/输出量的概率分布及数据对称性等均满足用GUM法评定测量不确定度的假设条件，即转基因成分测量不确定度的评定可以采用GUM法。

在表4的9个变量的不确定度中，不确定度贡献最大的是玉米内源基因zSSIIb绝对含量的不确定度(1.99)，转基因玉米59122的相对含量C的不确定度最小(1.82×10-4)，显示出检测质量高、数据可靠。虽然玉米内源基因zSSIIb绝对含量的不确定度最大，但其对测量结果影响非常小。因此，将玉米内源基因zSSIIb绝对含量的概率分布输入测量模型传递到转基因玉米59122的相对含量C这个最终输出量时，测量不确定度已经非常小。此外，部分浓度(Ct51～Ct53、Ct11#～Ct13#、Ct41～Ct43、Ct41#～Ct43#和Ct21～Ct23)的标准物质稀释贡献的不确定度分布在0.28～0.40之间。总体上，这些不确定度居于本次测量不确定度的中等水平，反映出转基因成分测量不确定度的主要因素为标准物质的梯度稀释，而影响梯度稀释的主要原因是微量移液器的使用，这与采用GUM法评定转基因成分测量不确定度报道的主要因素为微量移液[10]一致。本研究的结果表明，影响转基因成分测量的主要因素为微量移液，因此在转基因成分的定量分析中应将提高微量液体体积转移的准确度作为首要考虑对象。