基于三代測(cè)序技術(shù)的凡納濱對(duì)蝦(Litopenaeus vannamei)高通量候選基因關(guān)聯(lián)分析方法的建立*

李碧瀚 于 洋 劉桂嘉 羅 正 李富花①

(1. 中國(guó)科學(xué)院海洋研究所實(shí)驗(yàn)海洋生物學(xué)重點(diǎn)實(shí)驗(yàn)室 青島 266071; 2. 中國(guó)科學(xué)院大學(xué) 北京 100049)

水產(chǎn)動(dòng)物的經(jīng)濟(jì)性狀大多是由微效多基因控制的數(shù)量性狀, 如生長(zhǎng)速率、抗病性等(Yue, 2014;Gjedrem, 2015)。然而數(shù)量性狀易受環(huán)境等因素影響,導(dǎo)致了遺傳評(píng)估準(zhǔn)確率偏低?？剐孕誀钣捎诒硇投攘坷щy、遺傳力低, 加之與生長(zhǎng)或其他性狀之間存在負(fù)相關(guān), 使得利用傳統(tǒng)育種技術(shù)對(duì)抗病性狀進(jìn)行選育的效率較低(Cocket al, 2009)。隨著分子標(biāo)記, 特別是以單核苷酸多態(tài)性(SNP)標(biāo)記為主的基因分型技術(shù)的發(fā)展, 分子標(biāo)記輔助選擇(MAS)和基因組選擇(GS)等分子育種技術(shù)在水產(chǎn)動(dòng)物抗病育種中彰顯了巨大優(yōu)勢(shì), 比如個(gè)體評(píng)估的準(zhǔn)確性高, 可在繁殖周期的早期進(jìn)行選擇等。盡管分子育種技術(shù)優(yōu)勢(shì)大, 但是在大多數(shù)水產(chǎn)動(dòng)物的商業(yè)化育種中的應(yīng)用仍然較少。在水產(chǎn)動(dòng)物中開(kāi)展高效分子育種面臨著兩大難點(diǎn), 一者是缺少高密度、高精度的性狀連鎖標(biāo)記(Meuwissenet al,2004), 二者是對(duì)標(biāo)記和性狀之間關(guān)聯(lián)程度的評(píng)估有待優(yōu)化(VanRadenet al, 2009)。因此鑒定與性狀緊密相關(guān)的基因或基因位點(diǎn)是水產(chǎn)動(dòng)物分子育種亟需解決的問(wèn)題。

鑒定性狀相關(guān)位點(diǎn)的主要方法包括全基因組關(guān)聯(lián)分析(GWAS)和候選基因關(guān)聯(lián)分析(Xuet al, 2009)。其中, GWAS 已經(jīng)在人類(lèi)疾病相關(guān)基因的鑒定和農(nóng)作物、畜禽經(jīng)濟(jì)性狀相關(guān)SNP 鑒定中得到廣泛應(yīng)用, 而水產(chǎn)養(yǎng)殖物種中GWAS 研究起步較晚, 主要應(yīng)用于羅非魚(yú)、大西洋鮭魚(yú)、鯉魚(yú)和太平洋牡蠣等幾種研究基礎(chǔ)較好的水產(chǎn)物種(Houstonet al, 2020)。與GWAS 相比,候選基因關(guān)聯(lián)分析更加便捷、高效, 在魚(yú)類(lèi)、貝類(lèi)與蝦蟹類(lèi)等多種水產(chǎn)經(jīng)濟(jì)物種中得到了較好的應(yīng)用(Yue,2014; Yuet al, 2020)。隨著測(cè)序技術(shù)的不斷發(fā)展, 科研人員已成功破譯了多個(gè)水產(chǎn)動(dòng)物的全基因組(Lienet al, 2016; Zhanget al, 2019b; Houstonet al, 2020), 并構(gòu)建了高密度遺傳連鎖圖譜(Yuet al, 2015; Yueet al,2017), 篩選到大量分子標(biāo)記, 為候選基因關(guān)聯(lián)分析和性狀相關(guān)標(biāo)記的篩選奠定了重要基礎(chǔ)。

