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

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

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

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

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

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

1 材料與方法

1.1 實驗動物和測序數據

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

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

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

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

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

1.4 騾de novo SNP的識別

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

表1 De novo SNPs的基因型組合

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

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

使用CNVnator軟件及其推薦參數識別CNV[30]。使用以下標準過濾原始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個拷貝數變異區(qū)域(CNVR)。

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

1.6 基因富集分析

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

2 結果與分析

2.1 驢、馬和騾測序數據及比對

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

表2 參考基因組比對結果

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

2.2.1驢、馬和騾InDel和SNP識別

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

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

以驢基因組為參考,分別識別驢、馬和騾的高質量SNPs為3 212 499、23 549 224和21 870 390。馬(0.968 3%)和騾(0.899 3%)的SNP頻率明顯高于驢(0.132 1%)。實驗室已發(fā)表數據表明,以純血馬基因組為參考,分別識別驢、馬和騾的高質量SNPs為23 819 055、5 012 403和23 426 241;蒙古馬共識別5 012 403個SNPs,雜合SNP數高于純合,驢(0.950 1%)和騾(0.934 4%)的SNP頻率明顯高于馬(0.199 9%)[22]。以上結果表明,由于驢、馬和騾與參考基因組物種的親緣關系不同,導致與參考基因組不同物種識別的SNPs頻率遠高于相同物種(表3)。以純血馬基因組為參考識別的騾SNPs數略高于以驢基因組為參考的SNPs識別結果。由圖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個等位基因不同于親本。對上述識別的高質量SNP,進一步屏蔽掉重復序列和CNV區(qū)域的SNP。以純血馬基因組和驢基因組分別為參考,與親本相比,騾分別識別了555和419個denovoSNPs。與SNP識別趨勢一致,騾denovoSNP的識別可能存在馬參考基因組的偏向性。本研究的三成員家系中,騾denovoSNP的頻率是1.72×10-7～2.21×10-7(DenovoSNPs/全基因組序列)。

2.2.3騾denovoSNP相關基因注釋

為研究騾denovoSNP的潛在功能,使用ANNOVAR軟件,對SNP上下游5 kb范圍內的基因進行注釋。以純血馬基因組為參考,555個denovoSNPs共注釋658個相關基因。以驢基因組為參考,419個denovoSNPs共注釋540個相關基因。DenovoSNP均最高頻率的注釋于基因間區(qū)(表4)。以兩套基因組為參考注釋denovoSNP,均發(fā)現大部分相關基因直接參與機體的免疫功能。例如,純血馬參考基因組注釋的基因,包括MHC I類基因(如MHCX1和LOC100049798)和MHC II類基因(如DQA和DQB)、Src家族酪氨酸激酶(如FYN)、表皮生長因子受體基因(EGFR)以及癌癥相關基因(如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相關基因富集分析

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

圖3 De novo SNP相關基因的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相關基因功能分析

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

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

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

3 討 論

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

結合以純血馬基因組為參考識別騾denovoSNPs(555)的結果可知[22],騾denovoSNPs數遠高于用相同方法識別的純種黑猩猩的denovoSNPs數(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年通過實驗驗證的假設,即異種雜交個體具有高的突變率[40]。根據已有的文獻報道[39,41-42],推測騾基因組中較高的denovoSNP頻率,是由于馬和驢單倍體整合到騾基因組的應激引發(fā)的合子后快速突變所致。但本研究的數據不足以解釋其詳細分子機制。

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

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

4 結 論

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