皖岳黑豬基因組遺傳變異分析及特征SNPs挖掘

王靜琳,劉陽光,徐啟隆,陳 朔,鄧在雙,程詩雨,丁月云,鄭先瑞,殷宗俊,張曉東

(安徽農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院,合肥 230036)

全基因組重測(cè)序是在已知物種參考基因組的條件下,通過對(duì)比測(cè)序序列與參考基因組之間的差異,以獲得單核苷酸多態(tài)性位點(diǎn)(single nucleotide poly-morphisms, SNP)、短插入/刪除(insertion and deletion, InDel)、大片段的結(jié)構(gòu)變異(structural variation, SV)和拷貝數(shù)變異(copy number variation, CNV)等遺傳變異的方法。這些變異位點(diǎn)作為分子遺傳標(biāo)記,在人類復(fù)雜疾病、進(jìn)化、動(dòng)植物經(jīng)濟(jì)性狀和育種研究等方面具有重大意義[1]。已有多項(xiàng)研究通過全基因組重測(cè)序技術(shù)獲得遺傳變異以解析動(dòng)物基因組遺傳機(jī)制[2-6]。SNPs是研究動(dòng)植物品系遺傳變異的重要依據(jù)[7-8],其突變可能會(huì)改變基因的功能,進(jìn)而導(dǎo)致生物性狀改變甚至致死。SVs是基因組多樣性和基因表達(dá)變異的重要來源,對(duì)基因表達(dá)的影響較大[9-11]。CNV等結(jié)構(gòu)變異是哺乳動(dòng)物基因組變異的重要來源,覆蓋更寬的染色體區(qū)域,可能導(dǎo)致基因結(jié)構(gòu)的變化、基因調(diào)控的修飾、基因劑量的變化以及隱性等位基因的暴露,從而導(dǎo)致大的表型效應(yīng)[12]。

皖岳黑豬以大別山腹地安徽省岳西縣所產(chǎn)淮豬為母本,北京黑豬為第一父本,杜洛克豬為第二父本,采用雜交、回交制種及橫交固定方法,正在培育的淮豬和北京黑豬血統(tǒng)各占37.5%,杜洛克血統(tǒng)占25%的優(yōu)質(zhì)黑豬新品種,現(xiàn)已完成5個(gè)世代的持續(xù)選育并初步進(jìn)行了中間試驗(yàn),下一步擬申請(qǐng)國家新品種審定。目前關(guān)于皖岳黑豬基因組結(jié)構(gòu)變異方面的相關(guān)研究尚未見報(bào)道。

本研究利用皖岳黑豬全基因組重測(cè)序數(shù)據(jù),檢測(cè)基因組遺傳變異,揭示皖岳黑豬的群體遺傳結(jié)構(gòu),并與北京黑豬、杜洛克、大白豬、藍(lán)塘豬、民豬、深縣黑豬的SNP位點(diǎn)構(gòu)建數(shù)據(jù)集,利用挑選最大分類能力和機(jī)器學(xué)習(xí)方法進(jìn)行基因組特征SNP位點(diǎn)的挖掘,確定皖岳黑豬品種特征SNP位點(diǎn),為皖岳黑豬進(jìn)一步的選育開發(fā)提供科學(xué)依據(jù)。

1 材料與方法

1.1 樣品采集與DNA提取

在安徽省岳西縣徽名山皖岳黑豬保種場(chǎng)隨機(jī)選取代表全部血統(tǒng)且無血緣關(guān)系的24頭體重達(dá)110 kg 的皖岳黑豬,采集耳組織樣品,并使用磁珠法進(jìn)行DNA提取,用0.8%瓊脂糖凝膠(25 min,170 V)和納米滴分光光度計(jì)(NanoDrop 2000)評(píng)估DNA的質(zhì)量和濃度。

1.2 全基因組重測(cè)序分析

采用Illumina NovaSeq 6000平臺(tái)(Illumina,圣地亞哥,加利福尼亞州,美國)對(duì)24頭皖岳黑豬進(jìn)行全基因組重測(cè)序。此外,從公共數(shù)據(jù)庫下載22頭北京黑豬、22頭杜洛克、16頭大白豬、5頭藍(lán)塘豬、5頭民豬、39頭深縣黑豬的基因組數(shù)據(jù)作為對(duì)照。

