王茂元，賴銘勇，黃洪貴，吳妹英，田田，黃柳婷

（福建省淡水水產(chǎn)研究所，福建 福州 350002）

棘胸蛙Rana spinosa（俗稱石蛙、石鱗、石雞、石凍等）主要分布于中國(guó)南部以及越南北部丘陵山區(qū)[1]。棘胸蛙肉質(zhì)鮮美、口感爽滑、營(yíng)養(yǎng)豐富，具有較高的食用和醫(yī)用價(jià)值，是一種經(jīng)濟(jì)價(jià)值較高的名貴珍品。據(jù)《本草綱目》記載：“石蛙主治小兒癆瘦，疳疾最良，產(chǎn)婦尤佳”，市場(chǎng)需求量大[2]。由于人為捕殺和自然環(huán)境的惡化，野生棘胸蛙資源遭到嚴(yán)重破壞，已被中國(guó)物種紅色名錄列為易危（Vu）等級(jí)[3]。近年來(lái)，隨著棘胸蛙養(yǎng)殖技術(shù)逐步成熟，養(yǎng)殖業(yè)者為擴(kuò)大養(yǎng)殖規(guī)模，不斷進(jìn)行頻繁引種、近親繁殖、累代人工養(yǎng)殖等，存在產(chǎn)量下降、品質(zhì)退化、生物遺傳多樣性降低等諸多風(fēng)險(xiǎn)，嚴(yán)重危及棘胸蛙種質(zhì)資源。因此，開展棘胸蛙養(yǎng)殖群體遺傳多樣性研究，從分子水平分析其種群遺傳結(jié)構(gòu)特征和種群分化，對(duì)棘胸蛙的種質(zhì)資源保護(hù)與科學(xué)開發(fā)利用具有十分重要的意義。

測(cè)序基因分型技術(shù)（Genotyping-by-sequencing，GBS）[4]是近年來(lái)常見的單核苷酸多態(tài)性標(biāo)記（Single Nucleotide Polymorphism，SNP）開發(fā)最為有效、經(jīng)濟(jì)的簡(jiǎn)化基因組測(cè)序方法之一。其原理是采用限制性內(nèi)切酶加標(biāo)簽的方法，實(shí)現(xiàn)多樣本高通量平行測(cè)序[5]，獲取目的物種在基因組層面上的高密度SNP標(biāo)記。這些標(biāo)記可用于遺傳圖譜構(gòu)建[6,7]、種質(zhì)鑒定[8,9]、基因多樣性[10,11]及群體遺傳學(xué)[12,13]等多方面的研究。本研究通過(guò)利用GBS 簡(jiǎn)化基因組測(cè)序技術(shù)分析三明、龍巖和南平三個(gè)地區(qū)棘胸蛙群體的遺傳多樣性，旨在從分子水平揭示不同地域棘胸蛙的種群遺傳結(jié)構(gòu)，進(jìn)一步了解福建棘胸蛙的遺傳特征和種群分化，為棘胸蛙種質(zhì)資源保護(hù)和遺傳改良提供參考。

1 材料與方法

1.1 材料

2017 年10 月，于南平市正綠石蛙養(yǎng)殖有限公司取試驗(yàn)用棘胸蛙23 只、三明市龍?jiān)醇赝莛B(yǎng)殖場(chǎng)20 只以及龍巖市武平鑫旺質(zhì)勝農(nóng)業(yè)發(fā)展有限公司23 只，共計(jì)66 個(gè)樣品，平均體質(zhì)量為（143.79±26.55）g，均為當(dāng)?shù)匾吧苋鹤右淮? 齡蛙，分別命名為Rs-N、Rs-S 和Rs-L，采用低溫麻醉，使棘胸蛙進(jìn)入休眠狀態(tài)，活體解剖取后肢于95%乙醇中保存?zhèn)溆谩?/p>

1.2 方法

1.2.1 DNA 的提取與純化

采用SDS 法提取棘胸蛙后肢的肌肉組織全基因組DNA，在1%的瓊脂糖凝膠上100 V 電泳40 min，檢測(cè)樣本基因組DNA 的完整性。用核酸濃度測(cè)定儀（BioDrop-μLite）測(cè)定其濃度和純度，質(zhì)檢符合要求的高質(zhì)量DNA 樣本進(jìn)行后期的文庫(kù)構(gòu)建和測(cè)序。

1.2.2 GBS 文庫(kù)構(gòu)建

選擇質(zhì)檢合格的基因組DNA 樣本，采用限制性內(nèi)切酶ApekⅠ酶切基因組DNA。酶切后的片斷鏈接帶有內(nèi)切酶序列末端的普通接頭和探針接頭，再將不同樣品的鏈接產(chǎn)物集中，通過(guò)凝膠電泳回收180～480 bp 區(qū)間的DNA 片段，再進(jìn)行PCR 擴(kuò)增，用磁珠富集純法，構(gòu)建測(cè)序樣本文庫(kù)。

