楊 晴,鞏 靜,趙雪艷,朱曉東,耿立英,張傳生*,王繼英*

(1.河北科技師范學(xué)院動(dòng)物科技學(xué)院,秦皇島 066600;2.山東省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所 山東省畜禽疫病防治與繁育重點(diǎn)實(shí)驗(yàn)室,濟(jì)南 250100;3.農(nóng)業(yè)農(nóng)村部畜禽生物組學(xué)重點(diǎn)實(shí)驗(yàn)室,濟(jì)南 250100;4.棗莊黑蓋豬養(yǎng)殖有限公司,棗莊 277100)

隨著高密度芯片和測(cè)序技術(shù)的高速發(fā)展,生物全基因組范圍內(nèi)檢測(cè)出的標(biāo)記數(shù)量逐漸增多,大量的遺傳標(biāo)記信息使基因組遺傳變異分析得以更為準(zhǔn)確和精準(zhǔn)的實(shí)施。當(dāng)前,SNP芯片與測(cè)序技術(shù)已成為動(dòng)植物進(jìn)行遺傳變異信息分析工作的主要工具,被廣泛應(yīng)用于遺傳多樣性分析[1-2]、選擇信號(hào)檢測(cè)[3-4]、全基因組關(guān)聯(lián)分析[5-6]、基因組選擇等[7-8]。SNP芯片具有自動(dòng)化、成本低、效率高等優(yōu)點(diǎn),但其存在檢測(cè)位點(diǎn)較少、特異性強(qiáng)、無(wú)法發(fā)現(xiàn)新功能位點(diǎn)等缺點(diǎn)[9]。測(cè)序不受參考基因組的限制,甚至可通過(guò)提高測(cè)序深度來(lái)獲得所測(cè)樣本的全部遺傳變異信息,包括覆蓋低、中、高密度甚至全基因組范圍內(nèi)的所有已知或未知的SNP位點(diǎn)信息[10],但其存在數(shù)據(jù)量大、分析復(fù)雜、成本較高的缺點(diǎn)[11]。實(shí)際研究中,還是要根據(jù)研究目的選擇適當(dāng)?shù)姆中头椒?。雖然增加SNP檢測(cè)密度會(huì)提高分析結(jié)果的準(zhǔn)確度,但在實(shí)際應(yīng)用中,高密度分型會(huì)帶來(lái)高昂的經(jīng)濟(jì)成本,極大地限制了高密度標(biāo)記在基因組遺傳變異分析中的應(yīng)用空間,所以根據(jù)研究目的探索適當(dāng)?shù)臉?biāo)記密度和經(jīng)濟(jì)的分型方法,保證分析結(jié)果的準(zhǔn)確性,成為近年來(lái)SNP標(biāo)記分析的熱點(diǎn)研究?jī)?nèi)容之一。

已有報(bào)道顯示,低密度面板的基因組選擇,通過(guò)基因型填充等方法能夠達(dá)到中高密度,甚至測(cè)序數(shù)據(jù)相似的基因組預(yù)測(cè)精確程度,是一種低成本且高效的遺傳評(píng)估方法[12-14]。但是分析不同分型方法或不同SNP密度對(duì)全基因組遺傳變異分析結(jié)果是否存在影響且影響是否較大的報(bào)道仍較少。因此,本研究以35頭棗莊黑蓋豬的高密度SNP芯片數(shù)據(jù)和重測(cè)序SNP數(shù)據(jù)為基礎(chǔ),利用重測(cè)序信息構(gòu)建不同密度的SNP面板,以探究不同SNP分型方法和不同SNP密度對(duì)遺傳變異分析的影響,找到適用于遺傳變異分析的低成本、高效的分型方法和SNP密度,為今后豬及其他畜禽遺傳特性分析中適宜的基因分型技術(shù)和標(biāo)記密度的選擇提供重要參考。

1 材料與方法

1.1 樣本來(lái)源

本研究所用的35頭棗莊黑蓋豬均采自山東省棗莊黑蓋豬養(yǎng)殖有限公司,包括16頭母豬和19頭公豬。采集試驗(yàn)豬耳組織樣品存放于裝有75%酒精的2 mL凍存管內(nèi),放入-20 ℃低溫冰箱中保存?zhèn)溆谩?/p>

1.2 DNA提取與質(zhì)檢