使用fastp(v0.20.1)[13]過濾原始重測(cè)序數(shù)據(jù),刪除N含量超過該read堿基數(shù)的10%和低質(zhì)量(Q≤5)堿基數(shù)超過該條read堿基數(shù)的50% 的paired reads。利用BWA[14]軟件 將原始測(cè)序數(shù)據(jù)比對(duì)到豬參考基因組(SusScrofa11.1),使用GATK[15]軟件檢測(cè)SNP變異并對(duì)其進(jìn)行過濾,保留質(zhì)量得分>30,MQ RMS映射質(zhì)量>20,dp>5,覆蓋率>30%,最小等位基因頻率(MAF)>0.01的位點(diǎn)。使用SnpEff[16]軟件注釋SNP,統(tǒng)計(jì)SNP各變異類型發(fā)生的比例及在各染色體上的位置;DELLY[17]進(jìn)行SV檢測(cè),去除Y染色體上的變異并過濾掉檢測(cè)出的低質(zhì)量和<50 bp的SVs。ANNOVAR(v2019)[18]軟件對(duì)過濾后的變異進(jìn)行注釋并統(tǒng)計(jì)SV在基因組各區(qū)域的分布情況。使用CNVnator和CNVcaller共同檢測(cè)皖岳黑豬群體的CNV變異情況,剔除p1>0.01、Q0>0.5的CNV,將重疊至少1 bp的CNV合并為一個(gè)CNVR區(qū)域,只保留在3個(gè)及以上個(gè)體中檢測(cè)到的CNVR,用CNVnator檢測(cè),再用CNVcaller直接檢測(cè)得到群體的CNVR情況,取兩個(gè)軟件檢測(cè)出的CNVR交集作為皖岳黑豬群體的真實(shí)CNVR,將兩個(gè)軟件檢測(cè)出的CNVR重疊區(qū)域大于50%認(rèn)為是同一個(gè)CNVR,并合并為新的CNVR。

1.3 皖岳黑豬群體遺傳多樣性分析

使用Plink(v 1.90)[19]進(jìn)行遺傳多樣性分析。群體的遺傳多樣性是保證物種得以生存和進(jìn)化的關(guān)鍵因素。對(duì)皖岳黑豬保種群體遺傳多樣性分析主要包括有效群體大小(Ne)、群體的期望雜合度(He)和觀察雜合度(Ho)、群體多態(tài)性標(biāo)記比例(PN)、群體近交系數(shù)。有效群體大小反映了該群體在進(jìn)化過程中受到的遷移、雜交的程度[20]。雜合度代表了群體遺傳多樣性的豐富度,分為期望雜合度(He)和觀察雜合度(Ho),當(dāng)He>Ho時(shí),推測(cè)群體發(fā)生了近交或選擇;如果He

1.4 皖岳黑豬遺傳距離分析

基于SNP分型結(jié)果,可以計(jì)算樣本間的IBS距離,IBS (identity by descent)叫做狀態(tài)同源,是指兩個(gè)個(gè)體中共有的等位基因序列相同。使用Plink(v 1.90)軟件中的plink--file hapmap1--cluster--matrix-noweb 指令計(jì)算遺傳距離,所得文件在R語言中結(jié)果可視化。

1.5 皖岳黑豬ROH近交系數(shù)

在Plink(V1.90)中調(diào)用--homozyg--homozyg-density 50--homozyg-gap 100--homozyg-kb 500--homozyg-snp 50--homozyg-window-het 1--homozyg-window-snp--homozyg-window-threshold 0.05-out指令計(jì)算每個(gè)樣本的ROH長(zhǎng)度,用個(gè)體ROH長(zhǎng)度除以基因組總長(zhǎng)計(jì)算得出每個(gè)個(gè)體的近交系數(shù)(FROH)。

1.6 基因組特征庫構(gòu)建及特異SNP位點(diǎn)功能注釋

