永登七山羊全基因組選擇信號檢測分析

栗登攀,馬克巖,韓金濤,白雅琴,李討討,馬友記*

(1.甘肅農(nóng)業(yè)大學(xué)動物科學(xué)技術(shù)學(xué)院,蘭州 730070;2.永登縣動物疫病預(yù)防控制中心,蘭州 730300; 3.甘肅省畜牧技術(shù)推廣總站,蘭州 730030)

我國地方綿羊品種資源豐富,為切實挖掘新資源并加以利用,在全國第三次畜禽資源普查中發(fā)現(xiàn)了永登七山羊。永登七山羊在永登縣七山鄉(xiāng)各村落均有分布,該群體具有特征一致、特性明顯、遺傳性能穩(wěn)定等特點。研究發(fā)現(xiàn),永登七山羊與蘭州大尾羊、灘羊和岷縣黑裘皮羊之間分化明顯,獨立形成聚類[1]。

在動物漫長的馴化過程中,動物的行為習(xí)慣、體型外貌因受到自然選擇而不斷發(fā)生著巨大的變化[2-3]。使用選擇信號來進行動植物的適應(yīng)性進化研究越來越多[4-6]。選擇信號是指動物在進化過程中受長期的自然選擇和強烈的人工選擇兩種作用而遺留在基因組上的遺傳標(biāo)記,從基因組水平探討種群間的遺傳差異,有助于理解動物對選擇適應(yīng)的遺傳機制和挖掘重要經(jīng)濟性狀的候選基因,為種質(zhì)資源的合理利用、保存提供理論依據(jù)[7-8]。Pan等[9]通過對89個中國地方綿羊品種進行全基因重測序揭示了RXFP2基因與角的獨特形狀相關(guān)。Li等[10]通過選擇多種尾型的品種羊進行選擇信號掃描,利用窗口XP-CLR和π ratio的方法,篩選出一系列(SGCZ、GLIS3、FGF7、PAPPA2、MAP2K3等)與脂肪沉積與代謝相關(guān)的基因。Li等[11]收集分析了中國、巴基斯坦和尼泊爾的15個代表性家山羊品種共115個樣本以及4個野山羊樣本的全基因組數(shù)據(jù),篩選出調(diào)控山羊熱適應(yīng)、產(chǎn)絨和產(chǎn)奶的關(guān)鍵候選基因與SNP錯義位點。這些研究結(jié)果有助于人們深刻了解動物受自然選擇和人工選擇的影響,并為物種起源及馴化機制提供見解。目前,對于永登七山羊的研究僅存在于遺傳多樣性分析上,需要進一步從基因組水平篩選鑒定可能與永登七山羊重要經(jīng)濟性狀相關(guān)的功能基因。

本研究利用SLAF技術(shù)檢測4種綿羊群體全基因組范圍內(nèi)的SNPs。通過主成分分析、基因流事件以及候選基因挖掘等方法,旨在從分子水平證明永登七山羊群體與其他群體的基因交流情況,并對適應(yīng)性進化過程中的優(yōu)勢基因進行挖掘,從而為永登七山羊的保護利用提供重要的理論參考。

1 材料與方法

1.1 試驗材料

本研究基于前期針對地方品種(永登七山羊)的遺傳多樣性分析研究[1],選擇蘭州大尾羊(LFT)、灘羊(TAN)、岷縣黑裘皮羊(MBF)作為參考對象。如表1所示,每個品種采集10只健康的成年母羊,頸部靜脈采血3～5 mL,用5 mL EDTA抗凝管收集,置于-80 ℃冰箱保存?zhèn)溆谩?/p>

表1 樣品信息Table 1 Sample information

1.2 試驗方法

對抗凝血樣提取DNA,利用NanoDrop 2000和Qubit 2.0檢測DNA的純度和濃度。采用瓊脂糖凝膠電泳檢測基因組 DNA 的完整性。對DNA質(zhì)檢合格的樣品采用雙酶切的RAD建庫方式,將長度范圍414～464 bp的酶切片段定義為SLAF標(biāo)簽,對SLAF標(biāo)簽進行處理后在Illumina HiSeqTM 2500測序儀上進行簡化基因組測序。測序由北京百邁客生物科技有限公司完成。

