基于家系基因組測序?qū)︱吇蚪M結(jié)構(gòu)多樣性分析

任秀娟 蘇少鋒 李雅靜 藺雅楠 賈紫潔 丁文淇 白東義 李 蓓 杜 明 芒 來 趙一萍*

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 動物科學(xué)學(xué)院/馬屬動物研究中心,呼和浩特 010018;2.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)科學(xué)院,呼和浩特 010031)

馬(Equuscaballus)和驢(Equusasinus)屬于2個獨立的物種,其物種分離發(fā)生在約400～450萬年前[1]。馬屬動物的物種形成伴隨著快速的染色體重排,其重排速率為2.9～22.2次/百萬代[2-3]。馬(2n=64)和驢(2n=62)已經(jīng)進化形成了兩套結(jié)構(gòu)和功能完善的獨立基因組[4-5]。生理特征比較表明,驢具有更有效的能量代謝和更強的免疫力,而馬的反應(yīng)更靈敏和運動能力更強[6]。馬和驢雜交能夠產(chǎn)生后代(馬騾,2n=63,公驢和母馬的雜交后代;驢騾,2n=63,公馬和母驢的雜交后代)。盡管有母騾產(chǎn)駒的報道,但騾基本上不能自然交配產(chǎn)生后代。

根據(jù)Dobzhansky-Muller不相容定律(Dobzhansky-Muller incompatibility),馬和驢物種間進化分離的2個或多個等位基因之間不能有效的相互作用,可能導(dǎo)致騾的適應(yīng)性存在缺陷,表現(xiàn)為不育和溶血等[7-8]。例如,物種形成基因的反復(fù)適應(yīng)性進化[9]、細(xì)胞質(zhì)基因和核基因的沖突進化[10]以及X染色體的減速分裂驅(qū)動等[11],均可導(dǎo)致異種雜交不相容。劇烈的染色體結(jié)構(gòu)變異也可能是馬和驢雜交不相容的遺傳因素。例如,染色體易位可導(dǎo)致部分異種雜交個體的某些基因完全丟失,引起不相容[12]。Prdm9基因的拷貝數(shù)差異,導(dǎo)致小鼠(Musm.musculus×Musm.domesticus)品種間雄性雜交個體減數(shù)分裂重組失敗[13]。馬和驢屬于古老的物種,除了物種形成基因外,在長期的進化過程中,其基因組之間累積了更復(fù)雜的遺傳不相容。例如,受平衡選擇形成種間多態(tài)性的免疫相關(guān)基因,是許多物種間雜交障礙的主要原因。對小鼠[14]和硬骨魚[15]的研究表明,主要組織相容性復(fù)合物(Major histocompatibility complex,MHC)基因的種間多態(tài)性會降低F1雜交個體的適應(yīng)性。

騾如何協(xié)調(diào)雙親基因組的不平衡并保證其自身生存的分子機制仍然未知。大量研究表明,在胚胎發(fā)育早期,會發(fā)生合子后體細(xì)胞突變,同時異常的染色體結(jié)構(gòu)變異會觸發(fā)機體的自我糾錯過程。這些突變大部分會導(dǎo)致機體發(fā)生癌癥和一些罕見的發(fā)育障礙疾病[16-17]。然而,大部分存活的雜交個體伴隨著基因組的高頻突變。例如,對植物的研究表明,高度雜合個體的全基因組突變率高于純合個體[18-19]。異種雜交可能會激活大量轉(zhuǎn)座子的活性,從而促使更多結(jié)構(gòu)變異的形成。例如,金魚(Carassiusauratus)和鯉魚(Cyprinuscarpio)雜交,其F1雜交個體基因組發(fā)生了高頻率的點突變和大片段的結(jié)構(gòu)變異[20]。王艷欣[21]對馬騾和驢騾的膀胱等10個組織轉(zhuǎn)錄組的研究發(fā)現(xiàn),騾不同組織均發(fā)生了高頻率的結(jié)構(gòu)變異。