本研究7個(gè)品種133個(gè)個(gè)體共檢測(cè)出1 35 915個(gè)SNPs位點(diǎn)。將全部個(gè)體分為訓(xùn)練集與測(cè)試集,其中訓(xùn)練集共有106個(gè)個(gè)體,測(cè)試集27個(gè)個(gè)體。對(duì)測(cè)試集利用挑選最大分類能力方法[8]進(jìn)行信息SNP的選擇,并結(jié)合機(jī)器學(xué)習(xí)方法對(duì)挑選出的SNPs進(jìn)行特征選擇,以獲取皖岳黑豬群體的特征位點(diǎn)。利用R(v 3.6.2)中的ClusterProfiler (v 3.14.0)[22]對(duì)篩選出的位點(diǎn)進(jìn)行注釋,并對(duì)所選位點(diǎn)注釋到的基因進(jìn)行GO和KEGG的富集分析,以P<0.05 作為判斷差異顯著的標(biāo)準(zhǔn)。

2 結(jié) 果

2.1 皖岳黑豬基因組變異檢測(cè)

皖岳黑豬原始重測(cè)序數(shù)據(jù)質(zhì)控后共獲得1 463.13 Gb的原始數(shù)據(jù),每頭獲得數(shù)據(jù)量平均約為60.96 Gb,平均比對(duì)深度為23.7×,平均比對(duì)率為98.4%,重復(fù)reads比率1.88%。SNP變異檢測(cè)統(tǒng)計(jì)結(jié)果顯示,共獲得31 534 384個(gè)SNPs,其中12.82%不存在于dbSNP庫。在外顯子區(qū)域的注釋中,沉默突變?cè)谒型蛔冾愋椭姓急茸畲鬄?3.91%,錯(cuò)義突變占35.62%(圖1A)。SNP在染色體上的密度分布如圖1B所示。

A.SNP檢測(cè)中各類型突變的比例;B. SNP在染色體上的密度分布:X軸代表染色體長(zhǎng)度,Y軸代表18條常染色體和X染色體,不同的顏色表示1 Mb內(nèi)包含的SNP個(gè)數(shù)A. Proportion of each type of mutation in SNP detection; B. Density distribution of SNPs on chromosomes: X-axis represents the chromosome length,Y-axis represents 18 autosomes and X chromosome, different colors indicate the number of SNPs contained in 1 Mb圖1 皖岳黑豬種群的SNP特征Fig.1 SNP characteristics of Wanyue black pig population

本研究共獲得43 673個(gè)SVs,變異總長(zhǎng)為354.36 Mb,覆蓋了豬全基因組的14.5%,變異類型及長(zhǎng)度所占比例表明,檢出的缺失變異最多,插入變異最少(圖2A);SVs在基因組各區(qū)域的分布表明,有3.6%的SVs落在基因外顯子區(qū)域(圖2B)。

A.SV類型的變異長(zhǎng)度分布比例;B. 皖岳黑豬SV變異在基因組各區(qū)域分布情況A.Proportion of variation length distribution of SV types; B. Distribution of SV variants in various regions of the genome of Wanyue black pigs 圖2 皖岳黑豬種群SV特征Fig.2 SV characteristics of Wanyue black pig population

本研究共檢測(cè)到3 258個(gè)CNVRs,其中“Loss”型2 427個(gè),“Gain”型831個(gè)。檢測(cè)出的CNVR在染色體上的分布如圖3所示,可以看出兩種類型的CNVR在染色體上的分布都較為均勻,并且“Loss”型CNVR的基因組覆蓋度遠(yuǎn)大于“Gain”型。CNVR在基因組的基因間區(qū)分布最多(47.5%),其次是基因組的內(nèi)含子區(qū)域(20.5%),見圖4。

圖3 CNVR在染色體上的分布Fig.3 Distribution of CNVR on chromosomes

圖4 CNVR在基因組各區(qū)域的分布Fig.4 Distribution of CNVR in various regions of the genome

2.2 皖岳黑豬群體遺傳多樣性分析

皖岳黑豬群體的Ne為4.2,He為0.320,Ho為0.328,PN為0.788。He與Ho近乎相等且He小于Ho,這一結(jié)果表明皖岳黑豬群體遺傳多樣性豐度較低且有其他血緣的引入,這也映照了皖岳黑豬在育種過程中有杜洛克及北京黑豬血緣的加入。但其有效群體大小依舊偏低,這就提示我們,在今后的選育方案中要考慮新的特別是公豬血緣的加入。

2.3 皖岳黑豬遺傳距離

