基于微單體型分子標記的草魚親子鑒定方法

夏 雷 石米娟 張婉婷 段 攸 程瑩寅 吳 南 夏曉勤 ,

(1. 中國科學院水生生物研究所, 武漢 430072; 2. 中國科學院科學院大學, 北京 100149;3. 中國科學院種子創新研究院, 北京 100101)

魚類育種技術是水產養殖業可持續發展的重要保障, 目前已有大量的遺傳育種手段應用于魚類育種[1]。在各類經濟魚類的選育過程中, 所得到的子代群體往往數量龐大, 為了保證飼養環境條件的一致而將多個家系群體同池混養, 后期重建親本與子代之間的對應關系時, 需要使用親子鑒定技術[2]。

在通常情況下, 親子鑒定主要根據孟德爾遺傳定律, 以分子標記作為依據來判斷具體的子代與親代之間是否存在親子關系。目前親子鑒定主要采用微衛星(Simple Sequence Repeat, SSR)和單核苷酸多態性(Single Nucleotide Polymorphism, SNP)這兩種分子標記, 其中SNP在人類親子鑒定中應用較多[3],而SSR主要應用于水產養殖育種[2,4—6]。SSR具有信息含量高、多態性好的優點[7], 但標記本身的篩選過程比較繁瑣, 而且后續的親子鑒定實驗也較多地依賴于人力勞動, 樣本量較大時, 耗時長, 效率低;相比之下, SNP標記則具有不易發生變異和易于分析的優點, 但其多態性較低[8], 需要使用較多的SNP標記才可達到與SSR標記同等的效果。理論上SSR標記可以直接從測序數據與參考基因組比對后所得的插入缺失(INDEL)區域中進行篩選而獲得[9,10], 但直接獲得SSR的方法都受限于reads長度及堿基滑移方式, 對于重復次數相近或超過reads長度的SSR無法進行有效分型[10], 因此主要還是依靠傳統的實驗手段來檢測。單個SNP位點的分型也可以通過多種測序技術獲得, 2011年Davey等[11]對這些技術進行了總結, 包括重測序技術和RAD-seq(Restriction-site-associated DNA sequencing)等簡化基因組測序技術, 其中重測序技術可以獲得全基因組范圍內的SNP位點, 而簡化基因組測序技術只能獲得部分SNP位點, 但相較而言其成本較低。在水產動物的遺傳育種研究中, SSR分子標記主要用于親子鑒定[12—14]和群體遺傳多樣性分析[15]等工作, SNP分子標記則主要應用于遺傳連鎖圖譜構建[16]和GWAS分析[17]。

2013年Kidd等[18]結合上述兩種分子標記的優點, 提出了微單體型(Microhaplotype, MH)的概念,這是指長度在200 bp以內、可連鎖遺傳的SNP組合。作為一種分子標記, MH多態性高并且能夠穩定遺傳, 已應用于法醫行業[18—20]、人類群體分析[21,22]及祖先推斷[23,24]等。2018年, Garciafernandez等[25]以每個目標基因中的所有SNP位點作為一個單體型(Haplotype)分子標記, 即長片段的連鎖遺傳的SNP組合, 應用于鯛(Sparus aurata)的親子鑒定。另外微單體型分子標記也曾用于分析墨綠平鲉(Sebastes atrovirens)的親緣關系[26], 然而迄今沒有應用于魚類親子鑒定。

目前獲取微單體型或單體型分型的方法主要有兩種: 第一種方法是根據所獲得的各個SNP位點,用群體遺傳學的手段進行單體型推斷; 第二種方法是直接用個體基因組的測序數據的單體型組裝。用群體遺傳學推斷方法的常用軟件HaploView[27]、PHASE[28]、SHAPEIT[29]和Beagle[30]等, 它們對群體大小的依賴性較強, 對于大樣本群體可以獲得較準確的分型結果[31]。然而, 在魚類的遺傳育種中親本數量通常較小[4,15,32], 難以獲得適用于此方法的足夠數量樣品, 不容易做到準確分型?；谛蛄薪M裝的軟件主要有HapCUT2[33]、ReFHap[34]以及FLfinder[35]等, 目前僅用于單體型的推斷與分型, 而未應用于微單體型的分析。這種方法不受樣本數量限制, 僅對測序深度敏感, 測序深度不足時無法獲得準確的分型結果[36], 但隨著高通量測序成本的迅速下降, 測序深度已不再構成瓶頸問題。

