用蒙特卡洛方法評定轉基因玉米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]一致。本研究的結果表明，影響轉基因成分測量的主要因素為微量移液，因此在轉基因成分的定量分析中應將提高微量液體體積轉移的準確度作為首要考慮對象。