異種雜交個體發(fā)生的高頻突變,可能反映了物種間進化累積的差異突變導(dǎo)致的遺傳不相容。但是,雜交應(yīng)激觸發(fā)的快速突變,是否能緩沖其親本單倍體基因組之間的不平衡,從而提高雜交個體的適應(yīng)性尚不清晰。本研究基于馬屬動物三成員家系的全基因組Illumina測序數(shù)據(jù),旨在分析馬和驢雜交應(yīng)激觸發(fā)的騾基因組結(jié)構(gòu)多樣性,揭示馬和驢遺傳不相容以及影響騾適應(yīng)性的可能基因座,為進一步開展馬和驢異種雜交的遺傳基礎(chǔ)及其分子機制的研究提供候選遺傳位點。

1 材料與方法

1.1 實驗動物和測序數(shù)據(jù)

本研究的實驗動物包括1匹雌性蒙古馬、1匹雄性家驢和1匹雌性馬騾。3匹實驗動物來自于1個馬屬動物異種雜交家系。基因組序列由內(nèi)蒙古農(nóng)業(yè)大學(xué)馬屬動物研究中心產(chǎn)生和存儲。所有數(shù)據(jù)已提交NCBI,驢的項目登錄號:PRJNA205517[6],馬和騾的項目登錄號:PRJNA842856(SRA登錄號:SRR19427107和SRR19427108)[22]。

1.2 驢、馬和騾參考基因組的序列比對

以純血馬基因組(Equcab3.0)[5]和驢基因組(ASM1607732v2)[23]為參考,使用BWA(version 0.7.5a-r416)軟件的默認(rèn)參數(shù),將驢、馬和騾的高質(zhì)量reads比對到參考基因組[24]。使用SAMtools(version 0.1.19-44428cd)軟件獲得唯一比對,參數(shù)為“-q 30”[25]。使用Picard(version 1.93,http:∥sourceforge.net/projects/picard/)軟件標(biāo)記測序過程產(chǎn)生的潛在重復(fù)。

1.3 驢、馬和騾SNP和InDel的識別

使用SAMtools和GATK(version 3.5-0-g36282e4)軟件[26]對3個樣本分別識別SNP和InDel。2個軟件分析結(jié)果的交集用于后續(xù)分析。使用GATK的HaplotyCaller模塊識別SNP,參數(shù)為-standa_call_conf 30,-standa_emit_conf 10。采用VariantFiltration命令,根據(jù)(http:∥www.broadinstitute.org/gatk/guide/best-practices)推薦參數(shù)對獲得的SNP進行嚴(yán)格過濾。為提高SNP的陽性率,在上述嚴(yán)格過濾的基礎(chǔ)上,進一步過濾掉符合以下4個標(biāo)準(zhǔn)的SNP:(i)位于低復(fù)雜度或簡單重復(fù)區(qū)域;(ii)測序深度<4或>50;(iii)InDel周圍50 bp范圍內(nèi);(iv)Gap周圍10 bp范圍內(nèi),至少3個SNPs被識別。

1.4 騾de novo SNP的識別

騾denovoSNP是指針對參考基因組的一個堿基位點,騾基因組的該位點至少有1個等位基因不同于親本。騾denovoSNP基因型組合見表1。采用兩種方法分別識別denovoSNP,并取交集用于后續(xù)分析。首先,參考Roach的方法[27],利用bcftools軟件的query工具和vcftools軟件,按表1的基因型組合,識別騾denovoSNP。其次,使用VarScan軟件的“trio”命令,識別騾denovoSNP,參數(shù):--min-coverage 10,--min-var-freq 0.20,--p-value 0.05,--adj-var-freq 0.05,--adj-p-value 0.15[28]。將兩種方法獲得的denovoSNP取交集。本研究三成員家系的試驗設(shè)計,因為基因型的錯誤分配,可能會影響騾denovoSNP的準(zhǔn)確識別[29]。因此,將定位在拷貝數(shù)變異(CNV)區(qū)域和重復(fù)序列區(qū)域的SNP屏蔽掉。用CNVnator軟件[30]檢測CNV和用RepeatMasker軟件[31]檢測重復(fù)序列。用于后續(xù)分析的denovoSNP符合以下2個標(biāo)準(zhǔn):1)SNP位點的reads支持?jǐn)?shù)≥10;2)SNP等位位點的reads支持?jǐn)?shù)≥5[32]。

