利用全基因組重測序數(shù)據(jù)檢測8個鴨品種基因組拷貝數(shù)變異

林 燕,黃 敏,李秀金,張續(xù)勐,黃運茂,田允波,伍仲平*

(1.仲愷農(nóng)業(yè)工程學院動物科技學院,廣州 510225;2.浙江農(nóng)林大學動物科技學院·動物醫(yī)學院,杭州 311300)

隨著高通量測序技術的蓬勃發(fā)展,將高深度測序數(shù)據(jù)比對至參考基因組,可獲得大量的變異信息用于畜禽性狀分析。這些變異信息主要包括單核苷酸多態(tài)性(single nucleotide polymorphism,SNP)、插入缺失(insertion-deletion,InDel)以及拷貝數(shù)變異(copy number variation, CNV),其中CNV是常見的結構變異(structural variation, SV)。CNV是指與參考基因組(reference genome)相比,個體基因組中DNA片段從50 bp到5 Mb不等的缺失(deletions)或重復(duplications)[1-2]。由具有重疊片段的相鄰CNV構成的區(qū)域則稱為拷貝數(shù)變異區(qū)(copy number variation region, CNVR)。與單核苷酸多態(tài)性SNP相比,CNV覆蓋的基因組區(qū)域更廣闊,通過劑量效應、位置效應、縮并效應、基因融合或中斷和隱性或功能多態(tài)位點的暴露等作用機制引起染色體結構變異,是導致個體基因和表型變異的重要原因之一[3]。

關于CNV的研究起先主要聚焦在人類的某些遺傳缺陷疾病上,隨后人們在家養(yǎng)動物中也發(fā)現(xiàn)了與表型變異相關的CNV[4-5]。例如,Moller等[6]發(fā)現(xiàn)KIT基因座的串聯(lián)重復導致了豬顯性白毛色的形成;Chen等[7]在豬5號染色體鑒別到38.7 kb大小的CNV通過影響miR-584-5p的表達阻礙了靶基因MSRB3表達從而導致豬耳變大;Zheng等[8]研究發(fā)現(xiàn),AHR基因的拷貝數(shù)變異與豬產(chǎn)仔數(shù)有關;Wright等[9]發(fā)現(xiàn),SOX5基因第一內含子的拷貝數(shù)變異導致雞豆冠表型的形成;Yang等[10]研究發(fā)現(xiàn),HOXB7和HOXB8基因的拷貝數(shù)變異與北京油雞胡須性狀的形成有關;Lin等[11]發(fā)現(xiàn),SOX6基因部分區(qū)域的拷貝數(shù)變異有助于雞的肌肉生長;Weich等[12]研究發(fā)現(xiàn),KITLG基因上游6 kb序列的拷貝數(shù)增加使得家犬的被毛顏色更深、更均勻。

近年來,研究人員利用SNP芯片技術[13-14]、微陣列比較基因組雜交技術(array-based comparative genomic hybridization, aCGH)[15]和全基因組重測序技術(whole-genome sequencing, WGS)已在牛[16-17]、山羊[18-19]、綿羊[20]、豬[21]、狗[22-23]和雞[24-25]等家養(yǎng)動物中開展全基因組范圍內的CNVs檢測。與豬、雞等畜禽相比,鴨基因組CNVs研究相對較少。Skinner等[26]利用aCGH技術在北京鴨中鑒別到32個CNVs,其中5個為雞、火雞等鳥類中共享的保守CNVs;章雙杰等[27]利用重測序數(shù)據(jù)檢測了婁門鴨、昆山麻鴨、巢湖鴨和高郵鴨的全基因組CNVs,發(fā)現(xiàn)婁門鴨和昆山麻鴨中存在1 059個共有的CNVs,并推斷婁門鴨與昆山麻鴨的血緣關系較近。張易[28]在潤州鳳頭白鴨和櫻桃谷鴨雜交F2群體中,利用重測序數(shù)據(jù)構建了鴨全基因組水平的CNV遺傳圖譜,并挖掘到2個與羽色性狀相關的潛在CNV位點;Xu等[29]在北京鴨與野鴨雜交的F2分離群體中,利用全基因組CNVs開展關聯(lián)分析,鑒別到6個與鴨椎骨數(shù)量變異相關的CNVs。