1.3 數(shù)據(jù)處理與分析

對測序得到的原始數(shù)據(jù)按照酶切效率、測序質(zhì)量和片段選擇進行數(shù)據(jù)質(zhì)量評估,得到各個樣品的原始序列片段。其中限制性內(nèi)切酶組合為Hpy166 II+EcoR V,采用讀長126 bp×2 作為后續(xù)的數(shù)據(jù)評估和分析數(shù)據(jù)。利用bwa[12]將原始數(shù)據(jù)比對到綿羊參考基因組(oar_v4.0)上,使用GATK[13]和samtools[14]兩種方法開發(fā) SNP,取SNP標(biāo)記交集并使用snpeff軟件[15]用于注釋變異。

1.4 主成分分析

使用EIGENSOFT[16]軟件包中的smart PCA程序基于SNP數(shù)據(jù),在PCA二維聚類圖[1]基礎(chǔ)上重新進行數(shù)據(jù)分析,得到4個綿羊群體的三維聚類情況。

1.5 基因流分析

利用Treemix軟件,使用全基因組范圍內(nèi)等位基因頻率數(shù)據(jù),推斷多個群體分裂和混群模式。使用R的OptM判斷最優(yōu)m值,m值取1～10。每個m值重復(fù)5次后,使用最優(yōu)m值進行 Treemix 分析,繪制群體最大似然樹;并在最大似然樹上用箭頭展示滲入的遷移事件。

1.6 選擇清除分析

基于SNP數(shù)據(jù),按照次要等位基因頻率(MAF:0.05)和位點完整度(INT:0.8)進行過濾,獲得高一致性SNP位點,按照100 kb的窗口、10 kb的步長對染色體進行選擇性清除區(qū)域檢測。使用PopGenome軟件包計算群體分化指數(shù)(Fst)、核苷酸多態(tài)性(π)和多樣性變化倍數(shù)(θπ ratio)等群體遺傳學(xué)指標(biāo)。在篩選選擇性清除區(qū)域時,分別以Fst和π ratio前5%分位數(shù)對應(yīng)的數(shù)據(jù)作為閾值,取該區(qū)域的交集為選擇性清除區(qū)域。

1.7 受選擇區(qū)域篩選及基因富集分析

針對選擇性清除區(qū)域進行相鄰區(qū)域窗口合并,統(tǒng)計合并后窗口的位置信息以及窗口內(nèi)基因數(shù)目,對應(yīng)的基因位置。使用Enesmbl數(shù)據(jù)庫的BioMart功能,利用功能分類和代謝途徑信息對候選基因進行基因功能注釋。利用DAVID數(shù)據(jù)庫對基因進行 GO(Gene Ontology)和 KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析,并分別生成富集分析圖表。

2 結(jié) 果

2.1 SNP檢測

本項目共獲得 709 067 個SLAF標(biāo)簽,標(biāo)簽的平均測序深度為 17.63x,其中,多態(tài)性SLAF標(biāo)簽有 445 768 個,共獲得 1 658 596 個SNPs標(biāo)記。根據(jù)SLAF在染色體上的分布,繪制SLAF在染色體上的分布圖(圖1)。

2.2 主成分分析

根據(jù)主成分分析推斷出群體之間聚類模式的差異。由4個綿羊群體的主成分分析(圖2)可明顯看出,PC1、PC2和PC3的貢獻率分別為3.86%、3.43%和3.39%;LFT群個體最為聚集,MBF群體最為分散;MBF、QS和TAN群個體可獨立分群。

圖2 4個綿羊群體的主成分分析Fig.2 Principal component analysis of 4 sheep populations

2.3 基因流分析