皖岳黑豬群體IBS遺傳距離在0.157 4～0.287 3之間,平均為(0.243 5±0.022 2)。皖岳黑豬群體IBS距離矩陣的結(jié)果如圖5所示。皖岳黑豬個(gè)體間的IBS遺傳距離總體較遠(yuǎn),呈中等程度的親緣關(guān)系,部分個(gè)體間的IBS遺傳距離較近,存在較高的親緣關(guān)系,可能是皖岳黑豬在選育過程中出現(xiàn)了近交現(xiàn)象,在今后的選育中一定要多加注意,避免其發(fā)生近交衰退現(xiàn)象。

圖5 IBS遺傳距離矩陣Fig.5 IBS genetic distance matrix

2.4 皖岳黑豬ROH近交系數(shù)

利用Plink軟件計(jì)算皖岳黑豬個(gè)體的ROH,共檢測(cè)出973個(gè)ROH片段,平均每個(gè)個(gè)體中有40個(gè)ROH片段且個(gè)體總長(zhǎng)度為14.6～178.0 Mb,平均ROH長(zhǎng)度為(61±9.4)Mb。當(dāng)前整個(gè)群體的近交系數(shù)平均值為(0.025±0.004)。基于ROH的近交系數(shù)FROH的分布如圖6所示。

圖中心的白色圓點(diǎn)代表該群體FROH的中位數(shù),中間黑色長(zhǎng)方形方框的上緣和下緣分別為群體 FROH的上四分位數(shù)和下四分位數(shù)。小提琴圖的寬窄表示群體FROH的概率密度分布,小提琴圖越寬的部分表示處于該水平的樣本數(shù)目越多,反之則越少The white dot in the center of the figure represents the median of FROH in the group, and the upper and lower edges of the black rectangle in the middle are the upper and lower quartiles of the group FROH, respectively. The width of the violin plot indicates the probability density distribution of the population FROH, and the wider the violin plot, the larger the number of samples at that level, and vice versa圖6 基于ROH的近交系數(shù)分布圖Fig.6 Distribution of the inbreeding coefficient FROH based on ROH

2.5 基因組特異SNPs位點(diǎn)

通過機(jī)器學(xué)習(xí)方法最終選出33個(gè)SNP位點(diǎn)作為皖岳黑豬群體的特征位點(diǎn)。用訓(xùn)練集對(duì)選出的SNP位點(diǎn)進(jìn)行驗(yàn)證,其主成分分析結(jié)果如圖7所示,可以看出這33個(gè)SNPs位點(diǎn)可以很好的將皖岳黑豬個(gè)體與其他6個(gè)品種分開,用測(cè)試集測(cè)試這33個(gè)位點(diǎn)區(qū)分皖岳黑豬的準(zhǔn)確性,準(zhǔn)確率達(dá)到95%以上,證明這33個(gè)SNPs位點(diǎn)可以作為皖岳黑豬的特征庫。

圖7 基于10個(gè)特征SNPs的7個(gè)豬品種的主成分分析Fig.7 Principal component analysis of 7 pig breeds based on 10 characteristics SNPs

2.6 基因功能注釋

對(duì)位點(diǎn)進(jìn)行Ensembl數(shù)據(jù)庫檢索,33個(gè)SNPs位點(diǎn)映射到15個(gè)基因Ensemble號(hào),標(biāo)記到11個(gè)基因(表1),分別是SUSD4、GPC6、TENM2、NELL1、TMSB10、CHD1L、SLC41A2、KHDRBS2、CRISP1、NECTIN1、GRIK4。還有4個(gè)基因號(hào)未被標(biāo)記,但是同樣也參與生物過程。GO和KEGG進(jìn)行富集分析,對(duì)這些基因進(jìn)行功能評(píng)估,在84個(gè)GO項(xiàng)和5個(gè)KEGG途徑中發(fā)現(xiàn)顯著(FDR<0.05)的基因富集性,GO結(jié)果顯示,這些基因富集到的生物過程(biological process, BP)為 39 個(gè),細(xì)胞組分(cellular component, CP)為 34 個(gè),分子功能(molecular function, MF)為 11個(gè)(圖8)。主要參與成骨細(xì)胞活性(NELL1)、免疫系統(tǒng)過程(SUSD4)、細(xì)胞膜突觸(NECTIN1)等。KEGG分析結(jié)果顯示,這些基因富集到了粘附連接(adherens junction,NECTIN1)、谷氨酸突觸(glutamatergic synapse,GRIK4)、細(xì)胞粘附分子(cell adhesion molecules,NECTIN1)、單純皰疹病毒1型感染(herpes simplex virus 1 infection,NECTIN1)、神經(jīng)活性配體-受體相互作用(neuroactive ligand-receptor interaction,GRIK4)這5條通路(圖9)。