鑒于目前關于多品種鴨全基因組CNVs及品種特異性CNVRs的研究鮮有報道,本研究利用從美國國家生物技術信息中心(National Center for Biotechnology Information, NCBI)公共數(shù)據(jù)庫中下載的包括家鴨和野鴨在內的8個品種共78個個體的全基因組重測序數(shù)據(jù)開展鴨全基因組CNVs檢測,篩選不同鴨品種特有CNVRs,并鑒別與鴨重要經(jīng)濟性狀潛在相關的CNVs,進一步豐富人們對鴨基因組CNVs的了解,為解析CNVs對鴨經(jīng)濟性狀的影響提供前期研究基礎。

1 材料與方法

1.1 全基因組重測序數(shù)據(jù)下載

本研究從NCBI公共數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov)下載了包括櫻桃谷鴨(Cherry Valley duck, CV, n=8)、北京鴨(Beijing duck, BD, n=8)、楓葉鴨(Maple Leaf duck, ML, n=8)、金定鴨(Jinding duck, JD, n=8)、山麻鴨(Shan Partridge duck, SP, n=8)、紹興鴨(Shaoxing duck, SX, n=8)、高郵鴨(Gaoyou duck, GY, n=8)、綠頭野鴨(Mallard, MD, n=22)等8個品種共78個個體的全基因組重測序數(shù)據(jù),所有個體的測序深度均在6.42×(NCBI檢索號:PRJNA419832)[30]。

1.2 全基因組重測序數(shù)據(jù)序列比對

獲得下載的鴨全基因組重測序原始數(shù)據(jù)后,首先使用FastQC軟件(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)對原始數(shù)據(jù)進行質控過濾,然后利用BWA v0.7.17軟件[31]將過濾后的測序數(shù)據(jù)與從Ensembl網(wǎng)站下載的北京鴨參考基因組(CAU_duck1.0)進行比對組裝,最后使用GATK v.4.2.3軟件[32]將組裝好的文件轉換成BAM文件,用于后續(xù)全基因組CNVs檢測。

1.3 CNV檢測和CNVR定義

本研究分別使用CNVnator[33]和CNVcaller[34]軟件進行全基因組CNVs檢測,二者分別采用50 bp和800 bp的窗口大小(bin size)。首先將CNVnator軟件檢測結果中同一品種所有個體存在至少1 bp重疊的同類型CNVs進行合并,得到每個品種的CNVRs,再與CNVcaller軟件檢測結果中存在至少1 bp重疊的同類型CNVRs進行合并,得到每個品種的CNVRs。為了展示CNVRs在基因組上的分布情況,利用R語言ggbio包[35]對鴨全基因組CNVRs遺傳圖譜進行可視化。

1.4 品種特異性CNVRs篩選

本研究將從單個品種中檢測到的CNVRs定義為品種特異性CNVRs,即利用shell語言將上述得到的CNVRs進行品種間比較,篩選在其他品種中不存在的某個品種特有CNVRs,將其作為品種特異性CNVRs。

1.5 CNVRs基因注釋與功能富集分析