然而對(duì)水產(chǎn)動(dòng)物候選基因進(jìn)行關(guān)聯(lián)分析仍然缺乏高效、低成本的分型方法, 制約了相關(guān)技術(shù)的應(yīng)用。目前, 主要的基因分型方法包括基因芯片、SNaPshot 法、Taqman 分型、MassArray 飛行質(zhì)譜法、PCR 產(chǎn)物一代測(cè)序法、二代靶向測(cè)序法等。芯片法需要單獨(dú)定制芯片, 造價(jià)高; SNaPshot 法、Taqman 分型、MassArray 飛行質(zhì)譜法可以對(duì)候選基因進(jìn)行分型,但是只能分析已知的SNP 位點(diǎn), 且單個(gè)標(biāo)記的分型成本仍然較高; 一代測(cè)序分型是指通過(guò)PCR 產(chǎn)物直接測(cè)序分型, 優(yōu)勢(shì)在于精確度高、速度快, 但是其分型結(jié)果受DNA 序列結(jié)構(gòu)影響大、讀長(zhǎng)有限(<1 000 bp),且測(cè)序通量低、成本高; 二代靶向測(cè)序指通過(guò)探針捕獲或者PCR 擴(kuò)增目標(biāo)片段, 之后對(duì)該片段進(jìn)行高通量測(cè)序, 該技術(shù)能夠發(fā)掘新的SNP 位點(diǎn), 且價(jià)格相對(duì)較低, 但是對(duì)于長(zhǎng)的序列片段捕獲探針設(shè)計(jì)成本高,且PCR 建庫(kù)過(guò)程復(fù)雜。近年來(lái), 隨著三代測(cè)序技術(shù)的發(fā)展, 基于CCS 的測(cè)序技術(shù)準(zhǔn)確率不斷提高, 加大測(cè)序深度至15×后, HIFI 數(shù)據(jù)準(zhǔn)確率可達(dá)99.3% (Eidet al, 2009), 且由于三代測(cè)序具有長(zhǎng)度長(zhǎng)(平均5 000 bp), 無(wú)偏好性、準(zhǔn)確性高等特點(diǎn), 使得利用三代測(cè)序進(jìn)行基因分型成為可能(Rhoadset al, 2015), 尤其是利用其長(zhǎng)讀長(zhǎng)以及無(wú)系統(tǒng)偏差和無(wú)GC 偏好性的特點(diǎn),無(wú)需對(duì)目的序列進(jìn)行分段擴(kuò)增以及花費(fèi)大量精力組裝拼接測(cè)序結(jié)果, 因此在解決復(fù)雜基因組的基因分型方面具有明顯優(yōu)勢(shì)。

凡納濱對(duì)蝦(Litopenaeus vannamei)具有生長(zhǎng)速度快、抗逆能力強(qiáng)、肉質(zhì)鮮美且出肉率高等特點(diǎn), 是世界公認(rèn)的優(yōu)良對(duì)蝦養(yǎng)殖品種。然而, 隨著養(yǎng)殖環(huán)境的惡化, 細(xì)菌病、病毒病的頻發(fā), 嚴(yán)重制約了對(duì)蝦產(chǎn)業(yè)的健康可持續(xù)發(fā)展。其中, 帶有毒素質(zhì)粒的副溶血弧菌(Vibrio parahaemolyticus)所引發(fā)的對(duì)蝦早期死亡綜合征(Early Mortality Syndrome, EMS), 又稱為急性肝胰腺壞死癥(Acute Hepatopancreatic Necrosis Disease, AHPND) (Leeet al, 2015), 具有發(fā)病迅速、傳染性強(qiáng)、致死率高的特點(diǎn), 給對(duì)蝦養(yǎng)殖業(yè)造成了巨大損失?？共∑贩N的選育被認(rèn)為是解決病害問(wèn)題的根本途徑。然而由于抗病性狀的表型難以度量, 且抗病性狀的遺傳力通常較低, 受環(huán)境的影響較大, 利用傳統(tǒng)選育技術(shù)進(jìn)行選育的準(zhǔn)確性低。分子育種在抗病品種的選育中具有其獨(dú)特優(yōu)勢(shì), 抗病分子標(biāo)記和基因的篩選是進(jìn)行分子育種的基礎(chǔ), 目前有關(guān)對(duì)蝦抗病性狀相關(guān)基因和標(biāo)記的篩選進(jìn)展相當(dāng)緩慢, 缺乏高通量、低成本的基因分型技術(shù)是一個(gè)重要的限制因素。本團(tuán)隊(duì)前期通過(guò)對(duì)目標(biāo)序列進(jìn)行混池測(cè)序分型的方式, 篩選到了多個(gè)與對(duì)蝦抗弧菌性狀相關(guān)的SNP 標(biāo)記和基因, 其中LvPI3K的一個(gè)SNP 標(biāo)記與抗弧菌性狀呈顯著相關(guān)關(guān)系(Zhanget al, 2019a)。