表1 De novo SNPs的基因型組合

ANNOVAR軟件用于denovoSNP的功能注釋。對于定位到基因間區(qū)的SNP,僅保留上下游5 kb范圍內(nèi)的基因。

1.5 驢、馬和騾CNV檢測及注釋

使用CNVnator軟件及其推薦參數(shù)識別CNV[30]。使用以下標(biāo)準(zhǔn)過濾原始CNV:獲得性CNV(duplication)的RD(normalized read depth)>2,缺失性CNV(deletion)的RD<0.4,P≤0.05,q0≤0.5,Length>1 000 bp。將驢、馬和騾位置重疊≥1 bp的CNV,合并為1個拷貝數(shù)變異區(qū)域(CNVR)。

使用bedtools軟件對CNV進行基因注釋,標(biāo)準(zhǔn)是CNV的bed文件和基因組的gtf文件至少有1個堿基的重疊[33]。

1.6 基因富集分析

使用R軟件的clusterProfiler包,進行KEGG富集分析,純血馬參考基因組的參數(shù)設(shè)置為:organic=“ecb”,keyType=“KEGG”;驢參考基因組的參數(shù)設(shè)置為:organic=“eai”,keyType=“KEGG”。

2 結(jié)果與分析

2.1 驢、馬和騾測序數(shù)據(jù)及比對

經(jīng)嚴(yán)格質(zhì)控后存儲于內(nèi)蒙古農(nóng)業(yè)大學(xué)馬屬動物研究中心的驢、馬和騾的高質(zhì)量基因組數(shù)據(jù),分別為100.01、103.78和114.36 Gb[6,22]。將基因組數(shù)據(jù)分別比對到純血馬參考基因組(Equcab3.0)和驢參考基因組(ASM1607732v2)。以純血馬基因組為參考和以驢基因組為參考時,至少98%的基因組覆蓋深度大于1×。從表2可知,分別以馬和驢的基因組為參考,讀段的比對率不存在比對的偏向性。

表2 參考基因組比對結(jié)果

2.2 驢、馬和騾SNP和InDel檢測分析

2.2.1驢、馬和騾InDel和SNP識別

以純血馬基因組為參考,驢、馬和騾識別的高質(zhì)量InDels分別為2 110 786、402 533和1 816 910。以驢基因組為參考,驢、馬和騾識別的高質(zhì)量InDels分別為527 351、2 279 875和2 036 942。以純血馬基因組為參考時,驢和騾的InDel數(shù)明顯高于馬,這符合蒙古馬和純血馬的親緣關(guān)系更近,騾是驢和馬的雜交個體。同樣的原因,當(dāng)以驢基因組為參考時,馬和騾的InDel數(shù)明顯高于驢。由圖1可知,InDel基本均勻地分布于常染色體。

(a)、(b)和(c)表示以驢基因組為參考時InDel的常染色體密度分布。(d)、(e)和(f)表示以純血馬基因組為參考時InDel的常染色體密度分布。(a)和(d)表示驢;(b)和(e)表示馬;(c)和(f)表示騾。

以驢基因組為參考,分別識別驢、馬和騾的高質(zhì)量SNPs為3 212 499、23 549 224和21 870 390。馬(0.968 3%)和騾(0.899 3%)的SNP頻率明顯高于驢(0.132 1%)。實驗室已發(fā)表數(shù)據(jù)表明,以純血馬基因組為參考,分別識別驢、馬和騾的高質(zhì)量SNPs為23 819 055、5 012 403和23 426 241;蒙古馬共識別5 012 403個SNPs,雜合SNP數(shù)高于純合,驢(0.950 1%)和騾(0.934 4%)的SNP頻率明顯高于馬(0.199 9%)[22]。以上結(jié)果表明,由于驢、馬和騾與參考基因組物種的親緣關(guān)系不同,導(dǎo)致與參考基因組不同物種識別的SNPs頻率遠高于相同物種(表3)。以純血馬基因組為參考識別的騾SNPs數(shù)略高于以驢基因組為參考的SNPs識別結(jié)果。由圖2可知,SNP基本均勻分布于常染色體。

