人工模擬條件下環(huán)境DNA 宏條形碼技術(shù)的定量分析初探*

牟 銘 李 昂 趙新寧 柳淑芳① 莊志猛

(1. 上海海洋大學(xué)水產(chǎn)與生命學(xué)院 上海 201306；2. 中國水產(chǎn)科學(xué)研究院黃海水產(chǎn)研究所 青島海洋科學(xué)與技術(shù)試點國家實驗室海洋漁業(yè)科學(xué)與食物產(chǎn)出過程功能實驗室 山東 青島 266071)

環(huán)境 DNA 宏條形碼技術(shù)(environmental DNA metabarcoding)可以通過高通量測序?qū)Νh(huán)境 DNA(eDNA)進(jìn)行多物種分類鑒定。隨著eDNA 宏條形碼技術(shù)的應(yīng)用和發(fā)展，該技術(shù)對群落中單一或組合物種的定量評估逐漸顯現(xiàn)出優(yōu)勢。無論是水族箱或中觀實驗等小型生態(tài)系統(tǒng)(Minamotoet al; Moyeret al, 2014;Pilliodet al, 2013; Thomsenet al, 2012; Doiet al,2017)，還是江河湖泊等自然淡水系統(tǒng)(Doiet al, 2017;Lacoursiererousselet al, 2016)，甚至是開放式的海洋系統(tǒng)(Joet al, 2017)，研究人員都嘗試過使用eDNA宏條形碼技術(shù)測算種群相對數(shù)量。大量研究表明，特定環(huán)境中物種的個體數(shù)量和基于eDNA 宏條形碼技術(shù)獲得的高通量測序reads 數(shù)之間具有明顯的相關(guān)性。Evans 等(2016)設(shè)計了一個包含魚類和兩棲類的中觀實驗，證明水體中物種個體數(shù)或生物量相對豐度與eDNA 宏條形碼技術(shù)獲得的高通量測序reads 數(shù)呈正相關(guān)；Hanfling 等(2016)對比較了一個天然湖泊的長期監(jiān)測魚類豐度結(jié)果與eDNA 宏條形碼技術(shù)得到的高通量測序reads 數(shù)之間的相關(guān)性，結(jié)果表明二者呈正相關(guān)；Thomsen 等(2016)研究發(fā)現(xiàn)，當(dāng)將魚類的eDNA 高通量測序reads 數(shù)聚類到“科”的分類等級時，eDNA 宏條形碼技術(shù)測得的深海生境中的生物量與傳統(tǒng)拖網(wǎng)捕撈數(shù)據(jù)比較吻合?？梢灶A(yù)見，eDNA 宏條形碼技術(shù)具有快速測算群落中物種豐度的潛能，將成為資源保護(hù)和管理中具有應(yīng)用前景的調(diào)查工具。

然而，通過對已發(fā)表的相關(guān)文獻(xiàn)進(jìn)行梳理和歸納，發(fā)現(xiàn)當(dāng)前大多數(shù)利用eDNA 高通量測序結(jié)果的reads 數(shù)來評估自然環(huán)境中生物相對數(shù)量的研究并不能得到明確的量化關(guān)系結(jié)果(Limet al, 2016)。究其原因，是因為從生物到eDNA 以及從eDNA 再到高通量測序結(jié)果的2 個過程中的數(shù)學(xué)關(guān)系難以理清。其中，實驗室操作過程中的主要問題是：(1)自然水體中，不同濃度 eDNA 富集效率不同(Kellyet al, 2016)；(2)PCR 擴增過程中，通用引物對不同物種DNA 存在偏倚性(Elbrechtet al, 2015; Pinolet al, 2015)。假定水體中的eDNA 全部回收，且PCR 擴增時不存在引物偏倚性，這種理想狀態(tài)下的水體中eDNA 組成與其高通量測序reads 數(shù)是否完全呈線性關(guān)系？為了探究這個問題，本研究在實驗室可控條件下，選擇2 個同屬近緣種，對其DNA 樣品進(jìn)行不同比例混合，模擬從自然水體中富集到的eDNA 復(fù)合樣品，既保證了樣品的回收率，又降低了引物偏倚的干擾。以此為模板，探究eDNA 宏條形碼技術(shù)檢測種群相對數(shù)量的準(zhǔn)確性，旨在為驗證eDNA 宏條形碼技術(shù)監(jiān)測水生生物資源量的可行性提供直接證據(jù)。