本研究在基于個體測序數據的單體型組裝的基礎上, 開發了一套基于微單體型分子標記的親子鑒定流程, 其核心為微單體型標記的獲取、分型與篩選, 所需要樣本更少且準確率高, 適用于已有參考基因組的二倍體水產動物。為評估新方法的效果, 我們以一個草魚群體的全基因組重測序數據為例, 以其親子鑒定結果與SSR鑒定結果的一致性為主要指標對新方法加以驗證和評估。

1 材料與方法

1.1 基于個體序列組裝的微單體型獲取與親子鑒定流程

我們通過個體單體型組裝軟件獲取微單體型,整個親子鑒定流程主要分為三個步驟(圖1): 第一步是獲取親本微單體型。通過與參考基因組進行比對, 獲取每個親本個體的所有SNP位點, 依據測序結果使用HapCUT2軟件進行親本個體的單體型組裝, 從中選取微單體型區域, 將所有個體的區域集合起來, 獲取在該基因片段上的微單體型區域。按照劃分好的微單體型區域, 對每個個體進行微單體型組裝及分型。然后對子代進行微單體型的分型。對于上面從親本中得到所有微單體型區域, 為每個子代個體進行微單體型組裝。將該區域上子代與親本的全部分型結果形成列表, 依據各分型的索引構建每個個體在該微單體型下的分型結果。最后根據全部子代分型結果, 用CERVUS 3.0軟件[37]進行親子鑒定。

圖1 微單體型獲取和親子鑒定流程Fig. 1 Pipeline for microhaplotype genotyping and paternity testing

1.2 微單體型區域(Microhaplotype region, MR)的劃分

用HapCUT2軟件處理親本的重測序數據, 獲取每個親本個體的MR, 即一個單體型的首尾SNP位點在參考基因組上的位置。集中所有親本個體的MR, 合并任何在位置上有重疊的MR, 構成一個序列較長的、信息量更高的MR總庫。

接下來將總庫中的每一個MR逐一與各親本的MR進行比較。由于總庫MR是合并所有親本中相重疊MR所得, 那么, 每一個總庫MR必然不短于單個親本中與之有重疊的MR, 如果其非重疊區域(超出的部分)在該親本序列中也存在雜合的SNP位點,那么這個MR將按重疊與非重疊區域拆分成兩個區域, 僅保留至少有2個SNP位點的區域作為新的MR,這樣可增加后續分析結果的準確性。最終所得的MR總庫將用于后續分型研究。

1.3 微單體型分型及分析方法

依據MR總庫可對各親本進行微單體型的分型, 即針對總庫中的每個MR確定各親本的具體微單體型序列。二倍體親本在每個MR最多只可能有2種微單體型, 這可以根據HapCUT2的分型結果篩選, 也可以參考該親本與參考基因組比對所得的VCF文件來確定。如果由于測序錯誤或數據量不足, 無法明確地鑒定某MR的1—2種微單體型序列,則舍棄該MR。剩下的各有效MR中所得到微單體型序列都可僅保留各SNP位點的堿基, 然后計算各MR在親本群體中的信息熵, 公式如下:

式中,n表示該MR中微單體型序列的種類數,pi表示第i種微單體型在親本群體中出現的頻率,H表示該MR在親本群體中的信息熵。

使用與親本微單體型分型的方法, 也可以在每個子代中獲得各MR的微單體型序列, 用于后續的親子鑒定工作。

1.4 基于群體遺傳學方法的單體型分子標記獲取與分析

傳統上單體型標記是通過群體遺傳學分析推斷而獲得的, 因此本研究也使用此方法進行了單體型分子標記的分型。分析過程比較簡單: 首先通過SHAPEIT軟件將每個基因構建一個單體型, 再依據各分型出現順序標注每個個體的分型情況。

1.5 使用SSR標記進行親子鑒定的方法