1.2.3 測(cè)序

構(gòu)建好的文庫(kù)用Agilent 2100 Bioanalyzer 和QPCR 的方法進(jìn)行質(zhì)控，通過(guò)質(zhì)控的文庫(kù)進(jìn)行測(cè)序。利用第二代測(cè)序技術(shù)（Next-Generation Sequencing，NGS），在Illumina HiSeq4000 測(cè)序平臺(tái)上進(jìn)行雙末端150 bp 測(cè)序。

1.3 數(shù)據(jù)統(tǒng)計(jì)分析

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

由測(cè)序所得的原始數(shù)據(jù)（Raw data）進(jìn)行質(zhì)控過(guò)濾，包括去除含有接頭序列的雙末端序列；去除N的含量超過(guò)該條序列長(zhǎng)度比例的10%時(shí)的雙末端序列；去除低質(zhì)量的序列（質(zhì)量值Q≤5 的堿基占整條序列的50%以上的被定義為低質(zhì)量序列）；去掉無(wú)法依據(jù)接頭區(qū)分樣品的序列等。過(guò)濾后的數(shù)據(jù)進(jìn)行后續(xù)分析。

1.3.2 聚類與比對(duì)

原始數(shù)據(jù)經(jīng)過(guò)過(guò)濾后獲得的數(shù)據(jù)為有效數(shù)據(jù)（Clean data）。將每個(gè)樣本的有效數(shù)據(jù)分別與擬參考基因組非洲爪蟾Xenopus laevis 進(jìn)行BWA 0.5.7 比對(duì)分析[14,15]，bwa mem 使用參數(shù)：-M-R-R"@RG；ID：wHAIPI016410-33；PL：Illumina；PU：150211_I19 1_FCC6L7EANXX_L5_wHAIPI016410-33；LB：wHAI PI016410-33；SM：TEST1；CN：BGI"all.chrs.con.faread1.fq.gz read2.fq.gz ＞aligned_reads.sam。統(tǒng)計(jì)比對(duì)完成后，應(yīng)用samtools eference{samtools1}軟件[16,17]根據(jù)基因組上的坐標(biāo)對(duì)SAM 文件進(jìn)行排序，并將其轉(zhuǎn)化成BAM 文件。

1.3.3 SNP 檢測(cè)

應(yīng)用GATK 檢測(cè)SNPs，應(yīng)用BGI 自主開發(fā)軟件對(duì)其進(jìn)行注釋和統(tǒng)計(jì)[18,19]，SNP 過(guò)濾參數(shù)為："QD＜2.0||FS＞60.0||MQ＜40.0||MQRankSum＜-12.5||Read-PosRankSum＜-8.0"。SNP 檢測(cè)的過(guò)程如下：（1）將所有樣本的比對(duì)文件bam 放到一個(gè)list 文件里面；（2）使用GATK 得到所有樣本的SNP；（3）使用GATK過(guò)濾得到的變異結(jié)果，選取可靠的變異結(jié)果。

1.3.4 群體遺傳學(xué)分析測(cè)定

基于各樣品的SNP 使用p-distance 方法計(jì)算各樣品的遺傳距離矩陣[20]，用NJ 法構(gòu)建群體間的系統(tǒng)發(fā)育樹。應(yīng)用Admixture 軟件基于最大似然法對(duì)66 份棘胸蛙樣品進(jìn)行主成分分析。

2 結(jié)果與分析

2.1 總DNA 的檢測(cè)

用1%瓊脂糖凝膠電泳檢測(cè)棘胸蛙樣本的總DNA（圖1）。結(jié)果顯示：基因組DNA 條帶單一清晰可見，無(wú)拖尾，表明基因組DNA 提取質(zhì)量良好，無(wú)明顯降解、無(wú)RNA 和蛋白質(zhì)污染，能夠滿足后期建庫(kù)測(cè)序要求。

2.2 數(shù)據(jù)分析

通過(guò)Illumina HiSeq 測(cè)序平臺(tái)，下機(jī)原始數(shù)據(jù)（Raw data）按照探針拆分、過(guò)濾后獲得有效數(shù)據(jù)（Clean data）共66.79 G。三區(qū)市66 個(gè)棘胸蛙樣本測(cè)序數(shù)據(jù)統(tǒng)計(jì)見表1。由表1 計(jì)算得出，66 個(gè)樣品的平均原始數(shù)據(jù)約為6.98 M，平均有效數(shù)據(jù)約為3.69 M，且Q20 均在97%以上，Q30 均高于93%，GC 含量平均為46.70%，表明測(cè)序數(shù)據(jù)質(zhì)量值高，能滿足后續(xù)信息分析。