本研究基于三代測(cè)序技術(shù), 在凡納濱對(duì)蝦中建立了一種基于三代靶向測(cè)序的候選基因關(guān)聯(lián)分析技術(shù), 并對(duì)前期篩選的抗弧菌性狀相關(guān)基因LvPI3K進(jìn)行了標(biāo)記發(fā)掘和關(guān)聯(lián)分析, 以篩選抗弧菌性狀相關(guān)標(biāo)記。本研究不僅為水產(chǎn)動(dòng)物提供了一種高效、低成本的候選基因關(guān)聯(lián)分析方法, 也為凡納濱對(duì)蝦的抗弧菌品種選育提供了有效的分子標(biāo)記。

1 材料與方法

1.1 實(shí)驗(yàn)材料

1.1.1 候選基因關(guān)聯(lián)分析材料的獲得 凡納濱對(duì)蝦遺傳材料來(lái)自渤海水產(chǎn)育種(海南)有限公司(海南,文昌), 分析用的群體為經(jīng)多代選育的多個(gè)弧菌抗性家系和敏感家系的雜交F2代群體。

隨機(jī)選擇500 尾個(gè)體, 利用副溶血弧菌進(jìn)行浸泡感染。感染所用的副溶血弧菌為從患病個(gè)體中分離出的副溶血弧菌菌株, 在含有2%氯化鈉的胰蛋白酶大豆肉湯(TBS)液體培養(yǎng)基中30 °C 恒溫培養(yǎng)過(guò)夜。經(jīng)對(duì)毒素質(zhì)粒PirA、PirB的PCR 擴(kuò)增檢測(cè)為陽(yáng)性, 確定其為致病性的副溶血弧菌(Hanet al, 2015)。采用血球計(jì)數(shù)板在光學(xué)顯微鏡下進(jìn)行計(jì)數(shù), 計(jì)算菌液濃度。設(shè)定弧菌浸泡感染濃度為5×106CFU/mL。弧菌感染實(shí)驗(yàn)持續(xù)8 d, 取最先死亡的96 尾蝦為敏感組樣品,最后存活的96 尾蝦為抗性組樣品, 記錄其存活和死亡時(shí)間, 并將樣品保存在–80 °C 直至DNA 提取。采用TIANGEN 植物DNA 提取試劑盒(TIANGEN, 北京)逐個(gè)提取DNA, 采用NanoDrop 1000 分光光度計(jì)(NanoDrop, 美國(guó))和瓊脂糖凝膠電泳以檢查DNA 的濃度和質(zhì)量。

1.1.2 驗(yàn)證群體材料的獲得 驗(yàn)證群體為經(jīng)多代選育的多個(gè)抗弧菌家系和敏感家系的雜交F2代群體,驗(yàn)證群體抗性與敏感材料的獲得方法與1.1.1 所示相同, 并且擴(kuò)大了驗(yàn)證群體樣本, 選擇弧菌感染后最先死亡的160 尾蝦為敏感材料, 最后死亡的160 尾蝦為抗性材料。DNA 提取方式同1.1.1。

1.2 LvPI3K 基因組DNA 擴(kuò)增

根據(jù)凡納濱對(duì)蝦基因組參考序列(Zhanget al,2019b), 獲得LvPI3K的基因組全長(zhǎng)序列, 共計(jì)7 378 bp, 使用primer 3 設(shè)計(jì)LvPI3K特異性擴(kuò)增引物。為實(shí)現(xiàn)在一個(gè)單分子實(shí)時(shí)測(cè)序單元(Single Molecule Real Time cell, SMRT cell)中區(qū)分每一個(gè)抗性個(gè)體和敏感個(gè)體的序列信息, 通過(guò)兩步PCR 的方式進(jìn)行擴(kuò)增。

第1 步PCR: 以提取的抗弧菌和敏感個(gè)體DNA為模板, 在LvPI3K的特異引物的5'端加上通用引物序列作為第一步PCR 的復(fù)合引物(圖1, 表1), 并將該引物的5'末端用5AmMC6 修飾封閉。采用PrimeSTAR GXL 聚合酶(TaKaRa, 日本)進(jìn)行擴(kuò)增, 擴(kuò)增體系為:50 ng 模板DNA, 10 pmol 引物F/R, 4 μL dNTP 混合物,10 μL 5×buffer, 1 μL GXL 酶, 滅菌水補(bǔ)充體系至50 μL。PCR 循環(huán)參數(shù)為: 98 °C 10 s, 60 °C 15 s, 68 °C 80 s, 共30 個(gè)循環(huán)。