基因流動是指位于不同種群之間或者在種群內(nèi)部產(chǎn)生的遺傳物質(zhì)的交流[17-18],既能抵消由遺傳變異和自然選擇引起的自我種群之間的分化,又能給整個大的群體中帶來新的變異。從圖3可以看出,TAN和LFT群體間存在較強的基因交流,QS和LFT群體間發(fā)生較弱的基因交流。

圖3 基因流動Fig.3 Gene flow

2.4 群體間受選擇區(qū)域篩選

將QS與LFT、TAN和MBF進行兩兩比較并將其余3個群體數(shù)據(jù)合為整體與QS進行比較。在QS與MBF之間,將Fst>0.273, π ratio>8.55的基因組交集區(qū)域定義為QS與MBF群體間選擇信號(圖4A),檢測到163個受選擇窗口,注釋到424個候選基因(圖5A);以同樣的方法,在QS與LFT之間篩選到127個受選擇窗口(圖4B),注釋到294個候選基因(圖5B);在QS與TAN之間篩選到133個受選擇窗口(圖4C),注釋到301個候選基因(圖5C)。

在QS與(MBF +LFT+TAN)(簡稱:MLT)之間,篩選到202個受選擇窗口,并繪制群體間選擇信號示意圖(圖4D);注釋到466個候選基因(圖5D)。另外,統(tǒng)計窗口中各個染色體的Fst和π ratio值,以及窗口內(nèi)各個基因的位置,繪制QS與MLT SNPs注釋表(表2)。

A.QS與MBF組;B.QS與LFT組;C.QS與TAN組;D.QS與MLT組 A.QS and MBF group;B.QS and LFT group;C.QS and TAN group;D.QS and MLT group圖4 群體間選擇信號示意圖Fig.4 Schematic of selection signatures between populations

A.QS與MBF組;B.QS與LFT組;C.QS與TAN組;D.QS與MLT組A.QS and MBF group;B.QS and LFT group;C.QS and TAN group;D.QS and MLT group圖5 兩種方法檢測的正選擇基因韋恩圖Fig.5 Venn diagram for the positively selected genes identified via two approaches

2.5 GO富集分析

對候選基因進行GO功能富集分析。在QS與MBF組間424個候選基因中,119個候選基因顯著富集于65個GO條目(P<0.05),其中包括9條細(xì)胞組分類,43條生物過程類和13條分子功能類,主要集中在細(xì)胞極性的建立或維持、肺泡發(fā)育、組蛋白H3-K4三甲基化、血管重塑等。在QS與LFT組294個候選基因中,119個候選基因顯著富集于79個GO條目(P<0.05),其中包括54條生物過程類,18條分子功能類和7條細(xì)胞組分類,主要富集在mRNA穩(wěn)定化、前置/后置模式規(guī)范、有絲分裂期板塊聚集等。在QS與TAN組301個候選基因中,90個基因顯著富集于41個GO條目(P<0.05),其中包括33條生物過程類和8條分子功能類,主要富集在煙酸和煙酰胺的代謝、金黃色葡萄球菌感染、雌性激素信號傳導(dǎo)途徑、初級膽汁酸的生物合成等。

另外,在QS與MLT組466個候選基因中,130個候選基因顯著富集在124個GO條目(P<0.05),其中包括96條生物過程類,14條細(xì)胞組分類和14條分子功能類,主要富集在多細(xì)胞進程、血液循環(huán)、發(fā)育生長調(diào)節(jié)、鈣離子轉(zhuǎn)運的正調(diào)節(jié)等過程。富集條目、通路名稱、P值以及對應(yīng)基因見表3。

表3 GO功能富集分析Table 3 GO function enrichment analysis

2.6 KEGG富集分析