(a)、(b)和(c)表示以驢基因組為參考時SNP的常染色體密度分布。(d)、(e)和(f)表示以純血馬基因組為參考時,SNP的染色體密度分布[22]。(a)和(d)表示驢;(b)和(e)表示馬;(c)和(f)表示騾。

表3 SNPs的統(tǒng)計

2.2.2騾denovoSNP

騾denovoSNP是指針對參考基因組的1個堿基位點,騾的該位點至少有1個等位基因不同于親本。對上述識別的高質(zhì)量SNP,進一步屏蔽掉重復(fù)序列和CNV區(qū)域的SNP。以純血馬基因組和驢基因組分別為參考,與親本相比,騾分別識別了555和419個denovoSNPs。與SNP識別趨勢一致,騾denovoSNP的識別可能存在馬參考基因組的偏向性。本研究的三成員家系中,騾denovoSNP的頻率是1.72×10-7～2.21×10-7(DenovoSNPs/全基因組序列)。

2.2.3騾denovoSNP相關(guān)基因注釋

為研究騾denovoSNP的潛在功能,使用ANNOVAR軟件,對SNP上下游5 kb范圍內(nèi)的基因進行注釋。以純血馬基因組為參考,555個denovoSNPs共注釋658個相關(guān)基因。以驢基因組為參考,419個denovoSNPs共注釋540個相關(guān)基因。DenovoSNP均最高頻率的注釋于基因間區(qū)(表4)。以兩套基因組為參考注釋denovoSNP,均發(fā)現(xiàn)大部分相關(guān)基因直接參與機體的免疫功能。例如,純血馬參考基因組注釋的基因,包括MHC I類基因(如MHCX1和LOC100049798)和MHC II類基因(如DQA和DQB)、Src家族酪氨酸激酶(如FYN)、表皮生長因子受體基因(EGFR)以及癌癥相關(guān)基因(如APC和PTEN、HRAS)[22]。驢參考基因組注釋的基因,包括MHC I類基因(如LOC106825028和LOC106843523)和MHC II類基因(如LOC106830318和LOC106834623)、Src家族酪氨酸激酶(如FYN)、表皮生長因子受體基因(EGFR)以及腫瘤抑制基因(如APC和PTEN)。

表4 De novo SNP的分布統(tǒng)計

2.2.4騾denovoSNP相關(guān)基因富集分析

對騾denovoSNP相關(guān)基因進行KEGG富集分析,發(fā)現(xiàn)大部分顯著富集(P<0.05)的通路參與機體的免疫應(yīng)答過程。以純血馬基因組和驢基因組為參考時,分別有43和27條通路被顯著富集,其中17條通路是相同的(圖3)。這些通路主要與機體的免疫過程相關(guān),包括在適應(yīng)性免疫應(yīng)答中起關(guān)鍵作用的抗原加工和提呈過程;機體排異相關(guān)過程:移植物抗宿主病、同種異體移植排斥;以及自身免疫性疾病通路:I型糖尿病、類風(fēng)濕性關(guān)節(jié)炎。

圖3 De novo SNP相關(guān)基因的KEGG富集

2.3 驢、馬和騾CNV特征分析

2.3.1驢、馬和騾CNV的檢測

以純血馬基因組為參考,驢、馬和騾分別獲得3 761、2 126和2 242個 CNVs(表5)。將CNVs合并后,共獲得5 178個CNVRs(表6)。騾特異性CNVRs總長度為2.15 Mb,平均長度為5.43 kb。馬和騾CNVRs總長度為7.64 Mb,平均長度為11.35 kb。驢和騾CNVRs總長度為5.41 Mb,平均長度為16.90 kb。驢、馬和騾CNVRs總長度為13.05 Mb,平均長度為26.74 kb。

表5 CNVs的統(tǒng)計