表1 用于兩步PCR 的寡核苷酸引物Tab.1 Information of the primers used for 2-step PCR

第2 步PCR: 引入每個(gè)樣本專(zhuān)屬的barcode, 以便在PacBio RSⅡ上進(jìn)行多樣品平行測(cè)序。以第1 步PCR 產(chǎn)物稀釋50 倍為模板, 引物為5′端帶有獨(dú)特barcode 序列的通用引物序列(圖1)。使用PrimeSTAR GXL 聚合酶(TaKaRa, 日本)進(jìn)行擴(kuò)增, 體系為: 1 μL稀釋50×的第一步PCR 的產(chǎn)物, 10 pmol 引物F/R,4 μL dNTP 混合物, 10 μL 5×buffer, 1 μL GXL 酶, 滅菌水補(bǔ)充體系至50 μL。PCR 循環(huán)參數(shù)為: 98 °C 10 s,60 °C 15 s, 68 °C 80 s, 共18 個(gè)循環(huán)。

圖1 兩步PCR 將barcode 整合進(jìn)擴(kuò)增產(chǎn)物Fig.1 Incorporation of barcodes into the PCR amplicon via a two-step PCR approach

1.3 SMRT 文庫(kù)的構(gòu)建及PacBio 測(cè)序

將192 個(gè)個(gè)體的第二輪PCR 擴(kuò)增產(chǎn)物等摩爾混合后, 樣品送至天津生物芯片有限公司進(jìn)行SMRT 文庫(kù)的構(gòu)建及PacBio 測(cè)序, 文庫(kù)構(gòu)建方式參考SMRT bell Express Template Prep Kit 2.0 的建庫(kù)方案。

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

依據(jù) PacBio 的數(shù)據(jù)分析流程, 使用 SMRT link9.0 對(duì)測(cè)序原始數(shù)據(jù)進(jìn)行過(guò)濾, 獲得CCS 的原始數(shù)據(jù)。之后使用Lima 對(duì)FASTA 格式的原始CCS 數(shù)據(jù)進(jìn)行barcode 拆分, 輸出FASTQ 格式文件, 使用DeepVariant 將該文件與參考基因組文件比對(duì), 并進(jìn)行變異位點(diǎn)的發(fā)掘(Poplinet al, 2018), 使用WhatsHap 構(gòu)建單倍型后(Martinet al, 2016), 再次使用DeepVariant 軟件根據(jù)構(gòu)建的單倍型發(fā)掘變異位點(diǎn),獲得結(jié)構(gòu)變異(圖2)。分型數(shù)據(jù)使用VCFtools 轉(zhuǎn)換為plink 格式, 之后結(jié)合抗性和敏感的表型數(shù)據(jù), 使用R軟件進(jìn)行Logistic 關(guān)聯(lián)分析并作圖, 使用Haploview進(jìn)行單倍型分析(Barrettet al, 2005)。

圖2 SMRT 測(cè)序數(shù)據(jù)分析流程圖Fig.2 Flowchart of the SMRT sequencing analysis

使用FastQC 與MultiQC 軟件對(duì)SMRT 測(cè)序數(shù)據(jù)進(jìn)行整體的質(zhì)量分析。使用SAMtools 及PLINK 軟件對(duì)FASTQ 文件進(jìn)行位點(diǎn)測(cè)序深度與reads 長(zhǎng)度的統(tǒng)計(jì), 并使用R 軟件進(jìn)行做圖。

1.5 對(duì)蝦抗弧菌性狀相關(guān)候選位點(diǎn)的驗(yàn)證

抗性群體和敏感群體等實(shí)驗(yàn)材料的選擇, 對(duì)于突變位點(diǎn)的篩選與關(guān)聯(lián)分析有著直接的影響。為了排除群體差異, 本研究選擇了另外的遺傳背景豐富的群體作為驗(yàn)證材料, 對(duì)三代測(cè)序的分型結(jié)果加以驗(yàn)證, 主要針對(duì)兩個(gè)與抗弧菌性狀極顯著相關(guān)(P<0.01)的位點(diǎn)(PI3K_5366, PI3K_2205), 基于SNP 的側(cè)翼序列設(shè)計(jì)引物, 在320 個(gè)驗(yàn)證個(gè)體中進(jìn)行PCR 擴(kuò)增和測(cè)序。使用PrimeSTAR GXL 聚合酶(TaKaRa, 日本)進(jìn)行擴(kuò)增, 體系組成參考1.2 所述, PCR 擴(kuò)增參數(shù)為:98 °C 10 s, 58 °C 15 s, 68 °C 45 s, 共30 個(gè)循環(huán)。PCR產(chǎn)物送青島睿博生物技術(shù)有限公司進(jìn)行Sanger 測(cè)序,基于測(cè)序峰圖確定每個(gè)樣品的基因型, 分型結(jié)果使用R 軟件進(jìn)行關(guān)聯(lián)分析。