對候選基因進行KEGG通路富集分析,并挑選前10個顯著的KEGG通路作富集氣泡圖(圖6)。結(jié)果顯示,QS群體與MBF、LFT、TAN品種相比較,分別有231條、222條和214條信號通路參與主要生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑。其中,QS與MBF比較組中顯著富集通路有15條,其中15個候選基因(GAD1、CYP11B1、PTAFR、CAMK2D、CDH1、BMP2等)可能涉及重要經(jīng)濟性狀,主要參與類固醇激素生物合成、葉酸的生物合成、雌激素信號傳導(dǎo)途徑、金黃色葡萄球菌感染、卵巢類固醇生成、花生四烯酸代謝、精氨酸和脯氨酸的代謝等通路。QS與LFT比較組中顯著富集通路有22條,涉及12個候選基因(ALDH1A1、GAD1、GRM1、PTAFR、CYP11B1等),主要集中在雌性激素信號傳導(dǎo)途徑、金黃色葡萄球菌感染、視黃醇代謝、丙氨酸、天門冬氨酸和谷氨酸的代謝、苯丙氨酸的代謝、色氨酸的代謝、類固醇激素的生物合成等。QS與TAN比較組中顯著富集通路有10條,涉及13個候選基因(為ALDH1A1、KMT2A、PARD6G、SETD1A、CAMK2D、GRM1等),主要集中在煙酸和煙酰胺的代謝、金黃色葡萄球菌感染、雌性激素信號傳導(dǎo)途徑、初級膽汁酸的生物合成、氮素代謝、粘著劑的交界處、賴氨酸的降解、視黃醇代謝、TGF-beta信號傳導(dǎo)途徑、軸突引導(dǎo)、嘧啶的代謝等。

A.QS vs MBF組;B.QS vs LFT組;C.QS vs TAN組;D.QS vs MLT組A.QS and MBF group;B.QS and LFT group;C.QS and TAN group;D.QS and MLT group圖6 KEGG 通路分析Fig.6 KEGG pathway analysis

QS與MLT比較組中,參與主要生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑有244條。其中PTAFR、UGT2A1、KRT15等基因主要參與類固醇激素生物合成、葉酸合成、雌激素信號通路、金黃色葡萄球菌感染、卵巢類固醇生成、甘露糖型O-聚糖生物合成、SNARE在膀胱運輸中的相互作用等。

值得注意的是,雌激素信號通路和金黃色葡萄球菌感染通路均顯著存在于4個組別中,并從中均發(fā)現(xiàn)V15、K38、KRT32、KRT35、KRT36這5個基因。

3 討 論

在動物遺傳學(xué)上,基于不同群體的SNP差異程度,可將各群體進行分群。本研究中,通過PCA分析,發(fā)現(xiàn)岷縣黑裘皮羊、永登七山羊和灘羊均可與其他群體分開,除永登七山羊外,其他群體與之前一些學(xué)者的研究結(jié)果一致[19-20]。

基因交流可以向其他群體中引入新的等位基因,進而影響群體遺傳多樣性,產(chǎn)生新的性狀組合,發(fā)生遺傳變異[21-23]。本研究發(fā)現(xiàn),灘羊和蘭州大尾羊之間發(fā)生過相對較強的基因交流,永登七山羊與蘭州大尾羊之間發(fā)生了較弱的基因交流。其原因是灘羊、蘭州大尾羊、永登七山羊均屬蒙古羊,雖然永登七山羊飼養(yǎng)環(huán)境相對封閉,但因為和蘭州大尾羊地理上接近,難免發(fā)生基因交流;也是自然選擇作用的表現(xiàn)。因此保護這些品種生物多樣性具有很大難度,可以采用等量留種的方式,盡量保持現(xiàn)有等位基因,防止近交衰退。