取樣本耳組織0.5 g左右,采用血液/細(xì)胞/組織基因組DNA提取試劑盒(DP304,TIANGEN公司,北京)進(jìn)行基因組DNA的提取。利用NanoDrop 2000和瓊脂糖凝膠電泳對(duì)DNA的濃度和質(zhì)量進(jìn)行檢測(cè),濃度>50 ng·μL-1,1.8

1.3 基因分型和質(zhì)控

使用CAUPorcineSNP50芯片(北京康普森生物技術(shù)有限公司)對(duì)35個(gè)個(gè)體進(jìn)行SNP分型,SNP檢出率平均為97.97%?；谌A大-MGISEQ-T7技術(shù)測(cè)序平臺(tái),利用雙末端測(cè)序(paired-end)的方法對(duì)35個(gè)個(gè)體進(jìn)行基因組重測(cè)序,平均測(cè)序深度為13X,Q20為98.18%。原始數(shù)據(jù)質(zhì)控后,使用BWA軟件[15]的BWA-MEN算法將質(zhì)控?cái)?shù)據(jù)與參考基因組(Ensembl Sus Scrofa11.1)進(jìn)行比對(duì),使用GATA[16]進(jìn)行重比對(duì),最后使用Samtools軟件[17]和Bcftools軟件[18]檢測(cè)基因組范圍內(nèi)的SNP。

使用Plink(V1.90)[19]對(duì)SNP芯片和重測(cè)序中的數(shù)據(jù)按如下標(biāo)準(zhǔn)進(jìn)行質(zhì)量控制,標(biāo)準(zhǔn)如下:1)僅保留位于常染色體上的SNP位點(diǎn);2)芯片數(shù)據(jù)刪除檢出率(call rate)<90%的SNP位點(diǎn),重測(cè)序數(shù)據(jù)刪除檢出率<95%的SNP位點(diǎn);3)刪除檢出率<90%的個(gè)體;4)刪除最小等位基因頻率(MAF)<0.05的SNP位點(diǎn)。

1.4 不同密度SNP面板的構(gòu)建

