山東小毛驢全基因組選擇信號(hào)檢測

陳建興，童家興，張孝忠，張向陽，王宇鑫，孫玉江,5

(1.赤峰學(xué)院 化學(xué)與生命科學(xué)學(xué)院，內(nèi)蒙古 赤峰 024000； 2.青島農(nóng)業(yè)大學(xué) 動(dòng)物科技學(xué)院，山東 青島 266109；3.內(nèi)蒙古東阿黑毛驢牧業(yè)有限公司，內(nèi)蒙古 赤峰 024328； 4.東阿阿膠股份有限公司，山東 聊城 252200；5.東營職業(yè)學(xué)院，山東 東營 257091)

山東小毛驢(SDL)過去主要產(chǎn)地在膠東半島、沂蒙山區(qū)和魯中平原[1]，舊稱膠東小毛驢，因現(xiàn)在中心產(chǎn)區(qū)在海陽等地，又稱海陽小毛驢[2]。山東小毛驢體質(zhì)外形與華北驢類同，屬小型驢，具有體型小、挽力大、耐粗飼、抗逆性強(qiáng)等特點(diǎn)[2]?，F(xiàn)代農(nóng)業(yè)機(jī)械化和交通業(yè)的快速發(fā)展，驢的役用地位迅速降低，使得山東小毛驢種質(zhì)資源迅速衰減，遺傳多樣性下降，育種潛力逐漸喪失。家畜遺傳資源問題是全球性生物資源問題的組成部分，與人類未來的生存與發(fā)展緊密相關(guān)[3]。驢產(chǎn)業(yè)又是我國的特色產(chǎn)業(yè)、民生產(chǎn)業(yè)和創(chuàng)新產(chǎn)業(yè)[4]，隨著我國社會(huì)經(jīng)濟(jì)發(fā)展和中國特色社會(huì)主義建設(shè)的不斷推進(jìn)，獨(dú)特的驢遺傳資源將會(huì)越來越珍貴，對山東小毛驢的種質(zhì)資源的保護(hù)開發(fā)將會(huì)日趨重要。

不同品種驢在體型外貌、生長性能、抗病力和適應(yīng)性等方面存在較大差異，這些差異是自然選擇與人工選擇共同作用于一些目標(biāo)基因定向選擇出來的結(jié)果，而選擇信號(hào)是與選擇相對應(yīng)的基因組信息，是選擇在基因組上留下的印跡，通常表現(xiàn)為受選擇的DNA片段或位點(diǎn)多態(tài)的降低或者基因的純合[5]。隨著單核苷酸多態(tài)性(SNP)芯片和第二代測序技術(shù)成本的不斷降低，利用選擇信號(hào)分析來揭示引起畜禽性狀表型差異的遺傳機(jī)制的相關(guān)研究日益普遍。呂世杰等[6]利用選擇性清除方法篩選郟縣紅牛和中國荷斯坦奶牛2個(gè)品種間差異的基因組區(qū)域，通過與動(dòng)物數(shù)量性狀座位(QTL)數(shù)據(jù)庫中牛繁殖性狀相關(guān)QTLs進(jìn)行比對，認(rèn)為CFDP1、CFDP2和FAM204A基因可優(yōu)先作為牛繁殖性狀相關(guān)候選基因。PETERSEN等[7]使用固定指數(shù)分析尋找了馬的基因組的選擇信號(hào)，研究結(jié)果表明，MSTN、ECA11和DMRT3基因受到了選擇，另外通過關(guān)聯(lián)分析發(fā)現(xiàn)，MSTN內(nèi)含子的1個(gè)SNP和啟動(dòng)子的1個(gè)InDel與肌纖維比例有關(guān)。樊英智[8]通過對我國6個(gè)不同類型的家驢品種群體(共57頭)進(jìn)行全基因組混池重測序，分析了不同類群基因組間的選擇信號(hào)，找到與驢品種表型如體尺、毛色等性狀有關(guān)聯(lián)的11個(gè)候選基因。