1 材料與方法

1.1 實驗方案設(shè)計

1.1.1 實驗樣品的選擇 鑒于通用引物對不同物種的擴增效率有所不同，本研究實驗樣品選取2 個同屬近緣物種用于制備模擬eDNA 樣品，以確保將PCR過程中因引物對物種DNA 的偏移引起的實驗誤差降低到最小程度。

實驗樣品選擇了2 種常見對蝦：凡納濱對蝦(Penaeus vannamei)和墨吉對蝦(Penaeus merguiensis)，它們同為節(jié)肢動物門(Arthropoda)、軟甲綱(Malacostraca)、十足目(Decapoda)、對蝦科(Penaeidae)、對蝦屬(Penaeus)，是較為常見的同屬近緣種。

1.1.2 模板DNA 的提取 采用傳統(tǒng)的酚-氯仿-異戊醇方法分別提取2 種對蝦肌肉組織DNA，經(jīng)純化后使用超微量核酸蛋白分析儀(GE，美國)檢測濃度，稀釋DNA 濃度至20 ng/μL (±1%)，-20℃保存?zhèn)溆谩?/p>

1.1.3 模擬eDNA 的制備 人工模擬eDNA 復(fù)合樣品，共設(shè)7 個實驗組：將凡納濱對蝦和墨吉對蝦的DNA 模板分別按照1∶1、1∶2、1∶4、1∶8、1∶10、1∶50 和1∶100 的比例均勻混合，且混合后DNA 終濃度一致。設(shè)陽性對照2 組：分別選用凡納濱對蝦和墨吉對蝦DNA 樣品。設(shè)陰性對照1 組：用純水代替DNA 樣品。每組樣品設(shè)3 個重復(fù)。

將上述制備的eDNA 復(fù)合樣品作為eDNA 宏條形碼檢測的模板DNA。

1.2 通用引物的選取

考慮到后續(xù)高通量測序及DNA 條形碼數(shù)據(jù)庫信息資源的豐富度，本研究選取了Leray 等(2013)推薦的長度為313 bp 的線粒體COⅠ基因序列作為待測eDNA 目標(biāo)片段，該片段長度適中、數(shù)據(jù)資源豐富、物種鑒定特異性好，且引物通用性強。通用引物序列見表1，由華大基因有限公司合成。

表1 COⅠ基因PCR 擴增引物信息Tab.1 Primer sets of COⅠ used for PCR amplification

1.3 PCR 擴增

PCR 擴 增 體 系：2×RapidTaqMaster Mix(Vazyme) 25 μL，上下游引物(mlCOⅠintF/jgHCO2198)(10mmol/L)各2 μL，DNA 模板5 μL，無菌水補齊總體積至50 μL。

為了增強擴增特異性，使用Touch down PCR 反應(yīng)程序：95℃預(yù)變性10 min；95℃變性30 s，62℃退火30 s，72℃延伸30 s，16 個循環(huán)；95℃變性30 s，50℃退火30 s，72℃延伸30 s，25 個循環(huán)；72℃延伸3 min；4℃保存。

用1%瓊脂糖凝膠電泳檢測PCR 產(chǎn)物質(zhì)量，選用DL2000 DNA marker，目標(biāo)條帶約為300 bp。挑選電泳條帶單一且明亮的PCR 產(chǎn)物用于高通量測序。

1.4 測序及結(jié)果注釋

將所得PCR 產(chǎn)物送至青島歐易生物科技有限公司，進(jìn)行高通量測序。采用Illumina Miseq 平臺測序。