基于重測(cè)序檢測(cè)的SNP位點(diǎn),利用R語(yǔ)言CVrepGPAcalc包(https://github.com/SmaragdaT/CVrep/)構(gòu)建不同密度的SNP面板[20],依據(jù)SNP芯片密度共設(shè)計(jì)了3個(gè)梯度,分別為34K、340K和3 400K。面板的構(gòu)建有兩種方法,第一種是在整個(gè)基因組中隨機(jī)抽樣來(lái)選擇SNP,第二種是根據(jù)特定步長(zhǎng)的物理距離均勻的選擇SNP。其中,34K面板選擇兩種方法分別進(jìn)行構(gòu)建,340K和3 400K面板均采用第二種方法進(jìn)行構(gòu)建。

1.5 數(shù)據(jù)分析

使用Plink(V1.90)計(jì)算群體的最小等位基因頻率(minor allele frequency, MAF)、觀察雜合度(observed heterozygosity,HO)、期望雜合度(expected heterozygosity,HE)、群體內(nèi)遺傳距離等遺傳多樣性指標(biāo),使用Plink(V1.90)將數(shù)據(jù)格式轉(zhuǎn)化為vcf格式,再利用vcf2phylip和Phylip通過(guò)鄰接法(neighbor-joining, NJ)構(gòu)建系統(tǒng)發(fā)生樹(shù)[21-22],最后利用FigTreev1.4.4軟件(http://tree.bio.ed.ac.uk/software/figtree/)將計(jì)算結(jié)果可視化。使用Plink(V1.90)計(jì)算狀態(tài)同源距離(identity by descent distance, IBS距離),隨后計(jì)算個(gè)體間遺傳距離(1-IBS距離),并利用BioLadder在線軟件(https://www.bioladder.cn/web/#/chart/6)繪制個(gè)體間遺傳距離熱圖。

使用R語(yǔ)言CMplot軟件包對(duì)SNP在染色體上的分布進(jìn)行可視化,使用R語(yǔ)言detectRUNS軟件包[23]對(duì)基因組進(jìn)行長(zhǎng)純合片段(runs of homogeneity,ROH)檢測(cè)并計(jì)算各分組內(nèi)的群體內(nèi)近交系數(shù)(FROH),參數(shù)設(shè)置[24-26]為:SNP密度最小為每1 000 kb必須有1個(gè)SNP;連續(xù)兩個(gè)SNPs的間隔最大為1 000 kb;滑窗大小為50個(gè)SNPs;ROH滑窗中允許有1個(gè)SNP位點(diǎn)為雜合;ROH滑窗中允許有5個(gè)SNPs位點(diǎn)缺失;滑動(dòng)窗口重疊比例至少為5%;ROH最少個(gè)數(shù)為40個(gè)SNPs。

2 結(jié) 果

2.1 SNP分型與質(zhì)控

利用CAUPorcineSNP50 芯片和基因組重測(cè)序?qū)?5頭棗莊黑蓋豬進(jìn)行基因組SNP檢測(cè),分別獲得了43 832個(gè)和31 437 418個(gè)SNPs位點(diǎn)。芯片的SNP檢出率平均為0.979 8,重測(cè)序的檢出率平均為0.997 0。各質(zhì)控條件下芯片和重測(cè)序數(shù)據(jù)SNP位點(diǎn)的詳細(xì)剔除數(shù)量見(jiàn)表1。經(jīng)過(guò)數(shù)據(jù)質(zhì)控后,芯片和測(cè)序數(shù)據(jù)剩余位點(diǎn)的比例分別為78.69%和65.76%。

表1 SNP質(zhì)控結(jié)果匯總Table 1 Summary of SNP quality control results

2.2 不同密度SNP面板的構(gòu)建

通過(guò)質(zhì)控標(biāo)準(zhǔn)的芯片SNP位點(diǎn)個(gè)數(shù)為34 494個(gè)。依據(jù)芯片密度(34K)設(shè)置梯度,以重測(cè)序數(shù)據(jù)為“原材料”構(gòu)建不同密度SNP面板。芯片和各密度SNP面板的SNP位點(diǎn)數(shù)目、MAF和相鄰SNP間距詳見(jiàn)表2?？梢钥闯?芯片標(biāo)記MAF均值為0.292,高于測(cè)序各組標(biāo)記的MAF均值(0.244～0.245)。密度同為34K的3組相比,芯片SNP間距均值最大(70 809.82 bp),均勻34K的次之(65 819.90 bp),隨機(jī)34K的最小(63 359.20 bp)。但是,隨機(jī)34K組SNP間距的標(biāo)準(zhǔn)差最大(80 185.61 bp),遠(yuǎn)高于芯片(57 626.21 bp)和均勻34K(1 771.16 bp)。綜合來(lái)看,芯片的SNP位點(diǎn)在染色體上的分布均勻度介于隨機(jī)34K和均勻34K之間。與圖1密度分布圖所示結(jié)果一致。不同密度測(cè)序SNP面板(均勻34K、均勻340K和均勻3 400K)相比較,均勻34K的SNP間距均值約為均勻340K的10倍,基本與構(gòu)建面板時(shí)采用的步長(zhǎng)大小(10×)相一致,標(biāo)準(zhǔn)差大小隨SNP密度的增加而減小。

表2 芯片和各測(cè)序面板SNP數(shù)目、最小等位基因頻率和間距Table 2 SNP number, MAF and space of adjacent SNPs of array and sequencing panels

2.3 遺傳多樣性和遺傳距離分析

利用芯片和各測(cè)序SNP面板的SNP標(biāo)記分析棗莊黑蓋豬的遺傳多樣性結(jié)果見(jiàn)表3?？梢钥闯?利用芯片SNP標(biāo)記分析的HO、HE、遺傳距離均高于測(cè)序各組,利用各測(cè)序面板SNP標(biāo)記分析的HO、HE、遺傳距離基本相同,特別是均勻分布的3組SNP(34K、340K和3400K)的遺傳多樣性指標(biāo)更為接近。圖2展示了使用芯片和測(cè)序各組數(shù)據(jù)分析的35頭棗莊黑蓋豬樣本間遺傳距離矩陣熱圖,與表3結(jié)果一致,芯片與隨機(jī)34K及均勻分布SNP組間的差別最為明顯。

A. 芯片;B. 隨機(jī)34K;C. 均勻34K。矩陣中每一個(gè)小方格代表樣本兩兩之間的遺傳距離值,該值越大越接近紫色,越小越接近黃綠色A. Array; B. Random 34K; C. Even 34K. Each small square in the matrix represents the genetic distance value between two samples, the larger the value, the color is closer to purple, and the smaller the value, the color is closer to yellow-green圖2 樣本間遺傳距離熱圖Fig.2 Heat map of genetic distance between samples

表3 芯片和各測(cè)序面板遺傳多樣性參數(shù)值Table 3 Values of genetic diversity analyzed based on array and sequencing panels

2.4 系統(tǒng)發(fā)生樹(shù)

利用芯片和各測(cè)序SNP面板的SNP標(biāo)記構(gòu)建了棗莊黑蓋豬群體鄰接法系統(tǒng)發(fā)生樹(shù),詳見(jiàn)圖3。系統(tǒng)發(fā)生樹(shù)是表示個(gè)體間親緣關(guān)系的樹(shù)狀圖,相同分支上的個(gè)體具有相近親緣關(guān)系,為同一個(gè)家系。可以看出,基于芯片和各測(cè)序SNP面板的SNP標(biāo)記構(gòu)建的系統(tǒng)發(fā)生樹(shù)均將35頭棗莊黑蓋豬劃分為3大分支,每個(gè)大分支又可進(jìn)一步細(xì)分成1～3個(gè)小分支。仔細(xì)對(duì)比分支上的個(gè)體,芯片與隨機(jī)34K、芯片與3組均勻SNP數(shù)據(jù)均存在一定的差別,而3個(gè)均勻分布的SNP數(shù)據(jù)(34K、340K和3 400K)構(gòu)建的系統(tǒng)發(fā)生樹(shù)基本一致。

A. 芯片;B. 隨機(jī)34K;C. 均勻34K;D. 均勻340KA. Array; B. Random 34K; C. Even 34K; D. Even 340K圖3 鄰接法構(gòu)建的系統(tǒng)發(fā)生樹(shù)Fig.3 Phylogenetic trees constructed by neighbor-joining method

2.5 基于ROH的基因組近交系數(shù)分析

利用芯片和各測(cè)序SNP面板的SNP標(biāo)記分析了棗莊黑蓋豬ROH和基因組近交系數(shù),詳見(jiàn)表4?？梢钥闯?芯片與隨機(jī)34K相比,芯片檢測(cè)的ROH數(shù)目少(723vs. 784),但ROH長(zhǎng)度大(14.86 Mbvs. 12.85 Mb),二者的FROH相近(0.125vs.0.127);均勻34K與隨機(jī)34K相比,均勻34K數(shù)據(jù)檢測(cè)到ROH數(shù)目更多(789vs. 784),長(zhǎng)度更大(13.51 Mbvs. 12.85 Mb),FROH近交系數(shù)更高(0.134vs. 0.127)。3個(gè)均勻分布的數(shù)據(jù)組相比,隨著標(biāo)記密度增加,檢測(cè)的ROH數(shù)目逐漸增多,ROH長(zhǎng)度逐漸降低,估計(jì)的FROH近交系數(shù)也逐漸增加。

表4 芯片和各測(cè)序面板ROH及基因組近交系數(shù)值Table 4 ROH and genomic inbreeding coefficients based on array and sequencing panels

3 討 論

單核苷酸多態(tài)性(SNPs)是人類(lèi)和其他動(dòng)物可遺傳的變異中最常見(jiàn)的一種,在基因組中廣泛存在,作為第三代分子標(biāo)記在畜禽遺傳多樣性分析、選擇信號(hào)檢測(cè)、全基因組關(guān)聯(lián)分析、基因組選擇等方面發(fā)揮著重要作用。近來(lái)的研究表明,人類(lèi)基因組上SNP總數(shù)可達(dá)3 800萬(wàn)個(gè)[27],目前已鑒定出的豬SNP已經(jīng)超過(guò)四千余萬(wàn)個(gè)[28]?；蚪M測(cè)序可以獲得所測(cè)樣本的全部SNP信息,因此,WGS數(shù)據(jù)有望可以用來(lái)更好地估計(jì)個(gè)體之間的真實(shí)關(guān)系[29]。SNP芯片僅包含了鑒定出的SNP位點(diǎn)的一個(gè)子集,SNP芯片的覆蓋率和密度適當(dāng)?shù)那闆r下,在估計(jì)基因組關(guān)系、遺傳多樣性分析等方面與測(cè)序技術(shù)一樣有價(jià)值[30]。