以驢基因組為參考,驢、馬和騾分別獲得3 572、7 812和6 020個CNVs(表5)。將CNVs合并后,共獲得8 967個CNVRs(表6)。騾特異性CNVRs總長度為3.77 Mb,平均長度為4.50 kb。馬和騾CNVRs總長度為15.14 Mb,平均長度為7.93 kb。驢和騾CNVRs總長度為3.00 Mb,平均長度為7.79 kb。驢、馬和騾CNVRs總長度為17.09 Mb,平均長度為11.39 kb。

2.3.2騾特異性CNVR相關(guān)基因功能分析

騾特異性CNVR是指相對參考基因組的某片段,馬和驢無拷貝數(shù)變異且只有騾發(fā)生的CNVR。以純血馬基因組為參考時,396個騾特異性CNVRs,共注釋226個基因。以驢基因組為參考,839個騾特異性CNVRs,共注釋435個基因。

KEGG富集分析顯示,以純血馬基因組為參考和以驢基因組為參考分別富集了93和145條pathways,其中共同富集的通路有66條。顯著富集(P<0.05)的通路分別是42和13條。盡管一些通路的富集沒達到顯著水平,但是對于分析突變基因的相關(guān)功能仍然有意義。共同富集的66條pathways主要與機體的免疫過程相關(guān)(圖4),包括抗原加工和提呈(如EQMHCB2和LOC100049798)、MAPK信號通路(如MAP3K2)、JAK-STAT信號通路(如AOX1)、NF-κB信號通路(如CARD14和TRIM25)、cGMP-PKG信號通路(如PDE5A和SLC8A2)。癌癥相關(guān)過程,包括癌癥途徑(如APC2和SOS2)、Ras信號通路(如RAP1B)、Wnt信號通路(如NKD2)和癌癥中轉(zhuǎn)錄調(diào)控異常等。

Donkey表示以驢基因組為參考時的富集結(jié)果;Horse表示以馬基因組為參考時的富集結(jié)果。

3 討 論

本團隊已有研究表明,以純血馬基因組為參考,使用標(biāo)準(zhǔn)流程識別的驢和蒙古馬SNPs及其頻率與已有報道相符[22,34-36]。使用相同的方法,以驢基因組為參考,識別了驢2 153 364個雜合SNPs,該數(shù)值和已發(fā)表的研究數(shù)據(jù)(2 187 070,以驢Maral har基因組為參考)相近[6]。驢雜合SNP頻率(0.088 5%)略高于Wang等[23]的報道(0.077 98%)。以上結(jié)果說明,本研究識別的SNPs數(shù)據(jù)可信。以驢基因組為參考識別的騾SNPs和denovoSNPs均略低于以純血馬基因組為參考的結(jié)果,說明SNP的識別可能存在純血馬參考基因組的偏向性。

結(jié)合以純血馬基因組為參考識別騾denovoSNPs(555)的結(jié)果可知[22],騾denovoSNPs數(shù)遠高于用相同方法識別的純種黑猩猩的denovoSNPs數(shù)(45)[37]。騾denovoSNP的突變率(1.72×10-7～2.21×10-7)高于人類(0.82×10-8～1.70×10-8)[38-39]和馬(7.24×10-9)[1]的自然突變率。騾基因組更高的突變頻率,符合Duncan在1915年通過實驗驗證的假設(shè),即異種雜交個體具有高的突變率[40]。根據(jù)已有的文獻報道[39,41-42],推測騾基因組中較高的denovoSNP頻率,是由于馬和驢單倍體整合到騾基因組的應(yīng)激引發(fā)的合子后快速突變所致。但本研究的數(shù)據(jù)不足以解釋其詳細(xì)分子機制。

