基于三代測序技術的凡納濱對蝦(Litopenaeus vannamei)高通量候選基因關聯分析方法的建立*

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

(1. 中國科學院海洋研究所實驗海洋生物學重點實驗室 青島 266071; 2. 中國科學院大學 北京 100049)

水產動物的經濟性狀大多是由微效多基因控制的數量性狀, 如生長速率、抗病性等(Yue, 2014;Gjedrem, 2015)。然而數量性狀易受環境等因素影響,導致了遺傳評估準確率偏低。抗性性狀由于表型度量困難、遺傳力低, 加之與生長或其他性狀之間存在負相關, 使得利用傳統育種技術對抗病性狀進行選育的效率較低(Cocket al, 2009)。隨著分子標記, 特別是以單核苷酸多態性(SNP)標記為主的基因分型技術的發展, 分子標記輔助選擇(MAS)和基因組選擇(GS)等分子育種技術在水產動物抗病育種中彰顯了巨大優勢, 比如個體評估的準確性高, 可在繁殖周期的早期進行選擇等。盡管分子育種技術優勢大, 但是在大多數水產動物的商業化育種中的應用仍然較少。在水產動物中開展高效分子育種面臨著兩大難點, 一者是缺少高密度、高精度的性狀連鎖標記(Meuwissenet al,2004), 二者是對標記和性狀之間關聯程度的評估有待優化(VanRadenet al, 2009)。因此鑒定與性狀緊密相關的基因或基因位點是水產動物分子育種亟需解決的問題。

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

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

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

本研究基于三代測序技術, 在凡納濱對蝦中建立了一種基于三代靶向測序的候選基因關聯分析技術, 并對前期篩選的抗弧菌性狀相關基因LvPI3K進行了標記發掘和關聯分析, 以篩選抗弧菌性狀相關標記。本研究不僅為水產動物提供了一種高效、低成本的候選基因關聯分析方法, 也為凡納濱對蝦的抗弧菌品種選育提供了有效的分子標記。

1 材料與方法

1.1 實驗材料

1.1.1 候選基因關聯分析材料的獲得 凡納濱對蝦遺傳材料來自渤海水產育種(海南)有限公司(海南,文昌), 分析用的群體為經多代選育的多個弧菌抗性家系和敏感家系的雜交F2代群體。

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

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

1.2 LvPI3K 基因組DNA 擴增

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

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

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

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

圖1 兩步PCR 將barcode 整合進擴增產物Fig.1 Incorporation of barcodes into the PCR amplicon via a two-step PCR approach

1.3 SMRT 文庫的構建及PacBio 測序

將192 個個體的第二輪PCR 擴增產物等摩爾混合后, 樣品送至天津生物芯片有限公司進行SMRT 文庫的構建及PacBio 測序, 文庫構建方式參考SMRT bell Express Template Prep Kit 2.0 的建庫方案。

1.4 數據分析

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

圖2 SMRT 測序數據分析流程圖Fig.2 Flowchart of the SMRT sequencing analysis

使用FastQC 與MultiQC 軟件對SMRT 測序數據進行整體的質量分析。使用SAMtools 及PLINK 軟件對FASTQ 文件進行位點測序深度與reads 長度的統計, 并使用R 軟件進行做圖。

1.5 對蝦抗弧菌性狀相關候選位點的驗證

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

2 結果

2.1 PacBio 測序數據分析

目的片段LvPI3K基因全長擴增子序列長度為7 378 bp (圖3a), 采用本實驗室構建的針對三代測序數據的專用分析流程, 對測序結果進行了全面分析。統計顯示有99.89%的位點測序深度大于20 000×, 位點平均測序深度為22 012×。對樣品的測序reads 總數進行統計, 結果顯示, 測序的192 個個體共產生了22 089 條reads, 其中有21 990 條(99.6%) reads 長度在目的片段長度周圍(7.4—7.6 kb) (圖3b)。通過FastQC與MultiQC 軟件對測序數據質量進行分析, 發現除極少數序列外, 三代測序數據均符合質控標準(Q30>90),并且全部樣品的平均測序深度為139×, 最高達902×,測序質量和測序深度均可以滿足后續變異發掘的要求(Edgeet al, 2019; Wenger et al, 2019)。

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

2.2 突變位點檢測及基因分型

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

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

2.3 候選基因關聯分析

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

圖5 LvPI3K 基因上與抗弧菌性狀呈極顯著相關(P<0.01)的標記及周邊標記的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 對蝦抗弧菌性狀相關位點(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

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

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

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

2.4 對蝦抗弧菌性狀相關候選基因位點的驗證

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

表4 PI3K_2205 和PI3K_5366 位點在驗證群體與三代測序分型群體中的基因分型情況對比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 三代測序技術的發展, 測序準確率日漸提高, 三代測序在基因組學和基因分型中的應用也越來越廣泛。本研究針對水產動物性狀相關標記高通量篩選的需求, 建立了基于PacBio SMRT 靶向測序的候選基因關聯分析方法, 并在凡納濱對蝦抗弧菌性狀相關標記的篩選中進行了應用, 證實了該方法是一種高通量、準確可靠的方法, 本研究為水產動物的性狀相關基因精細定位和分子育種研究提供了一種高效工具。

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

根據所建立的候選基因關聯分析方法, 結合團隊前期發掘的與LvPI3K緊密連鎖的抗弧菌性狀相關分子標記(Unigene19157_All__1806_348_223) (Zhanget al, 2019a), 本研究將LvPI3K基因作為凡納濱對蝦對抗弧菌性狀相關候選基因, 進一步鑒定與抗弧菌性狀相關的突變位點。在候選基因LvPI3K共計7 378 bp 的序列內發現91 個突變位點, 平均每81 bp出現一個SNP, 其中有95.6%位于內含子或者5'-UTR區域, 有4 個突變發生在外顯子區域, 由此推測, 非編碼區中的突變發生頻率要高于編碼區中的發生率。內含子區域的多態性對于性狀相關性的解釋較為困難, 目前主要被當作與抗病功能基因緊密連鎖的分子標記來應用。本研究中與抗弧菌性狀呈極顯著相關的PI3K_5366_T/C 位點發生在PI3K 的第六個外顯子區域, 其T>C 突變會導致mRNA 的ACU 突變為ACC,然而由于密碼子的簡并性, 氨基酸序列并沒有變化。但是同義突變依舊具有巨大的研究價值, 其可以從核酸和蛋白兩個方面影響基因的表達。在核酸層面,同義突變通過調控剪接位點序列進而產生多種剪接形式, 或者通過影響DNA 序列內的GC 含量, 從而影響序列穩定性等。在蛋白層面, 雖然同義突變并不改變其編碼的氨基酸序列, 但不同物種對于不同密碼子的偏好性不同, 從而影響mRNA 的翻譯效率。同義突變還會對蛋白質的折疊產生影響, 從而改變蛋白自身的理化性質。因此, 同義突變在分子功能層面并不是沉默突變, 可以通過參與基因表達調控的各個階段對基因功能產生影響(Plotkinet al, 2011)。

4 結論

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