表1 皖岳黑豬特征位點(diǎn)及注釋結(jié)果Table 1 Characteristic loci and annotation results of Wanyue black pig

圖8 皖岳黑豬特定SNPs基因的GO富集結(jié)果Fig.8 GO terms enriched by genes of specific SNPs in Wanyue black pig

圖9 皖岳黑豬特定SNPs基因的KEGG pathway 結(jié)果Fig.9 KEGG pathway enriched by genes of specific SNPs in Wanyue black pig

3 討 論

本研究檢測(cè)并注釋了皖岳黑豬全基因組遺傳變異和SNP位點(diǎn)。在外顯子區(qū)域檢測(cè)到了相對(duì)較多的沉默突變,由于密碼子的簡(jiǎn)并性,沉默突變的發(fā)生對(duì)蛋白質(zhì)的翻譯不造成影響[23],但最近一項(xiàng)重要的研究表明同義突變確實(shí)會(huì)導(dǎo)致蛋白質(zhì)折疊的改變,進(jìn)而損害細(xì)胞功能[24]。本研究中有3.6%的SV發(fā)生在基因組的外顯子區(qū)域,外顯子區(qū)域?qū)儆诨蚓幋a區(qū),稀有且相同的一些SV往往和疾病(包括癌癥)的發(fā)生相互關(guān)聯(lián)甚至還是其直接的致病誘因[25],該區(qū)域內(nèi)的變異對(duì)蛋白質(zhì)的翻譯起重要作用。本研究結(jié)果表明,皖岳黑豬具有豐富的變異信息,這些變異信息的潛在功能對(duì)研究皖岳黑豬種質(zhì)特性具有重要的作用。

當(dāng)前我國絕大多數(shù)畜禽保種場(chǎng)都面臨著群體有效大小偏低的問題,皖岳黑豬的有效群體大小相較于其他品種如北京黑豬(10.1)、民豬(8.1)[26]、撒壩豬(21)[27]要小很多,IBS遺傳距離(0.157 4～0.287 3)也相對(duì)于其他品種小。造成這一結(jié)果的原因一方面是皖岳黑豬屬于培育品種,在培育過程中比較封閉,且存在較嚴(yán)重的近交現(xiàn)象。該結(jié)果提示,在今后的育種過程中需要多加注意,積極引入本品種新的血緣,特別是公豬血緣。多態(tài)性標(biāo)記比例(PN, MAF>0.05)為0.788,表明該品種在培育過程中保留了其親本特異遺傳基因,獲得了較高的遺傳多樣性。

ROH是單個(gè)個(gè)體鑒定的純合基因型片段,由一系列相同的單倍型組成。ROH提供近親繁殖事件的年齡和起源的信息,可用于評(píng)估可靠的近親繁殖系數(shù)[28]。本研究中,皖岳黑豬平均ROH長(zhǎng)度為(61±9.4)Mb,低于青峪豬(131.39±121.33)Mb[29]和安慶六白豬255.19 Mb[30]的平均長(zhǎng)度。當(dāng)前整個(gè)群體的近交系數(shù)平均值為(0.025±0.004),與報(bào)道的恩施黑豬(0.069±0.06)[31]及青峪豬(0.055)[29]的近交系數(shù)相比更低。這可能是因?yàn)橥钤篮谪i在雜交培育過程中血統(tǒng)來源廣泛。