本研究取top 5%Fst和π ratio的交集用來確定基因組受選擇區(qū)域,并對篩選出來的候選基因進行GO和KEGG富集分析。采用Fst和π ratio的交集來進行篩選鑒定,能夠得到更加可靠的結(jié)果[24-25]。本試驗中,各比較組的受選擇窗口具有一定差異,且注釋到的基因均隨著受選擇窗口的數(shù)量而變化。通過SNPs注釋,發(fā)現(xiàn)了可能與永登七山羊重要經(jīng)濟性狀相關(guān)功能的20個候選基因,包括CDH1、V15、K38、KRT15、KRT32、KRT35、KRT36、BMP2、GRM1、ALDH1A1、CAMK2D、MASP1、LUC7L3、UGT2A1、PTAFR、PLOD1、SETD1A、PARD6G、KMT2A、SOD1。其中與綿羊生長發(fā)育有關(guān)的功能基因有BMP2、GRM1和ALDH1A1。在MBF與QS比較組中篩選出了BMP2基因,發(fā)現(xiàn)BMP2基因存在于Hippo信號通路、基底細(xì)胞癌、轉(zhuǎn)化生長因子-β(TGF-β)和細(xì)胞因子-細(xì)胞因子受體相互作用的信號通路中。BMP2基因影響羊生長性狀和脂肪酸組成[26-27],在肌肉發(fā)育中發(fā)揮重要作用[28]。與脂肪尾表型,脂質(zhì)代謝相關(guān),參與了羊尾脂肪沉積[29-31]。除此之外,BMP2可能是影響綿羊產(chǎn)仔數(shù)的候選基因[32]。在灘羊與永登七山羊比較組和蘭州大尾羊與永登七山羊比較組均篩選出GRM1和ALDH1A1基因。發(fā)現(xiàn)GRM1基因存在于雌激素信號通路、磷脂酶D信號通路、逆行內(nèi)源性大麻素信號、長時程增強、長時程抑制、間隙連接的信號通路,GRM1基因被認(rèn)為是影響斷奶后增益的重要候選基因,對綿羊生長和肉類生產(chǎn)性狀有一定影響[33-35],可以作為綿羊肉質(zhì)性狀的候選基因[36];又發(fā)現(xiàn)ALDH1A1基因存在于視黃醛代謝信號通路。ALDH1A1基因是綿羊脂尾進化關(guān)鍵基因[37],可能參與灘羊脂肪細(xì)胞分化調(diào)控,敲除灘羊尾部脂肪前體細(xì)胞中ALDH1A1基因能夠抑制脂肪細(xì)胞分化過程[38]。

另外,雌激素信號通路均顯著存在于4個組別中。雌激素介導(dǎo)的oar-miR-485-5p靶向PPP1R13B調(diào)控綿羊成肌細(xì)胞增殖,參與炎癥反應(yīng)、新陳代謝和干細(xì)胞存活[39-40]。雌激素與雌激素受體結(jié)合后可以調(diào)節(jié)細(xì)胞核中的基因轉(zhuǎn)錄,或者激活細(xì)胞質(zhì)中的激酶進而發(fā)揮其作用[41],且雌激素介導(dǎo)的信號通路與毛囊周期性生長和發(fā)育有關(guān)[42]。V15、K38、KRT32、KRT35、KRT36均作為角蛋白,廣泛存在于綿羊毛發(fā)中[43-44]。其中KRT36是影響羊毛彎曲性狀的潛在候選基因[45]。在毛球下部至中部,KRT32和KRT35基因在角質(zhì)層和纖維皮質(zhì)層均有表達[46]。KRT32在毛球中部至角質(zhì)生長層皮質(zhì)中廣泛表達,KRT35基因和毛囊透明蛋白參與初始分化的皮質(zhì)層結(jié)構(gòu)的建立,最終導(dǎo)致纖維彎曲[47-48]。本研究采用簡化基因組測序檢測永登七山羊選擇信號,注釋到大量候選基因。在候選基因功能顯著性富集分析中,功能基因揭示了永登七山羊在尾部脂肪沉積、肌肉發(fā)育和毛發(fā)生長等方面受到了選擇。

4 結(jié) 論

本研究通過主成分分析和基因交流揭示4個群體的進化關(guān)系,結(jié)果表明,永登七山羊群體與其他3個群體分化明顯,且與蘭州大尾羊發(fā)生基因交流;此外還篩選到與綿羊生長發(fā)育有關(guān)的功能基因BMP2、GRM1和ALDH1A1等。結(jié)果為挖掘永登七山羊生長發(fā)育等相關(guān)基因提供參考,為其遺傳資源保護與開發(fā)利用提供科學(xué)依據(jù)。