我們從草魚參考基因組[36]中篩選部分短序列重復片段作為候選SSR標記, 并設計引物對2尾父本(1#和M10)和3尾母本(2#, F13和F40)進行擴增,擴增程序為95℃預變性10min, 95℃變性45s, 55℃退火45s, 72℃延伸1min, 循環30次, 72℃延伸10min,4℃保存。使用ABI 3730對各標記的擴增產物進行毛細管凝膠電泳, 以條帶大小為分型依據判斷其多態性。最終選取5個多態性高的SSR標記進行親子鑒定, SSR標記的相關信息見表1。171尾子代個體也同樣進行目標SSR標記擴增及分型, 然后進行親子鑒定分析。

表1 用于親子鑒定的SSR標記信息Tab. 1 Information of the SSRs used for the paternity test

1.6 組裝微單體型標記使用的數據及其處理方法

我們以草魚基因組重測序數據(尚未發表)進行了微單體型組裝, 數據來自上述五尾親本及其171尾子代個體, 其中親本測序深度為30X, 子代測序深度為15X, 并且評估了微單體型標記應用于親子鑒定的結果與SSR鑒定結果的一致性。重測序數據經過了質量分析(FastQC)、與參考基因組的比對(BWA[38])等步驟, 最后利用SAMtools[39]和GATK 4.0[40]軟件獲取各樣本的SNP位點, 作為SNP信息集。

考慮到子代個體數并不多, 我們僅選擇了15個富含SNP的基因, 即adamts20(ADAM metallopeptidase with thrombospondin type 1 motif 20)、brca2(breast cancer 2)、dlc1(deleted in liver cancer 1)、gbp(guanylate binding protein)、lgals9(galectin 9)、lrp5(LDL receptor related protein 5)、meis2b(meis homeobox 2)、mrps23(mitochondrial ribosomal protein S23)、msi2(musashi RNA binding protein 2)、nos2b(nitric oxide synthase 2)、prtga(protogenin A)、rpz4(rapunzel 4)、snx14(sorting nexin 14)、thsd4(thrombospondin type 1 domain containing 4)和zmym4(zinc finger MYM-type containing 4), 以代表整個基因組來評估我們方法的效果, 這樣就大大減少了計算量。從這15個基因中SNP的密度分布可以看出, 在這些基因中確實存在一些SNP較為集中的片段(圖2)。微單體型由連鎖遺傳的SNP位點組合構成且這些SNP位點通常距離較近, 因此這些分布較為集中的SNP可以組合成潛在的微單體型。

1.7 親子鑒定標記組合(Paternity Test Marker Combinations, PTMC)篩選

為了獲得PTMC, 我們首先依據親本性別構建全部可能的親本對, 為降低連鎖遺傳帶來的影響,隨后依次挑選若干來源于不同基因的標記, 依據親本的分型情況, 可得到每一個親本對的后代在這些標記中的所有可能的分型組合。若存在某兩個親本對的子代分型有重合的情況, 則認為所選的標記組合無法進行親子鑒定; 反之, 若所有親本對的子代分型均唯一, 則認為所選標記組合可以用于親子鑒定, 并作為可用的親子鑒定標記組合(PTMC)輸出。通過此步驟可以得到所有可能的PTMC, 從中選用平均信息熵最高的一套PTMC用于親子鑒定可以提高結果的準確性。

基因組重組可能在子代中產生異于親本的分型類型, 從而妨礙親子鑒定。適當增加冗余標記可以應對這種分型矛盾。冗余標記的選擇可以采用一種簡單的策略, 即從與所有已選用標記位于不同基因(或染色體)上的標記中, 選擇信息熵最高的一個。重復該過程則可以選擇多個冗余標記。

2 結果

2.1 基于個體序列組裝的微單體型獲取

我們親本重測序數據對15個基因區域進行了微單體型的組裝, 統計了每個基因上微單體型區域和SNP位點的數量, 以及第一個SNP位點到最后一個SNP位點之間的距離(表2), 并計算信息熵作為評價微單體型區域信息量的標準。總體而言, 大部分微單體型區域的信息熵的值介于1和2之間(圖3A),但其分布范圍在各基因之中有較大的變化(圖3B),并且信息熵值也隨SNP數量的增長而呈上升趨勢(圖3C)。

2.2 用單個基因的MR進行的親子鑒定

親子鑒定使用的軟件為CERVUS 3.0[37], 參數為默認值。

圖2 15個基因中SNP分布密度Fig. 2 Distribution of SNPs in 15 genes

由于親本的測序深度高于子代, 有些從親本上鑒定的MR序列在子代中沒有被覆蓋, 只有在所有子代中出現的MR才是有效MR。我們首先嘗試將一個基因上的全部有效MR作為一套標記進行親子鑒定, 其結果與用SSR鑒定結果的一致性見表2。不同基因的鑒定能力差異很大, 其中adamts20、brca2、prtga和snx14鑒定的結果與SSR分析的一致性較高(＞95%), 而gbp、lrp5、meis2b和mrps23的一致性較低(＜50%)。親子鑒定的效果與基因中微單體型標記的數量或信息熵之和存在一定的關系, 但并非簡單的線性關系(R2=0.2302)。對于我們所測試的數據, 當基因中標記數或信息熵之和低于10時,親子鑒定的效果與之有正相關性, 但之后即達到飽和, 與SSR分析的一致性維持很高(圖4)。

2.3 MR組合的選擇與親子鑒定

如果不考慮基因重組, 根據5個親本所有標記的實際分型情況, 理論上僅需3個標記就能分辨出它們之間任何子代的親本, 而這樣的3標記組合總共有189個。從中選擇平均信息熵最高的一組作為標記進行親子鑒定, 結果在171個子代中, 與SSR鑒定結果不同的個體僅5個, 一致率達到97.08%。再逐次增加第1個和第2個冗余標記后, 不一致個體減少為2個和1個, 一致的比率進一步提高到98.83%和99.42%(表3)。

2.4 基于群體遺傳推斷的單體型分型結果

除了使用單個個體的測序數據進行微單體型組裝的方法外, 還可以利用群體遺傳學的方法進行單體型推斷, 從而獲取微單體型分型。為了比較這兩種方法在親子鑒定上的區別, 我們使用SHAPEIT軟件獲取了上述176尾個體在15個基因上的單體型分型情況。通過比較發現, 當某個基因中SNP數量過多時, 個體分型種類十分復雜, 不能用于親子鑒定, 因此我們選擇SNP數量較少的基因msi2嘗試進行單體型分型。在5尾二倍體親本中, 該基因的SNP總共出現了10種分型, 也就是說, 在5個親本中出現了10個不同的等位基因類型。在正常情況下,如果不考慮重組以及各種偏差, 任何兩個親本得到的后代都應該是這10個等位基因的組合, 可以被明確地鑒定其等位基因來自哪兩個親本。然而, 通過群體遺傳學推斷, 在171尾子代中共檢出42種分型,出現了34種親本并不存在的分型, 并且有一個母本(2#)的兩種分型還不包含在其中。

3 討論

3.1 標記獲取方法

目前標記的獲取主要有三種策略: (1)從已經報道的多態性標記位點中獲得[32,41]; (2)選取適當數量(通常小于100尾)的個體數據用于測試標記多態性[12,42,43]; (3)使用全部親本子代組成的群體進行標記篩選[25,44]。本研究所采用的方法為第二種, 并且考慮到重組和變異的頻率不可能很高, 子代的絕大部分標記均與其親本一致, 因此我們用于獲取標記的個體為全體親本, 不包括子代個體。在魚類育種研究中, 通常是很有限的親本產生大量的子代。在親子鑒定過程中, 我們只需要通過數量有限的親本獲取少量標記, 即可從大量的子代中獲取出必要的信息, 這種做法具有快速和低成本的優點, 非常有利于大規模的親子鑒定。

表2 各基因上的SNP位點與微單體型區域的數量以及親子鑒定準確率Tab. 2 SNP loci and MR in each gene and the accuracy of paternity testing

3.2 單體型/微單體型標記的分型方法比較

本研究比較了基于群體遺傳推斷的方法和基于個體序列組裝的方法。前者使用的算法大多為EM算法或其他概率算法, 核心為參數優化, 從理論本身而言, 需要大樣本量來進行推斷, 對于樣本量較少而SNP位點較多的群體分型結果不準確[31]。依據孟德爾分離定律與自由組合定律, 理論上絕大多數子代的分型應與其對應的親本相一致, 然而在本研究中該方法的分型結果與親本的實際分型情況相差很大, 因此我們沒有使用這些分型數據進行親子鑒定。出現這樣親子分型差異較大的情況, 說明對于這種方法來說, 包含171尾子代個體的群體樣本量依然太小, 導致進行概率計算時無法獲得準確的估計, 從而產生了錯誤的分型結果。可見使用群體遺傳學方法進行單體型分型對于樣本數據量有很高的要求, 在應用中有明顯的局限性。

圖3 基于個體序列組裝的微單體型分析結果Fig. 3 Analysis of the microhaplotypes assembled using the sequence data of individuals

圖4 親子鑒定一致性與MH標記數量(A)和各基因有效標記信息熵之和(B)的關系Fig. 4 Consistency between the SSRs used in the paternity test and MH marker numbers (A) or the sum of the informative indexes of the markers within each gene (B)

表3 所選用的MR標記與171尾子代個體的親子鑒定Tab. 3 MR markers adopted and paternity test results of 171 offspring

相比之下, 利用測序數據直接進行單體型或微單體型組裝的方法不僅更為準確地從親本中得到了各個標記的多種分型, 也能從子代中找到這些分型, 從而順利完成親子鑒定過程。這種方法的主要劣勢在于易受到單個位點測序深度和測序錯誤的影響[36]。在本研究中, 子代測序深度為15X, 在某些SNP位點上測序深度不足, 導致無法獲取某些微單體型的分型。不過, 在實際應用中這也不會成為一個問題, 因為在經通過親本的測序數據選定標記之后, 只需要在子代中擴增這些微單體型標記, 直接測序即可得到其準確的分型, 而不需要進行高成本的基因組重測序與序列組裝。此外, 目前靶向捕獲基因組聯合二代測序技術已經得到了廣泛應用[45],同樣可以用于在大量的子代群體中測序少量的標記序列。很顯然, 在高通量測序技術已經普及的背景下, 基于序列組裝的方法更加適合于魚類親子鑒定的工作。

3.3 MH在親子鑒定中的應用

本研究比較了用MH與SSR進行親子鑒定的結果, 并未得到完全的一致, 可能有兩方面的原因。首先, 某些個體測序深度不足, 導致無法分型, 這已在上文進行討論。其次, 各標記多態性不同, 導致親本區分度不同。從圖2和表2來看, 一致率低的基因往往只有極少或根本不存在SNP密集分布區域, 導致這些基因的絕大多數MH分型種類也較少,大量親本為純合子, 對各親本的區分度不足, 而目前常用的CERVUS 3.0軟件傾向于選擇純合子親本作為親子鑒定的輸出結果[46], 因此, 對于這類基因,MH親子鑒定得到的親本往往是純合子, 而與真實親本對可能有較大誤差, 最終影響親子鑒定結果的準確性。相反, 一致率高的基因大多存在若干SNP分布密集的區域且有效MR數較多, 即使單個基因近似連鎖遺傳, SNP密集的MH具有高多態性, 使得此基因MH標記總信息熵高, 可彌補其他MH區分度低的缺陷, 從而保留較高親本區分度。此外, 在這些MH上親本大多為雜合子, 因此可以極大提升親子鑒定效果, 最終親子鑒定一致率較高(圖4B)。最后, MH標記內部也可能存在的重組問題, 雖然這種可能性很小, 但一旦發生, 子代分型結果便無法與其真實親本對完全對應, 從而導致親子鑒定結果出現誤差, 這種情況在SSR標記中也存在。為盡可能減少上述因素的影響, 本研究使用了高多態的分子標記及冗余標記, 其中高多態的分子標記用于提升標記區分度, 冗余標記用于矯正測序深度和重組帶來的分型誤差, 最終親子鑒定一致率得到了提升(表3)。

MH標記應用于親子鑒定的優勢在于此類標記易于獲取、分型及篩選, 如果結合靶向測序就能夠以極低的成本高效而準確地完成子代的分型與親子鑒定; 其劣勢則在于鑒定能力嚴重依賴于MH標記的多態性, 而且使用全基因組重測序數據分析時,有些標記可能會受局部測序偏低的影響, 此外, 標記內重組也影響親子鑒定效果。因此, 我們的方法首先根據親本序列對目標區域的SNP密度進行評估, 篩選出高多態性的MH分子標記, 從而提升所用標記的分辨能力; 加入冗余標記則可以盡可能地減少測序深度低和標記內重組帶來的影響; 如果對目標基因或片段進行靶向測序, 就能夠完全排除局部測序深度低的干擾。

3.4 前景展望

雖然本研究的范例中僅從15個基因序列里篩選了親子鑒定的標記, 但對于有參考基因組的物種,該方法完全可以在全基因組范圍內進行微單體型標記的獲取與分型。通過計算這些微單體型信息熵值以篩選信息熵最高的微單體型標記, 提升標記的多態性和分辨能力, 可以用于個體鑒定、親子鑒定, 以及更廣泛的親緣關系鑒定, 甚至還可將其應用于經濟性狀的關聯分析中, 作為全基因組選擇育種的分子標記。