2 結(jié)果

2.1 PacBio 測(cè)序數(shù)據(jù)分析

目的片段LvPI3K基因全長(zhǎng)擴(kuò)增子序列長(zhǎng)度為7 378 bp (圖3a), 采用本實(shí)驗(yàn)室構(gòu)建的針對(duì)三代測(cè)序數(shù)據(jù)的專(zhuān)用分析流程, 對(duì)測(cè)序結(jié)果進(jìn)行了全面分析。統(tǒng)計(jì)顯示有99.89%的位點(diǎn)測(cè)序深度大于20 000×, 位點(diǎn)平均測(cè)序深度為22 012×。對(duì)樣品的測(cè)序reads 總數(shù)進(jìn)行統(tǒng)計(jì), 結(jié)果顯示, 測(cè)序的192 個(gè)個(gè)體共產(chǎn)生了22 089 條reads, 其中有21 990 條(99.6%) reads 長(zhǎng)度在目的片段長(zhǎng)度周?chē)?7.4—7.6 kb) (圖3b)。通過(guò)FastQC與MultiQC 軟件對(duì)測(cè)序數(shù)據(jù)質(zhì)量進(jìn)行分析, 發(fā)現(xiàn)除極少數(shù)序列外, 三代測(cè)序數(shù)據(jù)均符合質(zhì)控標(biāo)準(zhǔn)(Q30>90),并且全部樣品的平均測(cè)序深度為139×, 最高達(dá)902×,測(cè)序質(zhì)量和測(cè)序深度均可以滿足后續(xù)變異發(fā)掘的要求(Edgeet al, 2019; Wenger et al, 2019)。

圖3 LvPI3K 基因PCR 擴(kuò)增產(chǎn)物的瓊脂糖凝膠電泳圖(a)以及測(cè)序reads 長(zhǎng)度分布特征圖(b)Fig.3 Detection of PCR amplification products of LvPI3K by agarose gel electrophoresis (a) and length distribution of sequencing reads (b)

2.2 突變位點(diǎn)檢測(cè)及基因分型

使用變異調(diào)用軟件DeepVariant 在LvPI3K基因內(nèi)共篩選到91 個(gè)突變位點(diǎn), 最小等位基因頻率結(jié)果如圖 4 所示, 以最小等位基因頻率(Minor Allele Frequency, MAF)小于0.05 作為標(biāo)準(zhǔn)進(jìn)行稀有突變位點(diǎn)過(guò)濾, 最后得到65 個(gè)突變位點(diǎn)。這些位點(diǎn)中包括1個(gè)三等位SNP 位點(diǎn), 11 個(gè)indel, 其余均為二等位SNP位點(diǎn)。對(duì)SNP 的分布情況進(jìn)行分析, 發(fā)現(xiàn)內(nèi)含子區(qū)域內(nèi)的突變有36 個(gè), 在5'-UTR 區(qū)域的突變有24 個(gè), 其余4 個(gè)位于外顯子上。

圖4 質(zhì)控前篩選到的所有突變位點(diǎn)MAF 分布直方圖Fig.4 The minor allele frequency (MAF) distribution of the SNPs before quality control

2.3 候選基因關(guān)聯(lián)分析

為了降低由于不同關(guān)聯(lián)分析策略而導(dǎo)致的分析結(jié)果差異, 我們利用存活情況(閾值性狀)與存活時(shí)間(數(shù)量性狀)兩個(gè)指標(biāo)對(duì)篩選到的突變位點(diǎn)進(jìn)行關(guān)聯(lián)分析。經(jīng)過(guò)Logistic 關(guān)聯(lián)分析, 最終篩選到21 個(gè)與抗弧菌性狀顯著相關(guān)(P<0.05)的位點(diǎn)(表2), 其中有5 個(gè)SNPs 位于5′-UTR, 1 個(gè)SNP 位于外顯子區(qū)域, 12 個(gè)SNPs 位于內(nèi)含子, 3 個(gè)Indel 突變也均位于內(nèi)含子區(qū)域。在21 個(gè)抗性相關(guān)位點(diǎn)中, 有兩個(gè)位點(diǎn)PI3K_5366和PI3K_2205 與抗弧菌性狀呈極顯著相關(guān)(P<0.01)(圖5a, 5b)。

