徐嘉悅,施赫赫,王緒敏,高 翔,余細(xì)勇,于 軍,唐小江
(1.中國(guó)科學(xué)院北京基因研究所,北京 100101;2.廣東省醫(yī)學(xué)實(shí)驗(yàn)動(dòng)物中心,廣東佛山 528248;3.廣東貝格生物科技有限公司,廣東佛山 528100;4.南京大學(xué)模式動(dòng)物研究所,南京 210061;5.廣東省醫(yī)學(xué)科學(xué)院,廣州 510060)
融水小型豬是源自廣西苗寨的小黑豬,具有體型微小、性情溫馴、產(chǎn)仔力強(qiáng)、母性較好等特點(diǎn),具備培育成標(biāo)準(zhǔn)化實(shí)驗(yàn)用小型豬的有利特點(diǎn)。我們于2012年將該豬種源引到廣東三水繁育,建立了基礎(chǔ)種群和的系譜檔案,已繁殖出F1代和F2代,參照北京市地方標(biāo)準(zhǔn)[1]初步建立了融水小型豬的質(zhì)量控制標(biāo)準(zhǔn)。為了進(jìn)一步掌握其遺傳背景,本文對(duì)線粒體DNA進(jìn)行了分析。
融水小型豬由廣東貝格生物科技有限公司(廣東省醫(yī)學(xué)實(shí)驗(yàn)動(dòng)物中心小型豬基地)提供。用無(wú)菌肝素鈉抗凝采血管取融水小型豬外周靜脈血10 mL,用Qiagen公司的QlAamp DNA試劑盒按操作說(shuō)明書(shū)提取DNA。
采用全基因組鳥(niǎo)槍法的測(cè)序策略,使用Illumina基因組分析測(cè)序技術(shù)進(jìn)行測(cè)序研究。為了提高組裝的質(zhì)量,對(duì)融水小型豬選取了插入片段大小為300 bp和500 bp的paired-end文庫(kù)進(jìn)行末端測(cè)序,用基于 Brujin圖算法[2]的 SOAPdenovo軟件(http://soap.genomeics.org.cn)對(duì)測(cè)序得到的短片段進(jìn)行拼接組裝[3-4],得到融水小型豬的完整線粒體DNA序列。
為了確定融水小型豬的線粒體結(jié)構(gòu),首先以Susscrofa的線粒體序列為參考序列,采用BLAST軟件[4]完成了Susscrofa與融水小型豬的比對(duì),確定線粒體的編碼蛋白基因、rRNA和 tRNA,用軟件CodonW對(duì)編碼蛋白基因的密碼子使用情況進(jìn)行統(tǒng)計(jì)。另外用 tRNAscan-SE v.1.21軟件[5]對(duì)線粒體tRNA的二級(jí)結(jié)構(gòu)及各個(gè)tRNA的反義密碼子進(jìn)行預(yù)測(cè)。用clustal W[6]軟件將融水小型豬的編碼蛋白區(qū)序列與GenBank中的15頭國(guó)內(nèi)外的豬線粒體編碼蛋白區(qū)序列(表1)進(jìn)行比對(duì)。
用Dnasp5.10[7]軟件計(jì)算實(shí)驗(yàn)材料的單倍型、遺傳多樣性和遺傳分化的相關(guān)數(shù)值[8]?;谝陨系膯伪缎途垲?lèi)的分類(lèi),采用 Tajima’s D[9]和 Fu’s[10]兩個(gè)中性檢驗(yàn)?zāi)P蛯?duì)16個(gè)豬的單倍型分布進(jìn)行中性檢驗(yàn),其中種群遺傳格局及變異的檢測(cè)在Arkequin3.11[11]中 實(shí) 施。 使 用 分 子 變 異 分 析(AMOVA)的方法以10000次重復(fù)隨機(jī)抽樣單倍型重排后進(jìn)行顯著性檢驗(yàn),以此來(lái)評(píng)價(jià)種群遺傳變異水平與種群地理位置的相關(guān)性,使用統(tǒng)計(jì)量將遺傳差異分為3個(gè)不同層次的地理等級(jí)。用Network 4.6(http://www.fluxus-engineering.com/sharenet.htm)對(duì)16個(gè)豬種進(jìn)行聚類(lèi),并用MEGA 5.0[12]的最大似然法和鄰接法對(duì)16個(gè)豬種的編碼蛋白區(qū)進(jìn)行系統(tǒng)進(jìn)化分析。
融水小型豬與其他豬的結(jié)構(gòu)框架相同,都是由13個(gè)編碼蛋白,22個(gè)tRNA和2個(gè)rRNA組成(表2)。線粒體全長(zhǎng)為16888 bp。線粒體DNA的核苷酸組成為A>C>T>G。線粒體的GC含量總是低于AT含量。
融水小型豬由13個(gè)基因構(gòu)成,除基因ND6位于線粒體的輕鏈上之外,其余12個(gè)基因包括ND1、ND2、ND3、ND4、ND4L、ND5、COX1、COX2、COX3、CYTB、ATP6、ATP8都在線粒體的重鏈(表2)。這13個(gè)基因的總長(zhǎng)為11255 bp,其中基因ATP6和基因ATP8有42 bp的重疊區(qū)域,基因ND4L和基因ND4有7 bp的重疊區(qū)域,基因ND5和基因ND6有17 bp的重疊區(qū)域。位于融水小型豬重鏈的12個(gè)編碼蛋白基因的堿基組分為A>C >T>G,其中ND1、ATP6、COX3、ND3、ND4、ND5、ND2,CYTB8 個(gè)基因的堿基組份是一致的:A >C >T >G,而基因 COX1、COX2、ATP8、ND4L的堿基組份為A>T>C>G;另外,位于融水小型豬輕鏈的ND6基因的堿基組份為C>A>T>G。
在拼接的融水小型豬線粒體中,D-loop區(qū)的長(zhǎng)度為1254 bp,tRNA-Pro和tRNA-Phe排布于D-loop區(qū)的兩端。在D-loop區(qū)堿基組份為A>C>T>G,與線粒體DNA一致,AT的含量高于GC含量。此外,融水小型豬線粒體DNA中非編碼的基因間區(qū)長(zhǎng)度為654 bp,這些基因間區(qū)的主要原因是線粒體中tRNA簇的重排。其中最長(zhǎng)的一段基因間區(qū)是位于ND1和tRNA之間的區(qū)域,長(zhǎng)為139 bp。基因間區(qū)為高級(jí)別的系統(tǒng)進(jìn)化研究提供了可能。
表1 15個(gè)豬種線粒體DNA的基本信息及其在NCBI數(shù)據(jù)庫(kù)中的檢索號(hào)Tab.1 List of mitochondrial genome sequences of Swine breeds and their accession numbers as retrieved from NCBI database
表2 融水小型豬線粒體DNA的結(jié)構(gòu)特征Tab.2 Location of feature in the mitochondrial DNA of Rongshui miniature pig
線粒體中包括兩個(gè)核糖體 RNA,分別為12S rRNA and 16S rRNA(表1),12S rRNA的長(zhǎng)度為959 bp,16S rRNA的長(zhǎng)度為1512 bp。16S rRNA編碼于纈氨酸和亮氨酸之間,其堿基組份為A>T>C>G;12S rRNA編碼于丙氨酸和纈氨酸之間,其堿基組份為A>C>T>G。12S rRNA和16S rRNA的GC含量均在40%左右,同樣保持了GC含量低于AT含量這一特點(diǎn)。
對(duì)融水小型豬線粒體,我們預(yù)測(cè)了其22個(gè)轉(zhuǎn)運(yùn)RNA的二級(jí)結(jié)構(gòu),并確定了各轉(zhuǎn)運(yùn)RNA的反義密碼子(如表2所示)。這22個(gè)轉(zhuǎn)運(yùn)RNA的GC含量都低于AT含量,其中精氨酸的GC含量最低,約21%;甲硫氨酸的GC含量最高,約47%。
融水小型豬線粒體DNA的密碼子使用情況(表3)。包括終止密碼子在內(nèi),共有3772個(gè)密碼子,其中使用密碼子CUA的數(shù)目最多,為306個(gè),而AGG和CGG兩個(gè)密碼子沒(méi)有被使用。在密碼子的統(tǒng)計(jì)中,密碼子第一位核苷酸使用的頻率為A>C>T>G;而在密碼子第二位置,堿基使用頻率為T(mén)>C>A>G;在密碼子的第三個(gè)位置,堿基使用頻率為A>C>U>G。眾所周知,密碼子在不同位置使用頻率的差異是由于密碼子使用的偏好性造成的,上述的密碼子使用情況也為研究線粒體DNA的進(jìn)化提供了有用的信息[13]。
同義突變和非同義突變之間的比例(Ka/Ks值)可用于分析不同物種之間的的系統(tǒng)進(jìn)化關(guān)系[14]。通過(guò)對(duì)16個(gè)豬13個(gè)編碼蛋白區(qū)的所有序列進(jìn)行比對(duì),構(gòu)建了一個(gè)參考序列(Ref),分別計(jì)算每頭豬的各個(gè)編碼蛋白基因序列與Ref的Ka/Ks比值。結(jié)果如封2圖1所示,16個(gè)豬種在編碼蛋白區(qū)的Ka/KS值均小于1,說(shuō)明所有13個(gè)編碼蛋白基因都在承受著正選擇,其中基因COX1的選擇壓力最大,這意味著該基因正在承受著很強(qiáng)的純化選擇;而ND5的選擇壓力相對(duì)較小,說(shuō)明該基因是在進(jìn)化中舍棄一些較為不利的突變。
使用 MEGA軟件中的 Kimura’s two-parameter模型(K2P)計(jì)算16個(gè)豬種的編碼蛋白區(qū)的遺傳距離。融水小型豬與其他豬的遺傳距離最遠(yuǎn),堿基的平均遺傳距離為1.141。而對(duì)總體的遺傳距離統(tǒng)計(jì)中,物種間以恒定速率進(jìn)行核苷酸替換的標(biāo)準(zhǔn)差僅為0.0015?;谏鲜鲞z傳距離矩陣,通過(guò)MEGA軟件中的最大似然法和鄰接法對(duì)16頭豬的編碼蛋白區(qū)構(gòu)建系統(tǒng)進(jìn)化樹(shù)(圖2),從圖2可以看出,16個(gè)豬種可以分成四個(gè)類(lèi)群,即歐洲群(E群)、亞洲1群(A1群)、亞洲2群(A2群)、亞洲3群(A3群)。其中,E群主要是由歐洲豬種組成,其中意大利的野豬是相對(duì)古老的豬種;A1群是由中國(guó)、日本、韓國(guó)等亞洲豬種及Berkshire豬共同組成,在這一分支中,也可以證明豬的進(jìn)化是由亞洲起源后分化到歐洲的[15];A2群是由滇南7號(hào)、海南野豬和蘭嶼豬構(gòu)成,這一結(jié)論也驗(yàn)證了之前的研究,即云南豬種和臺(tái)灣地區(qū)蘭嶼豬是古老的豬種[15]。在系統(tǒng)進(jìn)化分析中最值得關(guān)注的是融水小型豬,可以看出他比蘭嶼豬更為古老。
系統(tǒng)進(jìn)化分析是對(duì)過(guò)去已經(jīng)發(fā)生的進(jìn)化時(shí)間的模擬分析,為了驗(yàn)證系統(tǒng)進(jìn)化分析結(jié)果的準(zhǔn)確性,用單倍型的方法,對(duì)16個(gè)豬種的進(jìn)化關(guān)系做了進(jìn)一步分析。用Dnasp5.10[7]分析16個(gè)豬種的線粒體DNA的編碼區(qū)序列,共得到14個(gè)單倍型,但單倍型的多樣性指數(shù)為0.596。將分析后的單倍型數(shù)據(jù)利用Network軟件作圖(封2圖3)。14個(gè)單倍型可以分成3大類(lèi):漢普郡豬、蘇塞克斯家養(yǎng)豬、長(zhǎng)白豬、蘇塞克斯豬、杜洛克豬、意大利野豬、歐洲伯克希爾豬、聚為一類(lèi)(類(lèi)群1);濟(jì)州8#豬、濟(jì)州10#豬、韓國(guó)野豬、梅山豬和臺(tái)灣地區(qū)野豬聚為一類(lèi)(類(lèi)群2);蘭嶼豬、滇南7號(hào)豬、海南野豬、融水小型豬聚為一類(lèi)(類(lèi)群3;其中類(lèi)群1主要是歐洲豬種,類(lèi)群2主要以亞洲的韓國(guó)和臺(tái)灣地區(qū)豬種為主,類(lèi)群3主要是中國(guó)4頭古老的豬種??梢钥闯鯪etwork作圖結(jié)果與系統(tǒng)進(jìn)化分析的結(jié)果大致是一致的,這也進(jìn)一步驗(yàn)證了融水小型豬是一個(gè)非常古老的豬種。
基于單倍型的分析結(jié)果,分別對(duì)3個(gè)類(lèi)群進(jìn)行了中性檢驗(yàn)(如表4)。類(lèi)群3不但具有較高的單倍型多樣性(h=0.899),也具有較高的核苷酸多樣性(Pi=0.03417),表明類(lèi)群3中的豬種仍具有豐富的遺傳多樣性,存在著較高的進(jìn)化潛力。對(duì)于類(lèi)群1和類(lèi)群2,單倍型多樣性較高,但是核苷酸多樣性很低。說(shuō)明這兩個(gè)類(lèi)群的豬種經(jīng)歷過(guò)近期擴(kuò)張或者進(jìn)化歷史很短[17]。本研究中類(lèi)群1和類(lèi)群3的Tajin a’s D與Fs值為負(fù)值且差異顯著,這表明類(lèi)群1和類(lèi)群3中的豬種在歷史上都發(fā)生了種群擴(kuò)張。
表3 融水小型豬密碼子Tab.3 Codon usage for Rongshui miniature pig
圖2 16個(gè)豬種的系統(tǒng)進(jìn)化分析(最大似然法和鄰接法)Fig.2 The phylogenetic analysis of 16 pig breeds(NJ,ML analyses)
表4 對(duì)16個(gè)豬種編碼蛋白基因的中性檢驗(yàn)和遺傳多樣性分析Tab.4 Neutrality test and genetic diversity for main haplogroups in coding region for 16 swine mitochondrial genomes
線粒體DNA存在于細(xì)胞質(zhì)線粒體機(jī)質(zhì)中的閉合環(huán)狀雙鏈超螺旋DNA分子,是獨(dú)立于細(xì)胞核染色體外的基因組,具有自我復(fù)制、轉(zhuǎn)錄和編碼功能,是母系遺傳,但同時(shí)受核 DNA的控制[2]。與核酸DNA相比,線粒體DNA具有分子結(jié)構(gòu)簡(jiǎn)單,以母性方式遺傳,核苷酸歧義大,進(jìn)化速度快等特點(diǎn)。因此,使用線粒體研究哺乳動(dòng)物群體遺傳多樣性和種內(nèi)、種間親緣關(guān)系等,一直受到很多學(xué)者的關(guān)注[15,18-20]。本研究采用新一代測(cè)序技術(shù),對(duì)融水小型豬進(jìn)行測(cè)序分析,拼接成完整的線粒體DNA序列,從分子遺傳學(xué)角度評(píng)估了融水小型豬的結(jié)構(gòu)質(zhì)特性、系統(tǒng)進(jìn)化關(guān)系及其遺傳多樣性。為融水小型豬的實(shí)驗(yàn)動(dòng)物化提供科學(xué)依據(jù)。
隨著新一代測(cè)序技術(shù)的發(fā)展,越來(lái)越多的學(xué)者對(duì)哺乳動(dòng)物遺傳機(jī)制的研究不再僅僅停留在核DNA水平上,而是更多的綜合細(xì)胞質(zhì)DNA進(jìn)行分析。線粒體DNA結(jié)構(gòu)簡(jiǎn)單、穩(wěn)定,是母系遺傳方式,在遺傳過(guò)程中未曾發(fā)生重組。
本文使用新一代測(cè)序技術(shù),用Illumina測(cè)序平臺(tái)對(duì)融水小型豬進(jìn)行從頭測(cè)序,使用SOAPdenovo拼接成完整線粒體序列并進(jìn)行注釋。融水小型豬的線粒體DNA成環(huán)形,總長(zhǎng)為16888 bp,由13個(gè)基因,22個(gè)tRNA和2個(gè)rRNA組成。其中GC的平均含量約為44%。線粒體的GC含量總是低于AT含量,這一結(jié)論與其他物種一致,如人的GC含量為44.4%[21],小鼠的 GC 含量為 26.6%[22]。本文結(jié)合公共數(shù)據(jù)庫(kù)發(fā)布的15個(gè)豬種線粒體序列,以線粒體的編碼蛋白區(qū)為研究對(duì)象,對(duì)豬種的遺傳多樣性進(jìn)行了研究,進(jìn)而對(duì)豬種的起源和進(jìn)化進(jìn)行探討。從系統(tǒng)進(jìn)化分析中我們可以得出,融水小型豬在系統(tǒng)進(jìn)化樹(shù)中成為一個(gè)新的獨(dú)立外群,我們認(rèn)為融水小型豬是一個(gè)古老于蘭嶼豬的傳統(tǒng)中國(guó)豬種。為了進(jìn)一步驗(yàn)證這個(gè)結(jié)論,我們對(duì)系統(tǒng)進(jìn)化樹(shù)中的三個(gè)類(lèi)群進(jìn)行了遺傳多樣性分析,Network的結(jié)果與系統(tǒng)進(jìn)化分析的結(jié)果是吻合的。由滇南7號(hào)、海南野豬、蘭嶼豬和融水小型豬構(gòu)成的類(lèi)群3中有4個(gè)單倍型,其具有較高的單倍型多樣性和核苷酸多樣性顯示出類(lèi)群3具有豐富的遺傳多樣性。
Tajin a‘s D與Fu’s Fs中性檢驗(yàn)常用以推測(cè)種群擴(kuò)張的歷史事件。如果Tajin a’s D與Fs值呈負(fù)值,且在統(tǒng)計(jì)學(xué)上達(dá)到顯著標(biāo)準(zhǔn),則說(shuō)明序列中含有比中性進(jìn)化模型中更多的核苷位點(diǎn)變化,可能預(yù)示著該種群曾經(jīng)有過(guò)擴(kuò)張的歷史。融水小型豬的Tajin a’s D與Fs值為負(fù)值且差異顯著,說(shuō)明了類(lèi)群3的古老豬種在長(zhǎng)期進(jìn)化過(guò)程中還在與其他的豬種還在發(fā)生基因交流和種群擴(kuò)張。
另外,在系統(tǒng)進(jìn)化分析時(shí),蘭嶼豬、滇南7號(hào)豬、海南野豬與融水小型豬分成了不同的分支,而在Network分析中成為一個(gè)類(lèi)群,出現(xiàn)這種現(xiàn)象的原因可能是融水小型豬在長(zhǎng)期的進(jìn)化過(guò)程中,發(fā)生了基因交流。中性檢驗(yàn)的結(jié)果不但從統(tǒng)計(jì)水平上驗(yàn)證了融水小型豬在豬線粒體發(fā)展中的祖先地位,而且還解釋了融水小型豬作為祖先豬種是如何與其他豬種發(fā)生基因流和種群擴(kuò)張的。遺傳多樣性分析進(jìn)一步提示融水小型豬是中國(guó)傳統(tǒng)豬種中一種古老的小型豬。
本文對(duì)融水小型豬DNA進(jìn)行了初步的分析,限于篇幅和GenBank的數(shù)據(jù),未與哥廷根小型豬、西藏小型豬、巴馬小型豬等其他實(shí)驗(yàn)用小型豬的線粒體DNA進(jìn)行比較,有待進(jìn)一步研究。
[1] 北京市質(zhì)量技術(shù)監(jiān)督局.實(shí)驗(yàn)用小型豬:微生物學(xué)等級(jí)及監(jiān)測(cè)DB11T 828.1-2011[S];寄生蟲(chóng)學(xué)等級(jí)及監(jiān)測(cè) DB11T 828.2-2011[S];遺傳質(zhì)量控制 DB11T 828.3 -2011[S];病理學(xué)診斷規(guī)范 DB11T 828.4-2011[S];配合飼料 DB11T 828.5-2011[S];環(huán)境及設(shè)施 DB11T828.6 -2011[S],2011-11-10發(fā)布,2012-03-01實(shí)施.
[2] Pevzner,P.A.,H.Tang,M.S.Waterman,An Eulerian path approach to DNA fragment assembly[J].Proc Natl Acad Sci U S A,2001,98(17):74-85.
[3] Angiero and Natalie.,Pigs Prove to Be Smart,if Not Vain.The New York Times,2009.
[4] Altschul,S.F.,Gish,W.,Miller,W.,et al.,Basic local alignment search tool[J].J Mol Biol,1990,215(3):403-410.
[5] Lowe,T.M.S.R.Eddy,tRNAscan-SE:a program for improved detection of transfer RNA genes in genomic sequence[J].Nucleic Acids Res,1997,25(5):95 -98.
[6] Thompson,J.D.,D.G.Higgins,T.J.Gibson,CLUSTAL W:improving the sensitivity ofprogressive multiple sequence alignment through sequence weighting,position-specific gap penalties and weight matrix choice[J].Nucleic Acids Res,1994,22(22):73-80.
[7] Rozas,J.,Sanchez-DelBarrio,J.C.,Messeguer,X.,et al.,DnaSP,DNA polymorphism analyses by the coalescent and other methods[J].Bioinformatics,2003.19(18):49 -67.
[8] Kimura,M.,Molecular evolutionary clock and the neutral theory[J].J Mol Evol,1987.26(1-2):24-33.
[9] Tajima,F(xiàn).,Statistical method for testing the neutral mutation hypothesis by DNA polymorphism[J].Genetics,1989,123(3):58-63.
[10] Fu,Y.X.,Statistical tests of neutrality of mutations against population growth,hitchhiking and background selection[J].Genetics,1997,147(2):91 -105.
[11] Excoffier L,L.L., Schneider S, Arlequin ver.3.0:an integrated software package for population genetics data analysis[J].Evolutionary Bioinformatics Online,2005,1:47 -50.
[12] Tamura,K.,Peterson,D.,Peterson,N.,et al.,MEGA5:molecular evolutionary genetics analysis using maximum likelihood, evolutionarydistance, and maximum parsimony methods[J].Mol Biol Evol.28(10):73 -91.
[13] Dubey,B., MeganathanP.R., HaqueI., DNA minibarcoding:an approach for forensic identification of some endangered Indian snake species[J].Forensic Sci Int Genet,5(3):18-24.
[14] Yang, Z. R. Nielsen, Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models[J].Mol Biol Evol,2000,17(1):32 -43.
[15] Chen,Huang.,H.L,Yang H.Y.,et al.Mitochondrial genome of Taiwan pig (Susscrofa) [J]. African Journalof Biotechnology,2011,10(13):2556-2561.
[16] Wu,G.S.,Yao,Y.G.,Qu,K.X.,et al.,Population phylogenomic analysis of mitochondrial DNA in wild boars and domestic pigs revealed multiple domestication events in East Asia[J].Genome Biol,2007,8(11):245 -253.
[17] Femandez,A.I.,Alves,E.,Ovilo,C.,et al.,Divergence time estimates of East Asian and European pigs based on multiple near complete mitochondrial DNA sequences[J].Anim Genet,42(1):86-88.
[18] Kijas,J.M.Andersson L.,A phylogenetic study of the origin of the domestic pig estimated from the near-complete mtDNA genome[J].J Mol Evol,2001,52(3):302 -8.
[19] Femandez,A.I.,Alves,E.,Ovilo,C.,et al.,Divergence time estimates of East Asian and European pigs based on multiple near complete mitochondrial DNA sequences[J].Anim Genet,42(1):86-89.
[20] Yang,S.,Zhang,H.,Mao,H.,et al.,The local origin of the Tibetan pig and additional insights into the origin of Asian pigs[J].PLoS One,6(12):218 -215.
[21] Anderson,S.,Bankier,A.T.,Barrell,B.G.,et al.,Sequence and organization of the human mitochondrial genome[J].Nature,1981,290(5806):45 -65.
[22] Bibb,M.J.,Van Etten,R.A.,Wright,C.T.,et al.,Sequence and gene organization of mouse mitochondrial DNA[J].Cell,1981,26(2 Pt 2):67 -80.