本研究中,35頭棗莊黑蓋豬基因組重測(cè)序共檢測(cè)到3 143.7萬(wàn)個(gè)SNPs位點(diǎn),經(jīng)過(guò)數(shù)據(jù)質(zhì)控后,測(cè)序數(shù)據(jù)剩余位點(diǎn)的比例(65.76%)小于芯片數(shù)據(jù)(78.69%),這與基因組重測(cè)序檢測(cè)到的SNPs中含有大量的(6 976 769個(gè),占位點(diǎn)總數(shù)的22.19%)低MAF位點(diǎn)(MAF<0.05)有關(guān)。與本研究結(jié)果一致,Wang 等[31]、Eynard等[29]在對(duì)大約克豬、荷斯坦牛的基因組測(cè)序數(shù)據(jù)分析中也發(fā)現(xiàn)基因組測(cè)序包含了20%左右的低MAF(MAF<0.05)位點(diǎn)。與基因組測(cè)序相比,芯片基因組SNP在設(shè)計(jì)過(guò)程中,優(yōu)先選擇測(cè)序樣本中發(fā)現(xiàn)的高M(jìn)AF的SNP位點(diǎn)[32]。本研究所用的CAUPorcineSNP50 芯片整合現(xiàn)有重要經(jīng)濟(jì)功能基因公開(kāi)報(bào)道的候選位點(diǎn),并加入部分地方豬種全基因組重測(cè)序鑒定的特有SNP綜合優(yōu)化研制而成,所以該芯片SNP位點(diǎn)平均MAF值(0.292)高于各測(cè)序面板(0.244～0.245)。