2.3 比對(duì)與SNP 檢測(cè)

以非洲爪蟾為參考基因組，應(yīng)用BWA 軟件將所有的測(cè)序數(shù)據(jù)比對(duì)到參考基因組上。其基因組大小為2 734 636 505 bp，有效基因組大小為2 408 724 787 bp（參考序列中不含N），參考基因組GC 含量為34.34%。所有樣本的比對(duì)率介于58.40%和60.76%之間。

對(duì)3 個(gè)地區(qū)棘胸蛙所獲得的SNP 進(jìn)行統(tǒng)計(jì)匯總（表2），為確保鑒定得到的SNP 的質(zhì)量，對(duì)檢測(cè)到的SNP 進(jìn)行過(guò)濾：（1）異置信度與質(zhì)量深度的比值≥2；（2）使用Fisher's 精確檢驗(yàn)來(lái)檢測(cè)鏈偏倚現(xiàn)象而得到的Fhred 格式的p 值，本試驗(yàn)設(shè)置為≤60；（3）均方根值比對(duì)質(zhì)量值≥40；（4）變異非參檢驗(yàn)、秩和檢驗(yàn)的零-均值（Zscore）和比對(duì)序列質(zhì)量值的比值≥-12.5；（5）變異非參檢驗(yàn)、秩和檢驗(yàn)的零-均值（Zscore）和比對(duì)序列位置偏差的比值≥-8.0。最終經(jīng)過(guò)過(guò)濾，3 個(gè)地區(qū)的棘胸蛙樣本共獲得396 969個(gè)SNP 位點(diǎn)，其中龍巖群體獲得SNP 位點(diǎn)數(shù)最多，為139 169 個(gè)；三明群體最少，為123 437 個(gè)SNP 位點(diǎn)。經(jīng)過(guò)再次篩選檢測(cè)，3 個(gè)地區(qū)的棘胸蛙樣品共獲得8 615 個(gè)共有SNP 位點(diǎn)用于高級(jí)分析。

表3 3 個(gè)地區(qū)的棘胸蛙樣本基因分型結(jié)果統(tǒng)計(jì)Tab.3 Test result of genotypes of all samples spine frog Rana spinosa collected from three areas

由表3 可知，南平棘胸蛙群體（Rs-N）純合基因型的占比平均為99.77%，雜合基因型的占比平均為0.23%；三明棘胸蛙群體（Rs-S）純合基因型的占比平均為99.66%，雜合基因型的占比平均為0.34%；龍巖棘胸蛙群體（Rs-L）純合基因型的占比平均為99.65%，雜合基因型的占比平均為0.35%。3 個(gè)群體的純合基因型的占比都在99%以上，表明3 個(gè)群體種質(zhì)都比較純，遺傳多樣性不夠豐富；南平群體的純合性占比最高，雜合性最低，表明南平群體種質(zhì)更加純合，受到外源種質(zhì)的污染更低，種質(zhì)保存更好；龍巖群體雜合基因型占比較高，說(shuō)明其群體的種質(zhì)不夠純，但遺傳多樣性相對(duì)豐富。

2.4 群體分化分析

為探究3 個(gè)棘胸蛙群體之間的遺傳分化程度，基于各樣品的SNP 使用p-distance 方法計(jì)算各樣品的遺傳距離矩陣，使用TreeBest 軟件中的鄰接算法，經(jīng)過(guò)1 000 次迭代繪制系統(tǒng)NJ 發(fā)育樹（圖2）。

由圖2 可知：南平（Rs-N）23 個(gè)棘胸蛙群體中絕大部分個(gè)體聚成1 支，可以很好地與三明群體（Rs-S）和龍巖群體（Rs-L）區(qū)分開來(lái)，說(shuō)明與這兩個(gè)群體的進(jìn)化差異顯著，群體遺傳距離近，種質(zhì)比較純，更適合作為種質(zhì)選育材料；三明群體和龍巖群體個(gè)體間都不能完全聚成1 支，表現(xiàn)為群體中的個(gè)體相互摻雜聚成多個(gè)分支，呈現(xiàn)無(wú)規(guī)律分散，說(shuō)明這兩個(gè)地區(qū)的棘胸蛙群體各自進(jìn)化差異不顯著，群體種質(zhì)不如南平地區(qū)的種純。

2.5 群體主成分分析

根據(jù)個(gè)體水平的遺傳距離對(duì)3 個(gè)棘胸蛙群體進(jìn)行主成分分析（圖3）。由圖3 可知，3 個(gè)地區(qū)的棘胸蛙群體區(qū)分度不明顯，群體遺傳差異不顯著。