圖5 LvPI3K 基因上與抗弧菌性狀呈極顯著相關(guān)(P<0.01)的標(biāo)記及周邊標(biāo)記的P 值分布圖Fig.5 The P values of very significant markers related to the resistance of shrimp against V. parahaemolyticus (P<0.01) and surrounding markers in LvPI3K

表2 對(duì)蝦抗弧菌性狀相關(guān)位點(diǎn)(P<0.05)在抗性群體和敏感群體中的分布Tab.2 Distribution of the 21 SNPs (P<0.05) related to the resistance of shrimp against V. parahaemolyticus in resistant and susceptive groups

為進(jìn)一步了解抗性標(biāo)記之間的關(guān)系, 對(duì)鑒定到的SNP 標(biāo)記進(jìn)行了單倍型分析, 通常R2>0.8 的連鎖才被認(rèn)為是可信的單倍域(De Bakkeret al, 2005), 以此為標(biāo)準(zhǔn), 在LvPI3K基因中鑒定出了5 個(gè)具有高LD(D'>0.9,R2>0.8)的區(qū)域(圖6)。之后對(duì)這5 個(gè)緊密連鎖的單倍域與抗病性狀進(jìn)行關(guān)聯(lián)分析, 發(fā)現(xiàn)Block3、Block4 與Block5 的部分基因型在抗弧菌/敏感群體存在顯著差異(P<0.05) (表3)。尤其是由PI3K_5031_A/T和PI3K_5366_T/C 兩個(gè)位點(diǎn)構(gòu)成的Block5, 基因型AT/TC 在群體中的累計(jì)頻率達(dá)0.977, 可以作為凡納濱對(duì)蝦弧菌抗病性狀分子育種中一個(gè)潛在可用的單倍型。

表3 不同單倍型在對(duì)蝦抗弧菌群體和敏感群體中的分布情況Tab.3 Distributions of haplotypes in resistant and susceptive shrimp groups against V. parahaemolyticus

圖6 LvPI3K 基因上標(biāo)記間連鎖不平衡分析圖Fig.6 Linkage disequilibrium of markers in LvPI3K gene

2.4 對(duì)蝦抗弧菌性狀相關(guān)候選基因位點(diǎn)的驗(yàn)證

在320 個(gè)驗(yàn)證樣本中, 分別對(duì)兩個(gè)最顯著的位點(diǎn)PI3K_2205 和 PI3K_5366 進(jìn)行了驗(yàn)證, 其中位點(diǎn)PI3K_2205 為插入缺失突變, 位點(diǎn)PI3K_5366 位于第6 個(gè)外顯子上, 為同義突變。Sanger 測(cè)序分型結(jié)果顯示PI3K_2205 在驗(yàn)證群體中存在插入缺失, 但是關(guān)聯(lián)分析顯示其與抗弧菌性狀不相關(guān)。PI3K_5366 標(biāo)記的基因型與三代測(cè)序基因型分布一致(表 4), 其中,PI3K_5366 標(biāo)記的TT 基因型在敏感群體中出現(xiàn)的頻率顯著高于抗性群體, CC 基因型則反之; 進(jìn)一步對(duì)T、C 的基因頻率進(jìn)行統(tǒng)計(jì)的結(jié)果顯示, 在驗(yàn)證材料中, T、C 等位基因在抗性群體與敏感群體中出現(xiàn)的頻率存在極顯著差異(P<0.01), 等位基因型C 為抗性優(yōu)勢(shì)基因, T 則為敏感等位基因型。

表4 PI3K_2205 和PI3K_5366 位點(diǎn)在驗(yàn)證群體與三代測(cè)序分型群體中的基因分型情況對(duì)比Tab.4 Comparison of the genotyping on the sites PI3K_2205 and PI3K_5366 between verification populations and the third-generation sequencing population

3 討論

隨著PacBio 三代測(cè)序技術(shù)的發(fā)展, 測(cè)序準(zhǔn)確率日漸提高, 三代測(cè)序在基因組學(xué)和基因分型中的應(yīng)用也越來(lái)越廣泛。本研究針對(duì)水產(chǎn)動(dòng)物性狀相關(guān)標(biāo)記高通量篩選的需求, 建立了基于PacBio SMRT 靶向測(cè)序的候選基因關(guān)聯(lián)分析方法, 并在凡納濱對(duì)蝦抗弧菌性狀相關(guān)標(biāo)記的篩選中進(jìn)行了應(yīng)用, 證實(shí)了該方法是一種高通量、準(zhǔn)確可靠的方法, 本研究為水產(chǎn)動(dòng)物的性狀相關(guān)基因精細(xì)定位和分子育種研究提供了一種高效工具。

隨著高通量SNP 分型技術(shù)的發(fā)展, 利用全基因組關(guān)聯(lián)分析(GWAS)定位性狀關(guān)鍵基因成為了重要的技術(shù)手段, 在通過(guò)GWAS 獲得性狀相關(guān)QTL 后, 需要進(jìn)一步對(duì)QTL 進(jìn)行精細(xì)分析, 以鑒定與性狀緊密連鎖的分子標(biāo)記或者因果突變, 因此針對(duì)QTL 區(qū)域開(kāi)展候選基因關(guān)聯(lián)分析是精細(xì)定位的關(guān)鍵步驟。然而目前水產(chǎn)動(dòng)物中對(duì)性狀進(jìn)行精細(xì)定位的研究較少,一方面水產(chǎn)動(dòng)物GWAS 的研究起步較晚, 目前較多處于QTL 或者標(biāo)記定位階段; 另一方面缺乏針對(duì)靶向區(qū)域進(jìn)行高效關(guān)聯(lián)分析的方法。在凡納濱對(duì)蝦分子育種研究方面, 目前已鑒定出一些與生長(zhǎng)(Andriantahinaet al, 2012)、耐氨氮(Luet al, 2018)和WSSV 抗性(Liuet al, 2014)等性狀相關(guān)的QTL 或者SNP 標(biāo)記等, 并通過(guò)精細(xì)定位鑒定了PKC、MMD2、SRC 等生長(zhǎng)候選基因(Wanget al, 2019), 在這些研究中, 針對(duì)候選基因的關(guān)聯(lián)分析采用的是PCR 擴(kuò)增后一代測(cè)序或二代測(cè)序的方法。該方法需要針對(duì)目標(biāo)基因的外顯子區(qū)域分段設(shè)計(jì)引物, 分別擴(kuò)增后測(cè)序, 工作量大、費(fèi)用高, 對(duì)于序列比較長(zhǎng)的基因分析難度大。凡納濱對(duì)蝦基因組十分復(fù)雜, 是公認(rèn)的復(fù)雜基因組之一, 其重復(fù)序列約占基因組的23.93%, SSR 的平均長(zhǎng)度為72.21 bp, 是其他節(jié)肢動(dòng)物(20.11—31.91 bp)的兩倍以上(Zhanget al, 2019b)。重復(fù)序列在全基因組上出現(xiàn)頻率高以及長(zhǎng)度長(zhǎng)的特點(diǎn), 使得利用一代和二代測(cè)序技術(shù)在進(jìn)行對(duì)蝦基因分型面臨較大困難(Sulovariet al, 2019)。雖然對(duì)重復(fù)序列區(qū)域測(cè)序的難度較大, 但是這些區(qū)域通常顯示出較高的突變率,這些突變通常與疾病和進(jìn)化相關(guān)(Huddlestonet al,2017)。一代測(cè)序通常難以跨越連續(xù)的SSR 序列, 二代測(cè)序由于其讀長(zhǎng)限制, 需要對(duì)測(cè)序結(jié)果進(jìn)行拼接,而高重復(fù)序列又會(huì)影響拼接的準(zhǔn)確性。相較于一代和二代測(cè)序技術(shù)、三代測(cè)序兼具通量高、測(cè)序長(zhǎng)度長(zhǎng)的特點(diǎn), 通過(guò)長(zhǎng)片段靶向捕獲結(jié)合三代測(cè)序成為開(kāi)展凡納濱對(duì)蝦候選基因關(guān)聯(lián)分析的最佳選擇(Yueet al, 2017)。本研究建立了一套基于PacBio SMRT測(cè)序技術(shù)進(jìn)行候選基因關(guān)聯(lián)分析的實(shí)驗(yàn)流程和數(shù)據(jù)分析流程, 對(duì)長(zhǎng)片段靶向序列進(jìn)行全長(zhǎng)測(cè)序的同時(shí),也通過(guò)增加測(cè)序深度, 實(shí)現(xiàn)了對(duì)候選基因的精準(zhǔn)分型。三代測(cè)序數(shù)據(jù)的準(zhǔn)確性很大程度上受SMRT 建庫(kù)質(zhì)量的影響, 且擴(kuò)增子混合建庫(kù)本身也存在混合不均勻、擴(kuò)增子本身質(zhì)量不佳等問(wèn)題(DePristoet al,2011)。采用三代測(cè)序數(shù)據(jù)適配的分析軟件和流程,對(duì)于得到更精準(zhǔn)的基因分型結(jié)果至關(guān)重要。本研究采用IGV 等可視化工具, 對(duì)測(cè)序數(shù)據(jù)的連續(xù)性進(jìn)行了評(píng)估, 發(fā)現(xiàn)測(cè)序結(jié)果連續(xù)性較好, 在7 378 bp 的全長(zhǎng)中無(wú)測(cè)序?qū)е碌臄帱c(diǎn), 測(cè)序reads 連續(xù)性將大大減少后續(xù)單倍型分析的難度和工作量。