利用芯片SNP標(biāo)記分析的HO、HE、遺傳距離等遺傳多樣性各指標(biāo)值均高于測(cè)序各組,利用芯片SNP標(biāo)記構(gòu)建的系統(tǒng)發(fā)生樹(shù)與測(cè)序各組也存在較大不同,而測(cè)序各組SNP標(biāo)記分析的遺傳多樣性各指標(biāo)值基本相同,構(gòu)建的系統(tǒng)發(fā)生樹(shù)基本相似。本研究結(jié)果說(shuō)明,分型方法對(duì)遺傳多樣性、遺傳距離和系統(tǒng)發(fā)生樹(shù)分析存在影響。以往的研究表明,芯片SNP位點(diǎn)由于傾向于選擇高M(jìn)AF位點(diǎn)、位點(diǎn)群體代表性不全面等原因(即確定偏倚(ascertainment biases))會(huì)影響遺傳多樣性、群體分化、連鎖不平衡等分析的結(jié)果[33-35]。據(jù)此推測(cè),本研究中芯片與測(cè)序?qū)z傳距離分析結(jié)果的不同可能是由于芯片和測(cè)序標(biāo)記MAF差異所致。在測(cè)序方法下,不同SNP密度對(duì)遺傳多樣性、遺傳距離和系統(tǒng)發(fā)生樹(shù)分析結(jié)果影響較小,說(shuō)明3.4萬(wàn)個(gè)標(biāo)記已經(jīng)能充分滿足系統(tǒng)發(fā)生樹(shù)分析所需的標(biāo)記數(shù)量,增加標(biāo)記數(shù)目和增加數(shù)據(jù)運(yùn)算量并不能進(jìn)一步提高遺傳多樣性和系統(tǒng)發(fā)生樹(shù)的分析精確性。

4 結(jié) 論

本研究以重測(cè)序數(shù)據(jù)為“原材料”構(gòu)建了不同密度SNP面板,利用芯片和各測(cè)序SNP面板的SNPs標(biāo)記分析棗莊黑蓋豬的遺傳多樣性、系統(tǒng)發(fā)生樹(shù)和基因組近交系數(shù)。結(jié)果表明,利用芯片SNP標(biāo)記分析的HO、HE、遺傳距離等遺傳多樣性指標(biāo)值均高于各測(cè)序組,利用芯片SNP標(biāo)記構(gòu)建的系統(tǒng)發(fā)生樹(shù)與各測(cè)序組也存在較大不同,此外,芯片數(shù)據(jù)檢測(cè)出的ROH長(zhǎng)度較測(cè)序組大,基于ROH計(jì)算的近交系數(shù)偏小。各測(cè)序組的不同SNP密度對(duì)遺傳多樣性和系統(tǒng)發(fā)生樹(shù)分析結(jié)果影響較小,但對(duì)ROH及基于ROH計(jì)算的基因組近交系數(shù)影響很大。因此,在研究初期進(jìn)行試驗(yàn)設(shè)計(jì)時(shí),要根據(jù)研究目的選擇適宜的基因分型技術(shù)和標(biāo)記密度,以降低成本和提高結(jié)果的準(zhǔn)確性。