本研究利用shell語言將各品種特異性CNVRs在參考基因組中的位置提取出來,通過與北京鴨基因組注釋的GFF文件進行比對,對CNVRs所覆蓋的區(qū)域進行基因注釋。利用PANTHER 17.0軟件(http://www.pantherdb.org)對CNVRs內的基因進行基因本體Gene Ontology(GO)富集分析。

2 結 果

2.1 全基因組CNVRs檢測

本研究分別使用CNVnator和CNVcaller軟件對8個鴨品種共78個個體進行全基因組常染色體CNVs檢測。將兩者結果合并成CNVRs后,在8個鴨品種共檢測到7 550個CNVRs,總長度為16 111.2 kb,其中重復型7 098個,缺失型452個,平均長度為2 134 bp,覆蓋了鴨基因組(常染色體)的1.51%(表1),各品種CNVRs在整個鴨基因組的覆蓋度從0.15%到0.26%不等。在8個品種中,高郵鴨的CNVRs數(shù)量最多,達1 345個,其次為具有1 127個的山麻鴨,而綠頭野鴨具有的CNVRs數(shù)量最少,但具有數(shù)量最多的227個缺失型CNVRs,楓葉鴨中則未檢測到缺失型CNVRs。根據(jù)8個品種的CNVRs檢測結果繪制了鴨常染色體基因組CNVRs分布圖(圖1),可以看出CNVRs在染色體上分布不均勻。

表1 8個鴨品種拷貝數(shù)變異檢測結果Table 1 Descriptive statistics of copy number variant identified in 8 duck breeds

圖1 CNVRs在鴨基因組上的分布Fig.1 Distribution of CNVRs on duck genome

2.2 全基因組CNVRs分布特征

從長度分布情況看,CNVRs主要在1.6～2.1 kb之間,其中有38.07%的CNVRs長度在1.6～1.7 kb之間,41.40%的CNVRs長度在2.0～2.1 kb之間(圖2a)。從不同染色體上的分布情況看,CNVRs在鴨染色體上呈不均勻分布,CNVRs主要集中在1號和2號染色體上,占檢測總數(shù)量的40.53%,而在17號染色體中則未檢測到CNVRs(圖2b)。

圖2 鴨基因組CNVRs的長度和染色體分布Fig.2 CNVRs length and chromosome distribution in duck genome

2.3 品種特異性CNVRs

本研究在8個鴨品種中共篩選到4 304個只存于單個品種中的品種特異性CNVRs,其中重復型4 208個,缺失型96個,總長度為8 412 kb,占檢測CNVRs總長度的52.21%。在8個鴨品種中,高郵鴨的品種特異性CNVRs數(shù)量最多(n=772),其次是櫻桃谷鴨(n=621)。綠頭野鴨的品種特異性CNVRs數(shù)量雖然最少(n=301),但具有最多的品種特異性缺失型CNVRs(n=56)(表2)。

表2 8個鴨品種特異性CNVRsTable 2 The breed-specific CNVRs of 8 duck breeds

2.4 品種特異性CNVRs基因注釋及GO富集分析

本研究將在8個鴨品種中篩選到的4 304個品種特異性CNVRs所覆蓋的區(qū)域進行基因注釋,結果顯示,這些CNVRs共覆蓋了1 230個注釋基因,其中重復型CNVRs包含了1 183個基因,缺失型CNVRs包含了47個基因(圖3)。在8個品種中,櫻桃谷鴨品種特異性CNVRs包含的基因數(shù)量最多(n=305),其次為高郵鴨(n=250),綠頭野鴨最少(n=111)。

圖3 品種特異性CNVRs區(qū)域中的注釋基因數(shù)量Fig.3 Annotation genes in breed-specific CNVRs

為了分析這些基因的生物學功能,本研究利用PANTHER 17.0軟件對篩選到的所有品種特異性CNVRs覆蓋的基因進行了GO富集分析。結果顯示,這些基因主要富集在細胞過程、發(fā)育過程、免疫系統(tǒng)過程、細胞運動、代謝過程、多細胞生物過程、對刺激的反應、信號傳導及生長、繁殖等生物學功能上(圖4)。

圖4 品種特異性CNVRs基因GO富集分析Fig.4 GO enrichment analysis of breed-specific CNVRs genes

為鑒別與品種經(jīng)濟性狀相關特異性CNVRs,本研究分別對8個鴨品種特異性CNVRs覆蓋的基因進行了GO富集分析,共鑒別到38個與繁殖和生長相關的品種特異性CNVRs,其中,在櫻桃谷鴨中鑒別到2個與繁殖相關、1個與生長相關的CNVRs;在北京鴨中鑒別到1個與繁殖相關、1個與生長相關的CNVR;在楓葉鴨中鑒別到3個與生長相關、4個與繁殖相關的CNVRs;在金定鴨中鑒別到3個與繁殖相關、6個與生長相關的CNVRs;在山麻鴨中鑒別到1個與生長相關的CNVR;在紹興鴨中鑒別到1個與生長相關、3個與繁殖相關的CNVRs;在高郵鴨中鑒別到10個與繁殖相關的CNVRs;在綠頭野鴨中鑒別到2個與繁殖相關、1個與生長相關的CNVRs。這些CNVRs共覆蓋了20個與繁殖和生長相關的基因,其中與繁殖相關的基因14個,與生長相關的基因6個(表3)。

表3 8個鴨品種中與生長和繁殖潛在相關的品種特異性CNVRsTable 3 Breed-specific CNVRs potentially related to growth and reproduction in 8 duck breeds

3 討 論

本研究利用CNVnator和CNVcaller軟件,對從公共數(shù)據(jù)庫下載的8個鴨品種共78個個體的全基因組重測序數(shù)據(jù)進行了全基因組拷貝數(shù)變異檢測,并只保留兩個軟件檢測結果中存在至少1 bp重疊的同類型CNVRs,旨在消除假陽性結果對試驗的影響。此外,為了后續(xù)研究方便,本研究在合并CNVRs時,僅考慮只包含重復或缺失片段的CNVRs,未將同時包含重復和缺失片段的混合型CNVRs進行分析。共發(fā)現(xiàn)了7 550個CNVRs,總長度為16 111.2 kb,平均長度為2 134 bp,約占鴨基因組的1.51%,處于馬、豬、牛和雞中報道的0.8%～5.12%范圍內[36-39]。這7 750個CNVRs長度主要集中分布在1.6～2.1 kb之間,且在常染色體上呈不均勻分布,該結果與羊、雞等物種基因組的檢測結果一致[40-41]。

家鴨馴化后經(jīng)過長期的人工選擇,各項生產(chǎn)性能較其祖先綠頭野鴨均有明顯的改變,并根據(jù)人類需要選育出了不同經(jīng)濟類型的品種,如櫻桃谷鴨、北京鴨、楓葉鴨等肉用型品種,金定鴨、山麻鴨、紹興鴨等蛋用型品種以及蛋肉兼用型品種高郵鴨[42]。在家鴨的選育過程中,與經(jīng)濟性狀相關的遺傳變異逐漸在基因組中積累,其中包括拷貝數(shù)變異。因此,本研究為探究品種特異性CNVRs是否與其經(jīng)濟性狀有關,對各品種特異性CNVRs進行了篩選。結果在8個鴨品種中共篩選到4 304個品種特異性CNVRs,包含了1 230個注釋基因。通過對注釋基因開展GO富集分析,發(fā)現(xiàn)多數(shù)基因富集在細胞過程、發(fā)育過程、免疫系統(tǒng)過程等生物學功能上。此外,還有少數(shù)基因與生長和繁殖相關,而這些基因區(qū)域的拷貝數(shù)變異很可能與不同品種鴨特有的經(jīng)濟性狀有關。

為鑒別與不同品種鴨的重要經(jīng)濟性狀潛在相關的CNVRs,分別對各品種特異性CNVRs進行GO富集分析,結果在8個鴨品種中發(fā)現(xiàn)了13個與生長相關的品種特異性CNVRs,其中重復型12個,缺失型1個,這些CNVRs包含了SEMA3E、ANXA6、SEMA3C、ULK1、SLIT2、WWC2等6個基因,其中SLIT2基因已有多個研究報道與肉牛初生重、斷奶重、骨重以及肉雞體重、骨骼生長發(fā)育等生長性狀相關,但未報道與該基因的CNV有關[43-47]。此外,還在8個鴨品種中發(fā)現(xiàn)25個與繁殖相關的重復型CNVRs,包含了BUB1B、SPATA17、MEIOC、TUBGCP5、TUBGCP3、PDE3A、TDRD12、TRIP13、TDRP、TOP2B、FBXW11、PLCZ1、SLX4、TUBGCP6等14個基因。

高郵鴨是我國三大名鴨之一,以善產(chǎn)雙黃蛋而馳名中外。本研究在高郵鴨中篩選到10個品種特異性CNVRs,包含PDE3A、TUBGCP3、TOP2B、TRIP13、SPATA17等5個與繁殖相關的基因。大量研究表明,磷酸二酯酶3A(phosphodiesterase 3A,PDE3A)基因在哺乳動物卵母細胞發(fā)育成熟中起關鍵作用,該基因可通過調控環(huán)磷酸腺苷(cyclic adenosine monophosphate, cAMP)的降解從而促進卵母細胞的減數(shù)分裂并成熟[48-50]。而家禽的雙黃蛋主要是由2個卵泡同時成熟并排卵進入輸卵管而形成的,該過程與卵泡發(fā)育密切相關[51]。因此,推測高郵鴨PDE3A基因區(qū)域的拷貝數(shù)變異很可能與該品種產(chǎn)雙黃蛋的特有經(jīng)濟性狀有關。然而,目前關于PDE3A基因在家禽卵泡發(fā)育中的作用還未見文獻報道,有待進一步的試驗進行驗證。

4 結 論

本研究對NCBI數(shù)據(jù)庫中下載的8個鴨品種共78個個體的重測序數(shù)據(jù)進行全基因組拷貝數(shù)變異檢測,共發(fā)現(xiàn)了7 550個CNVRs,總長度為16 111.2 kb。這些CNVRs在鴨基因組呈不均勻分布,覆蓋了鴨基因組的1.51%。在8個鴨品種全基因組中共篩選4 304個潛在的品種特異性CNVRs,覆蓋1 230個注釋基因。通過基因功能GO富集分析,鑒別到38個可能與鴨生長和繁殖相關的CNVRs。這些結果對于進一步研究CNV與鴨重要經(jīng)濟性狀的關聯(lián)性有重要參考價值。