研究表明,用相同的方法識別不同品種馬的CNV,其數(shù)量變化范圍為從幾十到幾千[43]。而不同的方法對CNV識別的敏感度也有很大差別[44]。CNVnator軟件對CNV的識別具有高靈敏度(86%～96%),用該軟件對6個人深度測序(20×～32×)的illumina數(shù)據(jù)進行CNV分析,共識別了737～1 489個deletion CNVs[30]。使用CNVnator軟件,本研究識別的馬和驢deletion CNVs分別是1 543和3 223。與2014年Wang等[45]和2012年Doan等[46]的研究結(jié)果相比,本研究馬和驢CNVs數(shù)量相對較高,這可能與使用的分析方法、數(shù)據(jù)類型以及品種等有關(guān)。以驢基因組為參考識別的CNVs數(shù)高于以純血馬基因組為參考,該結(jié)果與韓紅梅[47]的研究結(jié)果趨勢一致,但與之相比本研究識別的騾CNVs數(shù)量更少,這可能是因為本研究對CNV的過濾更嚴(yán)格。與SNP相比,CNV結(jié)構(gòu)變異對機體具有更大的影響。馬和驢物種間的CNV以及相關(guān)基因間的遺傳距離,導(dǎo)致評估CNV對騾的影響非常復(fù)雜。

功能分析發(fā)現(xiàn)許多變異相關(guān)基因和機體的免疫過程相關(guān)。例如MHC基因家族,MHC基因是機體維持免疫穩(wěn)態(tài)和發(fā)揮適當(dāng)免疫功能的關(guān)鍵基因。與大多數(shù)哺乳動物一樣,馬屬動物MHC基因的特征是由病原微生物平衡選擇維持的極端多態(tài)性,和由不同地理環(huán)境導(dǎo)致的抗原結(jié)合位點的譜系特異性[48-50]。Masly等[12]對異種雜交小鼠的研究結(jié)果表明,MHC基因的種間多態(tài)性會降低其雜交后代的適應(yīng)性、存活率或繁殖成功率。許多在免疫過程中發(fā)揮關(guān)鍵作用的通路,在本研究中被富集,包括抗原加工和提呈過程、TCR信號傳遞的2個重要下游信號通路:MAPK信號通路和NF-κB信號通路[51]、以及在T細(xì)胞分化過程中具有重要調(diào)控作用的JAK-STAT信號通路[52]。王艷欣[21]對馬騾和驢騾不同組織轉(zhuǎn)錄組的研究發(fā)現(xiàn),許多結(jié)構(gòu)變異相關(guān)基因富集到了免疫相關(guān)通路。韓紅梅[47]對馬騾和驢騾基因組的研究發(fā)現(xiàn),許多結(jié)構(gòu)變異相關(guān)基因和機體免疫過程相關(guān)。騾部分結(jié)構(gòu)變異相關(guān)基因還參與機體的癌癥過程,例如,原癌基因(如HRAS)和腫瘤抑制基因(如APC和PTEN)的突變是各種癌癥發(fā)生的主要候選基因[53-55]。Liu等[20]對雜交魚的研究發(fā)現(xiàn),許多結(jié)構(gòu)變異相關(guān)基因和癌癥過程相關(guān)。綜上所述,這些突變可能反映了馬和驢之間進化累積的不相容位點,或者是反映了馬和驢雜交在基因組序列水平是不利于騾適應(yīng)性的[7]。另一方面,突變產(chǎn)生的結(jié)構(gòu)變異,可能通過改變酶活性或?qū)Ｒ恍?影響基因的正常轉(zhuǎn)錄及最終的功能效能[56]。因此,快速突變的SNP也可能作為一種“緩沖劑”來平衡親本單倍體基因組的不相容,從而提高騾的適應(yīng)性[7]。因此,在接下來的工作中,需要在更大樣本量、更高測序深度和更多試驗開展的基礎(chǔ)上,來提高SNP和CNV的陽性率和進一步驗證以上推測。

4 結(jié) 論

本研究分析了馬屬動物三成員家系之間結(jié)構(gòu)基因組的遺傳多樣性,識別了騾異常高的denovoSNP突變率,和騾的特異性CNV。大部分突變相關(guān)基因與機體的免疫過程相關(guān)。結(jié)果說明,馬和驢雜交作為一種應(yīng)激,可能觸發(fā)騾基因組發(fā)生高頻率的突變。另外,這些突變可能反映了馬和驢的遺傳不相容位點,也可能對騾的適應(yīng)性具有重要意義。這些結(jié)構(gòu)變異的識別,為進一步開展馬和驢異種雜交的遺傳基礎(chǔ)及其分子機制的研究提供候選遺傳位點。