皖岳黑豬經(jīng)過多年的人工選育,集合了來自親本的優(yōu)勢(shì)性狀,具有瘦肉率高,繁殖性能好,抗病能力強(qiáng),環(huán)境適應(yīng)性強(qiáng)等特點(diǎn)。因此,其在基因組水平上一定存在相適應(yīng)的變異,通過全基因組重測(cè)序技術(shù)在基因組水平上查找影響其表型和經(jīng)濟(jì)性狀的候選基因及功能突變位點(diǎn),經(jīng)過統(tǒng)計(jì)分析,篩選其分子遺傳標(biāo)記。選取不同地域分布的豬種,從地理隔絕方面來證實(shí)篩選位點(diǎn)的可靠性,降低基因重復(fù)性,以保證篩選位點(diǎn)最大程度的特異性。對(duì)篩選位點(diǎn)進(jìn)行基因注釋及功能富集分析,發(fā)現(xiàn)這些位點(diǎn)大部分位于內(nèi)含子區(qū)域,一般認(rèn)為內(nèi)含子區(qū)域突變對(duì)基因的功能不產(chǎn)生影響[32],但近幾年的研究表明,內(nèi)含子區(qū)域的變異對(duì)基因的調(diào)控作用及其致病風(fēng)險(xiǎn)要低于CDS區(qū)域和基因調(diào)控區(qū)域[33],但是位于第一個(gè)內(nèi)含子的SNP比其他內(nèi)含子中的SNP致病風(fēng)險(xiǎn)大。內(nèi)含子中SNP主要依靠影響剪切位點(diǎn)活性來影響基因功能。剪切位點(diǎn)的失活可能會(huì)影響翻譯,影響蛋白質(zhì)序列[34],內(nèi)含子中的SNP也有可能影響基因功能。基因功能富集結(jié)果顯示,這些基因在功能上可能與皖岳黑豬生長(zhǎng)速度快、環(huán)境適應(yīng)性強(qiáng)的品種特征相關(guān)。GPC6和NELL1基因與生長(zhǎng)相關(guān)。GPC6 (GLYPICAN 6)是糖基磷脂酰肌醇錨定的硫酸乙酰肝素蛋白聚糖家族的一員,與控制細(xì)胞生長(zhǎng)和細(xì)胞分裂有關(guān)[35-36]。NELL1基因編碼含有表皮生長(zhǎng)因子(EGF)樣重復(fù)序列的細(xì)胞質(zhì)蛋白,是一種對(duì)骨軟骨譜系高度特異性的有效生長(zhǎng)因子,并已證明了對(duì)骨骼的強(qiáng)效誘導(dǎo)[37-38]。NECTIN1編碼粘附蛋白作為偽狂犬病病毒(PRV)的糖蛋白D(gD)的受體,介導(dǎo)病毒進(jìn)入上皮和神經(jīng)元細(xì)胞[39-40]。KHDRBS2基因與大白豬[41]和梅古山羊[42]的繁殖性狀有關(guān)。皖岳黑豬生長(zhǎng)在高海拔、高坡度、高風(fēng)速、高濕度、水質(zhì)優(yōu)良、空氣負(fù)氧離子含量高的生態(tài)環(huán)境,適應(yīng)當(dāng)?shù)囟緡?yán)寒、夏季酷熱的氣候條件,因此形成了其獨(dú)特的種質(zhì)特性,具有較強(qiáng)的環(huán)境適應(yīng)能力。GRIK4基因的剪接和表達(dá)的改變有助于牦牛調(diào)節(jié)其行為認(rèn)知和神經(jīng)系統(tǒng)去適應(yīng)高海拔環(huán)境[43]。CHD1L基因編碼參與DNA修復(fù)的DNA解旋酶蛋白,與癌癥相關(guān)[44]。該基因也被發(fā)現(xiàn)在早期發(fā)育階段與細(xì)胞分裂有關(guān)[45],并通過與PARP1相互作用來調(diào)節(jié)干細(xì)胞的多能性[46],這意味著CHD1L在調(diào)節(jié)干細(xì)胞多能性及胚胎發(fā)育方面具有特定作用[47],可作為皖岳黑豬在胚胎發(fā)育方面研究的候選基因。

4 結(jié) 論

本研究揭示了皖岳黑豬全基因組遺傳變異信息,為今后對(duì)皖岳黑豬表型性狀遺傳機(jī)制的研究奠定了基礎(chǔ)。在后續(xù)選育過程中需要制定更加合理全面的育種方案以確保其遺傳多樣性。基于挑選最大分類能力方法初步獲得皖岳黑豬33個(gè)基因組特異SNPs位點(diǎn),通過注釋獲得這些位點(diǎn)的基因功能,將為未來皖岳黑豬種質(zhì)特性形成的分子機(jī)制研究提供科學(xué)參考。