按照97%相似性對非重復(fù)序列(不含單序列)進(jìn)行OTU(operational taxonomic units)聚類，在聚類過程中去除嵌合體(chimera，融合了2 個不同序列信息的序列)、異源雙鏈序列等干擾序列，得到OTU 聚類結(jié)果，然后對不同OTU 的代表序列進(jìn)行結(jié)果注釋。在結(jié)果注釋時，參照 NCBI 公共數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/)來進(jìn)行分類學(xué)信息的標(biāo)注，同時使用中國重要漁業(yè)生物DNA 條形碼信息平臺(http://www.fisherybarcode.cn)對注釋結(jié)果進(jìn)行輔助校正。

1.5 數(shù)據(jù)處理

統(tǒng)計每個樣品高通量測序的原始reads 數(shù)，以所有樣本中reads 數(shù)最低值為基準(zhǔn)，進(jìn)行抽平，使每個樣品測序結(jié)果的reads 數(shù)相同，用于后續(xù)數(shù)據(jù)分析。

在推導(dǎo)數(shù)學(xué)關(guān)系時，設(shè)：eDNA 樣品來自2 個物種，分別為A、B；實際參與PCR 反應(yīng)的底物模板DNA 的量為n；參與PCR 反應(yīng)的底物模板DNA 理論量為x；每個循環(huán)中，參與擴增的2 個物種底物模板DNA 比率為λ；擴增的循環(huán)數(shù)為c；PCR 產(chǎn)物的高通量測序reads 數(shù)為y。假定理想情況為：在PCR 擴增的不同循環(huán)中，參與擴增的2 個物種底物模板DNA比率固定；PCR 擴增全程都保持底數(shù)為2 的指數(shù)增長型(反應(yīng)物足夠)，則高通量測序reads 數(shù)可以反應(yīng)eDNA 樣品中2 個物種的相對量，即推得公式(1)：

經(jīng)變形得公式(2)：

2 結(jié)果與分析

2.1 PCR 擴增結(jié)果

以人工模擬的7 組eDNA 復(fù)合樣品及對照組樣品為模板，使用eDNA 宏條形碼通用引物進(jìn)行PCR 擴增，經(jīng)1%瓊脂糖凝膠電泳檢測，結(jié)果顯示，陽性對照和7 個實驗組均成功擴增出300 bp 左右的目標(biāo)片段，陰性對照無擴增產(chǎn)物(圖1)。

圖1 瓊脂糖凝膠電泳檢測結(jié)果Fig.1 Results of agarose gel electrophoresis

2.2 測序數(shù)據(jù)

2.2.1 高通量測序結(jié)果 本研究對7 個實驗組共21 個PCR 產(chǎn)物進(jìn)行了高通量測序。統(tǒng)計高通量測序獲得的序列長度分布，顯示99.80%以上序列的長度分布于301～350 bp 之間(圖2)。抽平后共計獲得有效序列674,205 條，序列平均長度為312.41 bp。

圖2 測序數(shù)據(jù)長度分布統(tǒng)計Fig.2 Result of distributional interval of the sequence length

2.2.2 分類學(xué)信息注釋及reads 數(shù)統(tǒng)計結(jié)果 將每個樣品reads 總數(shù)抽平后，對674,205 條序列進(jìn)行分類學(xué)信息注釋，保留有效序列(序列長度在301～350之間，且分類學(xué)信息注釋結(jié)果是墨吉對蝦和凡納濱對蝦)的reads。各實驗組高通量測序reads 數(shù)及其分類注釋結(jié)果詳見表2。實驗組I 的3 組平行樣共獲得94,774 條序列，每個樣品注釋出的墨吉對蝦和凡納濱對蝦平均reads 數(shù)分別為11,097 和20,494，二者的比值為13/24；實驗組Ⅱ的3 組平行樣共獲得95,208 條序列，二者reads 數(shù)比值為18/23；實驗組Ⅲ的3 組平行樣共獲得 95,572 條序列，二者 reads 數(shù)比值為77/73；實驗組Ⅳ的3 組平行樣共獲得95,250 條序列，二者reads 數(shù)比值為67/44；實驗組Ⅴ的3 組平行樣共獲得94,408 條序列，二者reads 數(shù)比值為43/24；實驗組Ⅵ的3 組平行樣共獲得94,640 條序列，二者reads數(shù)比值為 127/35；實驗組Ⅶ的 3 組平行樣共獲得95,526 條序列，二者reads 數(shù)比值為187/23。

表2 各實驗組獲得的序列數(shù)及分類學(xué)信息注釋Tab.2 Results of taxonomic information and numbers of reads

2.2.3 混合模板PCR 擴增的引物偏倚率 當(dāng)模板中2 個物種DNA 濃度比例為1∶1 (墨吉對蝦/凡納濱對蝦)時，2 個物種的模板DNA 含量相同，擴增過程中主要是引物對不同底物的偏好引起的擴增產(chǎn)物差別，此時主要是引物偏倚引起的reads 數(shù)差距。根據(jù)公式(2)變形得：

代入擴增循環(huán)數(shù)，由reads 平均數(shù)比值，計算得到每個循環(huán)中參與反應(yīng)的2 個物種模板DNA 的平均比率為66/67，換言之，PCR 擴增的每個循環(huán)中，平均偏倚率為1/67，約為1.5%，即該引物在同一PCR體系中，每次退火過程對墨吉對蝦DNA 模板的結(jié)合率比凡納濱對蝦DNA 模板大約少1.5%。

2.2.4 eDNA 組成與高通量測序reads 數(shù)的函數(shù)關(guān)系

以各實驗組設(shè)置的底物中墨吉對蝦DNA 與凡納濱對蝦DNA 的比值為橫坐標(biāo)，以各實驗組樣品高通量測序結(jié)果注釋的墨吉對蝦與凡納濱對蝦reads 平均數(shù)的比值為縱坐標(biāo)，繪制eDNA 組成與高通量測序reads 數(shù)的函數(shù)關(guān)系圖。在公式(2)的前提假設(shè)條件下，在高通量結(jié)果中，2 種對蝦的reads 數(shù)之比和初始PCR模板DNA 中2 種對蝦的DNA 比例應(yīng)呈現(xiàn)y=kx線性關(guān)系，對表2 中的7 個實驗組的reads 平均數(shù)比值和PCR 模板DNA 比值的2 組數(shù)據(jù)進(jìn)行一元線性回歸分析，最終得出eDNA 組成與高通量測序reads 數(shù)的函數(shù)關(guān)系：y=0.0716x+0.7043 (r2=0.9824)(圖3)。

圖3 2 個物種高通量測序結(jié)果的序列數(shù)比值與模板比例的關(guān)系Fig.3 The curve of ratio of two species'reads and ratio of template DNA

3 討論

利用eDNA 宏條形碼技術(shù)來評估水體中不同物種的相對數(shù)量，是該技術(shù)最具潛力的應(yīng)用前景之一。與傳統(tǒng)調(diào)查方法相比，eDNA 宏條形碼技術(shù)具有操作簡便、省時高效且對自然資源破壞小等特點。如何使用該技術(shù)來監(jiān)測水體中的生物量，一直是學(xué)界關(guān)注的焦點，但水體中的eDNA 難以完全富集，且富集到的eDNA 成分復(fù)雜；PCR 擴增時，引物對不同種類DNA模板的偏好還會導(dǎo)致擴增效率不同，這些問題亦是eDNA 宏條形碼技術(shù)推廣應(yīng)用的最大障礙。為運用eDNA 宏條形碼技術(shù)檢測環(huán)境生物相對數(shù)量而設(shè)計的中觀實驗，大都是在嘗試推導(dǎo)生物量與高通量測序結(jié)果reads 數(shù)的數(shù)學(xué)關(guān)系，主要包括環(huán)境中生物量與eDNA 量、環(huán)境中eDNA 量與提取的eDNA 量和提取的eDNA 模板與PCR 產(chǎn)物的高通量測序結(jié)果等3 個對應(yīng)關(guān)系，其中任何一個關(guān)系不明確，都將影響eDNA 宏條形碼技術(shù)的檢測結(jié)果，所以，必須逐步理清這3 種對應(yīng)關(guān)系的函數(shù)關(guān)系，進(jìn)而得到“高魯棒性”(模式體系固定且成熟，效果準(zhǔn)確且穩(wěn)定)的eDNA宏條形碼方法應(yīng)用體系。

分解困難是解決問題的有效方法之一，剖析DNA宏條形碼技術(shù)全過程，當(dāng)只關(guān)注富集到的eDNA 擴增及測序時，如果能創(chuàng)造一種eDNA 成分簡單、無引物偏好影響的理想條件，進(jìn)而探索eDNA 高通量測序reads 數(shù)和底物eDNA 含量之間的數(shù)學(xué)關(guān)系，便可以收獲“柳暗花明”的結(jié)果。

3.1 混合模板PCR 擴增的引物偏倚性

影響多物種eDNA 混合模板PCR 擴增效率的主要原因之一，便是PCR 過程中引物對不同模板DNA堿基組成存在偏好性(Kanagawa, 2003)。理論上，引物的每個脫氧核苷酸和模板DNA 鏈上同源區(qū)段的每個脫氧核苷酸可以實現(xiàn)一一對應(yīng)的關(guān)系，但實踐過程中通常難以實現(xiàn)，任何一個錯配位點的出現(xiàn)都會導(dǎo)致相應(yīng)的擴增效率降低，從而減少產(chǎn)物。為此，通用引物中加入了簡并堿基以降低錯配的影響。然而，在進(jìn)行復(fù)雜模板的PCR 擴增時，即使使用含簡并堿基的通用引物，引物結(jié)合能力也會因為不同模板DNA 引物結(jié)合區(qū)段序列的GC 堿基含量高低而不同(Acinaset al, 2005)，GC 含量高的引物結(jié)合力較強(Fonseca,2018)。為了提高模板DNA 引物結(jié)合區(qū)段堿基組成的一致性，從而降低引物偏倚的影響，本研究選擇較為常見的同屬近緣種凡納濱對蝦和墨吉對蝦作為實驗對象。同時，基于基因的進(jìn)化速率及其數(shù)據(jù)庫信息的豐富程度等因素，決定以mtDNA 中的COⅠ基因片段作為檢測的目的基因。

通過控制上述引物偏好的條件，最終高通量測序結(jié)果中引物的偏倚程度如何。研究分析發(fā)現(xiàn)，當(dāng)eDNA 模板中凡納濱對蝦和墨吉對蝦DNA 組成比例為1∶1 時，高通量測序獲得對應(yīng)物種reads 數(shù)比例為13/24。根據(jù)擴增循環(huán)數(shù)和公式(2)，推導(dǎo)此時參與擴增的2 種DNA 比約為66/67，即PCR 擴增的平均偏倚率為1/67，約1.5%。可見選擇凡納濱對蝦和墨吉對蝦制備簡單eDNA 混合模板產(chǎn)生的引物偏倚程度影響較小，符合最初的實驗設(shè)計預(yù)期。此時，2 個物種的模板濃度相同，高通量測序結(jié)果的偏倚主要是由于引物對不同模板的“偏好”程度不同產(chǎn)生的。在后續(xù)eDNA 宏條形碼技術(shù)的數(shù)量性研究探索中，設(shè)計實驗時可選用近源種作為實驗對象，從而降低擴增過程中引物對不同模板偏好性不同而導(dǎo)致的“偏倚”現(xiàn)象。

3.2 eDNA 組成與高通量測序reads 數(shù)的相關(guān)性

當(dāng)面對實際環(huán)境樣本時，如果直接應(yīng)用eDNA 宏條形碼技術(shù)分析高通量測序reads 數(shù)與待測樣本生物量的關(guān)系，從生物到eDNA 以及從eDNA 再到高通量測序結(jié)果的兩個過程中的數(shù)學(xué)關(guān)系則難以理清。鑒于此，本研究采取“化繁為簡”的策略，主要關(guān)注eDNA與高通量測序reads 數(shù)之間的關(guān)系。同時，考慮到復(fù)合底物的PCR 過程比較復(fù)雜(Harperet al, 2019)，最終注釋到每個物種的測序reads 數(shù)都受其反應(yīng)體系的影響，故直接比較樣本間的reads 數(shù)會產(chǎn)生一定偏差。充分考慮這些影響因素，本研究采用了樣本內(nèi)的2 個物種reads 數(shù)比值，用于確定其與模板eDNA 濃度的相關(guān)性。

統(tǒng)計分析實驗組的reads 平均數(shù)比值與PCR 模板eDNA 比值2 組數(shù)據(jù)，進(jìn)行一元線性回歸分析，最終得出eDNA 組成與高通量測序reads 數(shù)的函數(shù)關(guān)系：y= 0.0716x+ 0.7043 (r2=0.9824)。等式中的斜率可反應(yīng)出PCR 擴增過程中的偏倚程度(主要是引物偏倚)，當(dāng)擴增循環(huán)數(shù)相同，斜率越小，擴增過程中偏倚程度越大。等式中的截距主要是實驗過程中的誤差或者污染導(dǎo)致的。該等式可描述為：當(dāng)環(huán)境水樣中只存在墨吉對蝦和凡納濱對蝦的DNA 時，可通過高通量測序的reads 數(shù)推算出2 個物種的相對生物量。

在實際操作中，可能由于模板eDNA 的回收率存在差異，使之不能完全符合預(yù)設(shè)的比例；另外，PCR擴增時并不是全程以底數(shù)為2 的指數(shù)式增長，且存在其他不可避免的系統(tǒng)誤差，這些因素都會導(dǎo)致上述線性關(guān)系式出現(xiàn)一定偏離。但實驗結(jié)果已充分表明，高通量測序reads 數(shù)與初始模板eDNA 含量呈線性正相關(guān)。