通過對不同驢品種進(jìn)行基因組重測序，檢測不同驢品種群體間基因組的選擇信號(hào)的差異，有助于揭示品種群體的進(jìn)化歷史，了解重要表型性狀形成的遺傳基礎(chǔ)。目前，對山東小毛驢基因組選擇信號(hào)的檢測還未見報(bào)道。利用群體遺傳分化系數(shù)(Fst)和核苷酸多樣性比值(π ratio)方法，對山東小毛驢與德州驢的三粉類群(DZS)、德州驢的烏頭類群(DZW)、廣靈驢(GL)、華北驢(NC)等4個(gè)驢群體間的選擇信號(hào)差異進(jìn)行分析，以期篩選出山東小毛驢的選擇信號(hào)區(qū)域，探討山東小毛驢特性的形成原因，旨在為山東小毛驢的保護(hù)和利用提供參考。

1 材料和方法

1.1 供試動(dòng)物及樣品

本研究采用4個(gè)驢品種(5個(gè)群體)共60頭驢作為供試動(dòng)物(表1)。對所有樣本均采集頸靜脈血液10 mL于EDTA抗凝管中，置于-20 ℃冰箱短暫保存后于干冰保鮮盒中快遞至美吉生物公司(上海)進(jìn)行后續(xù)研究。

1.2 測序數(shù)據(jù)

利用超聲波將檢測合格的樣品基因組DNA片段化形成隨機(jī)片段，對片段化的DNA進(jìn)行末端修復(fù)、3′端加A、連接測序接頭后，再利用磁珠吸附富集400 bp左右的隨機(jī)片段，經(jīng)PCR擴(kuò)增形成測序文庫。構(gòu)建好的測序文庫通過Illumina HiSeqTM平臺(tái)進(jìn)行測序，測序策略為Illumina PE150。測序數(shù)據(jù)已上傳到NCBI的SRA數(shù)據(jù)庫，SRA號(hào)為SAMN14484743—SAMN1448 4802。

1.3 SNP、InDel檢測與注釋