根據(jù)所建立的候選基因關(guān)聯(lián)分析方法, 結(jié)合團(tuán)隊(duì)前期發(fā)掘的與LvPI3K緊密連鎖的抗弧菌性狀相關(guān)分子標(biāo)記(Unigene19157_All__1806_348_223) (Zhanget al, 2019a), 本研究將LvPI3K基因作為凡納濱對(duì)蝦對(duì)抗弧菌性狀相關(guān)候選基因, 進(jìn)一步鑒定與抗弧菌性狀相關(guān)的突變位點(diǎn)。在候選基因LvPI3K共計(jì)7 378 bp 的序列內(nèi)發(fā)現(xiàn)91 個(gè)突變位點(diǎn), 平均每81 bp出現(xiàn)一個(gè)SNP, 其中有95.6%位于內(nèi)含子或者5'-UTR區(qū)域, 有4 個(gè)突變發(fā)生在外顯子區(qū)域, 由此推測(cè), 非編碼區(qū)中的突變發(fā)生頻率要高于編碼區(qū)中的發(fā)生率。內(nèi)含子區(qū)域的多態(tài)性對(duì)于性狀相關(guān)性的解釋較為困難, 目前主要被當(dāng)作與抗病功能基因緊密連鎖的分子標(biāo)記來(lái)應(yīng)用。本研究中與抗弧菌性狀呈極顯著相關(guān)的PI3K_5366_T/C 位點(diǎn)發(fā)生在PI3K 的第六個(gè)外顯子區(qū)域, 其T>C 突變會(huì)導(dǎo)致mRNA 的ACU 突變?yōu)锳CC,然而由于密碼子的簡(jiǎn)并性, 氨基酸序列并沒(méi)有變化。但是同義突變依舊具有巨大的研究?jī)r(jià)值, 其可以從核酸和蛋白兩個(gè)方面影響基因的表達(dá)。在核酸層面,同義突變通過(guò)調(diào)控剪接位點(diǎn)序列進(jìn)而產(chǎn)生多種剪接形式, 或者通過(guò)影響DNA 序列內(nèi)的GC 含量, 從而影響序列穩(wěn)定性等。在蛋白層面, 雖然同義突變并不改變其編碼的氨基酸序列, 但不同物種對(duì)于不同密碼子的偏好性不同, 從而影響mRNA 的翻譯效率。同義突變還會(huì)對(duì)蛋白質(zhì)的折疊產(chǎn)生影響, 從而改變蛋白自身的理化性質(zhì)。因此, 同義突變?cè)诜肿庸δ軐用娌⒉皇浅聊蛔? 可以通過(guò)參與基因表達(dá)調(diào)控的各個(gè)階段對(duì)基因功能產(chǎn)生影響(Plotkinet al, 2011)。

4 結(jié)論

本研究開(kāi)發(fā)了一種適用于凡納濱對(duì)蝦等水產(chǎn)生物的基于三代測(cè)序的高通量靶向測(cè)序基因分型技術(shù),并發(fā)掘了LvPI3K基因上91 個(gè)多態(tài)性位點(diǎn), 其中有21個(gè)位點(diǎn)與對(duì)蝦的抗弧菌性狀呈顯著相關(guān)(P<0.05), 有2 個(gè)位點(diǎn)(PI3K_2205_GC/G, PI3K_5366_T/C)與抗弧菌性狀呈極顯著相關(guān)(P<0.01), 這些位點(diǎn)的開(kāi)發(fā)對(duì)于對(duì)蝦抗弧菌性狀的分子育種有重要指導(dǎo)意義。