3 討論

群體遺傳學(xué)及其分子系統(tǒng)進(jìn)化等研究的基礎(chǔ)是遺傳標(biāo)記技術(shù)。從20 世紀(jì)80 年代至今，DNA分子標(biāo)記技術(shù)快速發(fā)展，由最初的限制性片段長(zhǎng)度多態(tài)性（Restriction fragment length polymorphism，RFLP）[21]到隨機(jī)性擴(kuò)增DNA（Random amplified polymorphic DNA，RAPD）[22,23]，到現(xiàn)在應(yīng)用十分廣泛的微衛(wèi)星DNA（Microsatellite DNA）[24]以及單核苷酸多態(tài)性（Single nucleotide polymorphism，SNP）[25]。這些傳統(tǒng)的分子標(biāo)記方法工作量大，比較繁瑣和耗時(shí)。近些年，隨著高通量測(cè)序技術(shù)（Next generation sequencing，NGS）的飛速發(fā)展，簡(jiǎn)化基因組測(cè)序（Reduced-representation sequencing）技術(shù)避免了上述方法的不足，開始逐漸被應(yīng)用于非模式生物的研究中[26-28]，成為一種十分有前景的分子標(biāo)記技術(shù)。

棘胸蛙作為一種經(jīng)濟(jì)價(jià)值較高的易危蛙類，逐漸被人們所認(rèn)識(shí)。國(guó)內(nèi)外關(guān)于棘胸蛙的研究主要集中在分類與分布、形態(tài)學(xué)、繁殖與發(fā)育生物學(xué)、細(xì)胞遺傳學(xué)以及人工養(yǎng)殖技術(shù)等方面，而從分子生物學(xué)角度探討棘胸蛙種群遺傳結(jié)構(gòu)和系統(tǒng)發(fā)生的研究還不多，并且多停留在線粒體DNA 標(biāo)記技術(shù)[29,30]以及傳統(tǒng)的分子遺傳標(biāo)記技術(shù)[31-33]。本研究突破傳統(tǒng)分子標(biāo)記技術(shù)，首次將GBS 技術(shù)應(yīng)用于棘胸蛙遺傳多樣性分析，在3 個(gè)地區(qū)的棘胸蛙群體中開發(fā)獲得8 615 個(gè)高質(zhì)量SNP，其中以龍巖群體的SNP 最多，為3 589 個(gè)，南平群體的SNP 最少，為2 153 個(gè)，說(shuō)明不同地區(qū)棘胸蛙群體SNP 數(shù)存在較大差異，這可能與其棲息的生態(tài)環(huán)境或人工干預(yù)不同有關(guān)；在同一群體不同個(gè)體間SNP 也不相同，這為今后開展品種選育提供更多的遺傳信息。

本研究中的群體進(jìn)化樹表明，南平棘胸蛙群體主要聚在一起，與其地理分布一致，不存在與三明和龍巖個(gè)體聚類時(shí)混合交叉的現(xiàn)象，能很好地與三明和龍巖群體區(qū)分開。三明和龍巖群體在聚類時(shí)，個(gè)體混合交叉聚類現(xiàn)象明顯，不能很好地相互區(qū)分開。PCA 主成分分析表明，3 個(gè)群體的區(qū)分度均不明顯，雖然在聚類時(shí)呈現(xiàn)一定程度的分化，但未形成一定的地理隔離現(xiàn)象。

總之，本研究首次采用GBS 簡(jiǎn)化基因組測(cè)序技術(shù)，揭示了福建省南平市、三明市和龍巖市3 個(gè)地區(qū)棘胸蛙群體的遺傳多樣性與遺傳結(jié)構(gòu)關(guān)系。南平種群的遺傳多樣性已表現(xiàn)出雜合度低、遺傳多樣性下降等趨勢(shì)。系統(tǒng)發(fā)生樹以及PCA 主成分分析發(fā)現(xiàn)，南平種群與其他2 個(gè)種群的種群分化表現(xiàn)出一定程度的遺傳分化，但還未形成地理隔離。此次研究結(jié)果為制定有效的棘胸蛙種質(zhì)保護(hù)措施，更好地保護(hù)與利用棘胸蛙種質(zhì)資源提供了參考。但由于本次試驗(yàn)受采樣困難、覆蓋度不高等客觀原因限制，研究結(jié)果還不夠深入，下一步將提高采樣密度，結(jié)合核基因標(biāo)記技術(shù)，進(jìn)行更深入、精細(xì)的研究分析，以期獲得更為詳細(xì)的棘胸蛙群體遺傳信息資料。