本研究采用BWA軟件[9]將高質(zhì)量測序數(shù)據(jù)比對到參考基因組序列 (https://www.ncbi.nlm.nih.gov/genome/7038)上，利用GATK軟件[10]進(jìn)行比對后校正，并進(jìn)行SNP和Small InDel標(biāo)記的檢測；利用SnpEff軟件[11]和參考基因組的基因預(yù)測信息進(jìn)行變異功能注釋，得到SNP、InDel的功能注釋信息。

表1 供試驢信息

1.4 Fst與π ratio計(jì)算

Fst計(jì)算公式：Fst=(MSP-MSG)/[MSP+(n-1)MSG][12]。其中，MSP為群體間均方差，MSG為群體內(nèi)均方差，n為校正后平均樣本大小。Fst可用來評(píng)價(jià)群體間的分化程度，該值越接近1說明兩群體間分化程度越高，越接近于0說明兩群體間分化程度非常有限。對于高質(zhì)量SNP (次等位基因頻率maf不低于0.05，缺失率miss為0，樣本測序深度不低于5)，利用VCFtools軟件[13]計(jì)算了兩兩群體間的Fst值(2 Mb窗口，10 kb步長滑窗)。Π=∑ijxixjπij，其中，xixj分別代表第i個(gè)和第j個(gè)序列的對應(yīng)頻率，πij則為2個(gè)序列之間不同位點(diǎn)所占百分比。Π是遺傳變異的1個(gè)量化值，用于表征某一種群多態(tài)性的強(qiáng)弱，通常用于衡量種群內(nèi)或種群間的多樣性，該值不依賴于樣本大小[14]。每個(gè)群體的π ratio值也使用VCFtools軟件[13]進(jìn)行計(jì)算，同樣是2 Mb窗口，10 kb步長滑窗。

1.5 群體間選擇信號(hào)檢測

本研究采用基于Fst和π ratio的方法對5個(gè)驢群體間的選擇信號(hào)進(jìn)行檢測。Fst和π ratio分別選取閾值0.95和0.05(分位數(shù))，關(guān)聯(lián)Fst和π ratio提取相應(yīng)候選區(qū)域(取重疊區(qū)域)，并提取相應(yīng)區(qū)域內(nèi)的變異位點(diǎn)信息。采用選擇性清除方法，分別對SDL與其余4個(gè)驢群體進(jìn)行比較分析，檢測到相應(yīng)的選擇信號(hào)區(qū)域。選擇信號(hào)強(qiáng)弱的判定，根據(jù)VCFtools軟件篩選出有選擇性消除位點(diǎn)相應(yīng)的注釋，選擇注釋級(jí)別為HIGH而排除MODERATE和LOW的結(jié)果。對選擇信號(hào)區(qū)域中存在的基因，通過在GeneCards數(shù)據(jù)庫 (www.genecards.org) 或NCBI Gene (www.ncbi.nlm.nih.gov/gene/?term=)中查詢基因注釋來確定基因功能。

2 結(jié)果與分析

2.1 5個(gè)驢群體變異檢測結(jié)果

5個(gè)驢群體的全基因組重測序，共獲得740.57 G高質(zhì)量數(shù)據(jù)，平均每個(gè)樣本獲得了12.34 G的數(shù)據(jù)。將高質(zhì)量測序數(shù)據(jù)比對到參考基因組之后，總共獲得了10 096 033個(gè) SNP(SDL、DZS、DZW、NC、GL分別占72.03%、77.01%、71.62%、74.99% 和70.34%) 和1 311 358個(gè)InDel(SDL、DZS、DZW、NC、GL分別占75.50%、79.71%、74.46%、76.86%和73.48%)。整體來看，2種變異(SNP和InDel)的最大比例均出現(xiàn)在德州驢的三粉類群中，而最小比例出現(xiàn)在廣靈驢群體中(表2)。然而，雜合SNP和InDel數(shù)目以及π ratio的最大值都出現(xiàn)在山東小毛驢群體中，而最低值都出現(xiàn)在華北驢群體中。觀測雜合度和雜合SNP數(shù)、雜合InDel數(shù)一樣，最大值出現(xiàn)在山東小毛驢群體中，最低值出現(xiàn)在華北驢群體中。轉(zhuǎn)換顛換比(Ts/Tv)的最大值出現(xiàn)在山東小毛驢群體中，最低值出現(xiàn)在德州驢的三粉類群中(表2)。

表2 5個(gè)驢群體的變異和遺傳多樣性指數(shù)

2.2 群體間遺傳分化分析

通過計(jì)算Fst，群體間的遺傳分化程度能夠得到估量，計(jì)算結(jié)果見表3。Fst值從DZS和NC群體間的0.007 80到DZW和SDL群體間的0.011 460。顯然，各Fst值都非常低，接近于0(WRIGHT認(rèn)為[15]，該值處于0～0.05，表明群體間遺傳分化很小，可以不用考慮)，表明這5個(gè)驢群體間的分化水平極低。

表3 5個(gè)驢群體間遺傳分化系數(shù)

2.3 山東小毛驢與其他4個(gè)驢群體的選擇性清除結(jié)果

選擇性清除指新的有利突變會(huì)增加其頻率并固定下來，導(dǎo)致其相鄰核苷酸序列的差異下降或消除的過程[16]。為了檢測到可能的選擇性消除位點(diǎn)，采用基于Fst和π ratio的方法搜索尋找了驢的基因組中的高度固定的區(qū)域，篩選過程見圖1，選擇Fst值大于0.95分位閾值并且π ratio值小于0.05分位閾值的區(qū)域，篩選的具體結(jié)果見表4。

藍(lán)色點(diǎn)是篩選的候選區(qū)域，F(xiàn)st值大于0.95分位閾值并且π ratio值小于0.05分位閾值

由表4可知，山東小毛驢與其他4個(gè)驢群體共檢測到正向選擇區(qū)域95個(gè)，德州驢的烏頭類群最多，共30個(gè)，華北驢最少，只有17個(gè)。為進(jìn)一步了解這些信號(hào)選擇區(qū)域的功能，對落入選擇信號(hào)區(qū)域的39個(gè)基因及基因功能進(jìn)行了生物信息學(xué)分析。山東小毛驢與德州驢的烏頭類群比較，群體間信號(hào)選擇區(qū)落入的基因數(shù)最多，有13個(gè)，與華北驢群體比較信號(hào)選擇區(qū)落入的基因數(shù)最少，只有5個(gè)(表4)。經(jīng)統(tǒng)計(jì)，與德州驢的三粉類群相比，群體間檢測出強(qiáng)選擇信號(hào)候選基因5個(gè)，分別是Fsip1、AHNAK2、CTAGE2、CYP3A12、LOC106830441；與德州驢的烏頭類群相比，群體間檢測出強(qiáng)選擇信號(hào)候選基因5個(gè)，分別是NKG2DL1、KLK1E2、CTAGE2、FAM170A、LOC106823932；與廣靈驢相比，群體間檢測出強(qiáng)選擇信號(hào)候選基因6個(gè)，分別是NKG2DL1、AHNAK2、CTAGE2、FAM170A、LOC106823932、LOC106848008；與華北驢相比，群體間檢測出強(qiáng)選擇信號(hào)候選基因2個(gè)，分別是CTAGE2、FAM170A。落入強(qiáng)選擇信號(hào)的10個(gè)基因(表4中強(qiáng)選擇信號(hào)基因列中顯示的所有基因)，大多與免疫、生殖以及細(xì)胞作用和代謝相關(guān)。

表4 山東小毛驢與其他驢群體信號(hào)選擇區(qū)域

3 結(jié)論與討論

本研究采用基于Fst和π ratio的方法對山東小毛驢和德州驢的三粉類群、德州驢的烏頭類群、廣靈驢、華北驢等5個(gè)驢群體進(jìn)行了選擇信號(hào)分析。SNP和InDel的最大比例出現(xiàn)在德州驢的三粉類群中，這可能與德州驢自古以來都是優(yōu)秀大型驢種有關(guān)[17]，一直與各地驢種之間存在交流，不斷有優(yōu)秀的個(gè)體引入該群體，而SNP和InDel的最小比例出現(xiàn)在廣靈驢群體中，這與廣靈驢目前處于保種狀態(tài)是相一致的。因?yàn)樘幱诒７N狀態(tài)，很少有其他驢品種與之交流，必然導(dǎo)致群體近交系數(shù)增大，純合度提高，而群體多樣性以及SNP和InDel的比例下降。然而，雜合SNP和InDel數(shù)目以及π ratio的最大值都出現(xiàn)在山東小毛驢群體中。盡管山東小毛驢各驢體高明顯較德州驢小，而且毛色多為灰色，但很可能與其地理位置相近的德州驢有雜交，也可能與地緣位置較近的河北省的一些驢種有基因交流，而雜合SNP和InDel數(shù)目以及π ratio的最低值都出現(xiàn)在華北驢群體中，可能與華北驢距離各驢種都較遠(yuǎn)，而且也可能與采集的樣本來自偏遠(yuǎn)的蒙古族村子有關(guān)，與其余驢種很少有雜交，幾乎沒有外源基因?qū)朐撊后w。

遺傳分化通常是生物長期進(jìn)化的產(chǎn)物，受交配系統(tǒng)、生物歷史和基因流等生物特征的影響。WRIGHT[15]認(rèn)為，如果Fst<0.050 0，各群體間就幾乎沒有分化。本研究中，各驢群體間的Fst值介于 0.007 80～0.011 460，明顯低于具有混合交配系統(tǒng)或幾年壽命的物種的Fst值，表明各驢群體在遺傳上很相似。這也凸顯了對這些驢品種進(jìn)行保護(hù)的必要性，不及時(shí)進(jìn)行保護(hù)，現(xiàn)今猶存的一些群體特征會(huì)隨著群體間雜交而逐漸消失。

Fsip1是一種生精細(xì)胞特異性表達(dá)蛋白，也是精子鞭毛纖維鞘的成分，在成年動(dòng)物的睪丸組織有超高水平的表達(dá)，其次為中樞神經(jīng)系統(tǒng)也有一定表達(dá)，已知參與組裝AKAP4輔助調(diào)節(jié)蛋白激酶A (PKA)。LIU等[18]發(fā)現(xiàn)，F(xiàn)sip1在HER2過表達(dá)型乳腺癌中能夠與HER2直接結(jié)合調(diào)控乳腺癌細(xì)胞的生長和侵襲。隨后，LIU等[19]發(fā)現(xiàn)，F(xiàn)sip1能夠通過誘導(dǎo)自噬，減少線粒體生成及增強(qiáng)和激活A(yù)MPK途徑來調(diào)節(jié)三陰性乳腺癌細(xì)胞的耐藥性。AHNAK2是AHNAK2基因編碼的1個(gè)較大的核蛋白，該蛋白質(zhì)可能通過與鈣通道蛋白相結(jié)合在鈣信號(hào)調(diào)節(jié)中發(fā)揮重要作用[20]。NKG2DL1是NKG2D這種免疫受體蛋白的配體，該受體因?yàn)榫哂锌共《竞涂鼓[瘤功能近些年引起了廣泛關(guān)注[21]。KLK1E2基因編碼的是絲氨酸蛋白酶的1個(gè)亞基，該酶具有廣泛的生理功能，越來越多的證據(jù)表明其參與癌癥發(fā)生，有些分子具有作為癌癥和其他疾病新的生物標(biāo)志物的潛力。有研究表明，KLK1在流感病毒感染早期就干預(yù)抗病毒防御，調(diào)節(jié)流感感染的嚴(yán)重程度，慢性阻塞性肺病患者KLK1表達(dá)降低可能導(dǎo)致流感惡化[22]?；蛟S受到選擇的這些基因?qū)τ谏綎|小毛驢適應(yīng)嚴(yán)酷環(huán)境的能力至關(guān)重要，這些基因才會(huì)在山東小毛驢群體內(nèi)逐漸固定下來。這暗示山東小毛驢可能在免疫相關(guān)過程中經(jīng)歷了一定的選擇作用，對揭示山東小毛驢免疫性狀的遺傳機(jī)制具有一定的意義。這也說明相對于其他驢種，山東小毛驢在抗病力強(qiáng)、繁殖力強(qiáng)這些優(yōu)良性狀上經(jīng)歷了較強(qiáng)的人工選擇。

本研究利用5個(gè)驢群體共計(jì)60個(gè)個(gè)體的全基因組重測序數(shù)據(jù)，對山東小毛驢群體選擇信號(hào)進(jìn)行了檢測分析，檢測到選擇信號(hào)的區(qū)域共計(jì)95個(gè)，落入這些區(qū)域并且選擇信號(hào)強(qiáng)的候選基因有10個(gè)。相比于其他4個(gè)驢群體，山東小毛驢在免疫、生殖等性狀上經(jīng)歷了人工選擇。

致謝：感謝內(nèi)蒙古阿魯科爾沁旗太極天驢集團(tuán)有限公司的鐘勇先生和舒蕾先生在采集華北驢樣本時(shí)的大力支持；感謝國家級(jí)廣靈驢保種場姜正廣先生和許增孝先生在采集廣靈驢樣本時(shí)的大力支持；感謝東阿阿膠股份有限公司嵇傳良先生在采集德州驢樣本時(shí)提供的大力支持；感謝山東海陽小毛驢保種繁育基地由松利先生在采集山東小毛驢樣本時(shí)提供的大力支持。