林 燕,黃 敏,李秀金,張續(xù)勐,黃運(yùn)茂,田允波,伍仲平*
(1.仲愷農(nóng)業(yè)工程學(xué)院動物科技學(xué)院,廣州 510225;2.浙江農(nóng)林大學(xué)動物科技學(xué)院·動物醫(yī)學(xué)院,杭州 311300)
隨著高通量測序技術(shù)的蓬勃發(fā)展,將高深度測序數(shù)據(jù)比對至參考基因組,可獲得大量的變異信息用于畜禽性狀分析。這些變異信息主要包括單核苷酸多態(tài)性(single nucleotide polymorphism,SNP)、插入缺失(insertion-deletion,InDel)以及拷貝數(shù)變異(copy number variation, CNV),其中CNV是常見的結(jié)構(gòu)變異(structural variation, SV)。CNV是指與參考基因組(reference genome)相比,個體基因組中DNA片段從50 bp到5 Mb不等的缺失(deletions)或重復(fù)(duplications)[1-2]。由具有重疊片段的相鄰CNV構(gòu)成的區(qū)域則稱為拷貝數(shù)變異區(qū)(copy number variation region, CNVR)。與單核苷酸多態(tài)性SNP相比,CNV覆蓋的基因組區(qū)域更廣闊,通過劑量效應(yīng)、位置效應(yīng)、縮并效應(yīng)、基因融合或中斷和隱性或功能多態(tài)位點(diǎn)的暴露等作用機(jī)制引起染色體結(jié)構(gòu)變異,是導(dǎo)致個體基因和表型變異的重要原因之一[3]。
關(guān)于CNV的研究起先主要聚焦在人類的某些遺傳缺陷疾病上,隨后人們在家養(yǎng)動物中也發(fā)現(xiàn)了與表型變異相關(guān)的CNV[4-5]。例如,Moller等[6]發(fā)現(xiàn)KIT基因座的串聯(lián)重復(fù)導(dǎo)致了豬顯性白毛色的形成;Chen等[7]在豬5號染色體鑒別到38.7 kb大小的CNV通過影響miR-584-5p的表達(dá)阻礙了靶基因MSRB3表達(dá)從而導(dǎo)致豬耳變大;Zheng等[8]研究發(fā)現(xiàn),AHR基因的拷貝數(shù)變異與豬產(chǎn)仔數(shù)有關(guān);Wright等[9]發(fā)現(xiàn),SOX5基因第一內(nèi)含子的拷貝數(shù)變異導(dǎo)致雞豆冠表型的形成;Yang等[10]研究發(fā)現(xiàn),HOXB7和HOXB8基因的拷貝數(shù)變異與北京油雞胡須性狀的形成有關(guān);Lin等[11]發(fā)現(xiàn),SOX6基因部分區(qū)域的拷貝數(shù)變異有助于雞的肌肉生長;Weich等[12]研究發(fā)現(xiàn),KITLG基因上游6 kb序列的拷貝數(shù)增加使得家犬的被毛顏色更深、更均勻。
近年來,研究人員利用SNP芯片技術(shù)[13-14]、微陣列比較基因組雜交技術(shù)(array-based comparative genomic hybridization, aCGH)[15]和全基因組重測序技術(shù)(whole-genome sequencing, WGS)已在牛[16-17]、山羊[18-19]、綿羊[20]、豬[21]、狗[22-23]和雞[24-25]等家養(yǎng)動物中開展全基因組范圍內(nèi)的CNVs檢測。與豬、雞等畜禽相比,鴨基因組CNVs研究相對較少。Skinner等[26]利用aCGH技術(shù)在北京鴨中鑒別到32個CNVs,其中5個為雞、火雞等鳥類中共享的保守CNVs;章雙杰等[27]利用重測序數(shù)據(jù)檢測了婁門鴨、昆山麻鴨、巢湖鴨和高郵鴨的全基因組CNVs,發(fā)現(xiàn)婁門鴨和昆山麻鴨中存在1 059個共有的CNVs,并推斷婁門鴨與昆山麻鴨的血緣關(guān)系較近。張易[28]在潤州鳳頭白鴨和櫻桃谷鴨雜交F2群體中,利用重測序數(shù)據(jù)構(gòu)建了鴨全基因組水平的CNV遺傳圖譜,并挖掘到2個與羽色性狀相關(guān)的潛在CNV位點(diǎn);Xu等[29]在北京鴨與野鴨雜交的F2分離群體中,利用全基因組CNVs開展關(guān)聯(lián)分析,鑒別到6個與鴨椎骨數(shù)量變異相關(guān)的CNVs。
鑒于目前關(guān)于多品種鴨全基因組CNVs及品種特異性CNVRs的研究鮮有報(bào)道,本研究利用從美國國家生物技術(shù)信息中心(National Center for Biotechnology Information, NCBI)公共數(shù)據(jù)庫中下載的包括家鴨和野鴨在內(nèi)的8個品種共78個個體的全基因組重測序數(shù)據(jù)開展鴨全基因組CNVs檢測,篩選不同鴨品種特有CNVRs,并鑒別與鴨重要經(jīng)濟(jì)性狀潛在相關(guān)的CNVs,進(jìn)一步豐富人們對鴨基因組CNVs的了解,為解析CNVs對鴨經(jīng)濟(jì)性狀的影響提供前期研究基礎(chǔ)。
本研究從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]。
獲得下載的鴨全基因組重測序原始數(shù)據(jù)后,首先使用FastQC軟件(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)對原始數(shù)據(jù)進(jìn)行質(zhì)控過濾,然后利用BWA v0.7.17軟件[31]將過濾后的測序數(shù)據(jù)與從Ensembl網(wǎng)站下載的北京鴨參考基因組(CAU_duck1.0)進(jìn)行比對組裝,最后使用GATK v.4.2.3軟件[32]將組裝好的文件轉(zhuǎn)換成BAM文件,用于后續(xù)全基因組CNVs檢測。
本研究分別使用CNVnator[33]和CNVcaller[34]軟件進(jìn)行全基因組CNVs檢測,二者分別采用50 bp和800 bp的窗口大小(bin size)。首先將CNVnator軟件檢測結(jié)果中同一品種所有個體存在至少1 bp重疊的同類型CNVs進(jìn)行合并,得到每個品種的CNVRs,再與CNVcaller軟件檢測結(jié)果中存在至少1 bp重疊的同類型CNVRs進(jìn)行合并,得到每個品種的CNVRs。為了展示CNVRs在基因組上的分布情況,利用R語言ggbio包[35]對鴨全基因組CNVRs遺傳圖譜進(jìn)行可視化。
本研究將從單個品種中檢測到的CNVRs定義為品種特異性CNVRs,即利用shell語言將上述得到的CNVRs進(jìn)行品種間比較,篩選在其他品種中不存在的某個品種特有CNVRs,將其作為品種特異性CNVRs。
本研究利用shell語言將各品種特異性CNVRs在參考基因組中的位置提取出來,通過與北京鴨基因組注釋的GFF文件進(jìn)行比對,對CNVRs所覆蓋的區(qū)域進(jìn)行基因注釋。利用PANTHER 17.0軟件(http://www.pantherdb.org)對CNVRs內(nèi)的基因進(jìn)行基因本體Gene Ontology(GO)富集分析。
本研究分別使用CNVnator和CNVcaller軟件對8個鴨品種共78個個體進(jìn)行全基因組常染色體CNVs檢測。將兩者結(jié)果合并成CNVRs后,在8個鴨品種共檢測到7 550個CNVRs,總長度為16 111.2 kb,其中重復(fù)型7 098個,缺失型452個,平均長度為2 134 bp,覆蓋了鴨基因組(常染色體)的1.51%(表1),各品種CNVRs在整個鴨基因組的覆蓋度從0.15%到0.26%不等。在8個品種中,高郵鴨的CNVRs數(shù)量最多,達(dá)1 345個,其次為具有1 127個的山麻鴨,而綠頭野鴨具有的CNVRs數(shù)量最少,但具有數(shù)量最多的227個缺失型CNVRs,楓葉鴨中則未檢測到缺失型CNVRs。根據(jù)8個品種的CNVRs檢測結(jié)果繪制了鴨常染色體基因組CNVRs分布圖(圖1),可以看出CNVRs在染色體上分布不均勻。
表1 8個鴨品種拷貝數(shù)變異檢測結(jié)果Table 1 Descriptive statistics of copy number variant identified in 8 duck breeds
圖1 CNVRs在鴨基因組上的分布Fig.1 Distribution of CNVRs on duck genome
從長度分布情況看,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
本研究在8個鴨品種中共篩選到4 304個只存于單個品種中的品種特異性CNVRs,其中重復(fù)型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
本研究將在8個鴨品種中篩選到的4 304個品種特異性CNVRs所覆蓋的區(qū)域進(jìn)行基因注釋,結(jié)果顯示,這些CNVRs共覆蓋了1 230個注釋基因,其中重復(fù)型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
為了分析這些基因的生物學(xué)功能,本研究利用PANTHER 17.0軟件對篩選到的所有品種特異性CNVRs覆蓋的基因進(jìn)行了GO富集分析。結(jié)果顯示,這些基因主要富集在細(xì)胞過程、發(fā)育過程、免疫系統(tǒng)過程、細(xì)胞運(yùn)動、代謝過程、多細(xì)胞生物過程、對刺激的反應(yīng)、信號傳導(dǎo)及生長、繁殖等生物學(xué)功能上(圖4)。
圖4 品種特異性CNVRs基因GO富集分析Fig.4 GO enrichment analysis of breed-specific CNVRs genes
為鑒別與品種經(jīng)濟(jì)性狀相關(guān)特異性CNVRs,本研究分別對8個鴨品種特異性CNVRs覆蓋的基因進(jìn)行了GO富集分析,共鑒別到38個與繁殖和生長相關(guān)的品種特異性CNVRs,其中,在櫻桃谷鴨中鑒別到2個與繁殖相關(guān)、1個與生長相關(guān)的CNVRs;在北京鴨中鑒別到1個與繁殖相關(guān)、1個與生長相關(guān)的CNVR;在楓葉鴨中鑒別到3個與生長相關(guān)、4個與繁殖相關(guān)的CNVRs;在金定鴨中鑒別到3個與繁殖相關(guān)、6個與生長相關(guān)的CNVRs;在山麻鴨中鑒別到1個與生長相關(guān)的CNVR;在紹興鴨中鑒別到1個與生長相關(guān)、3個與繁殖相關(guān)的CNVRs;在高郵鴨中鑒別到10個與繁殖相關(guān)的CNVRs;在綠頭野鴨中鑒別到2個與繁殖相關(guān)、1個與生長相關(guān)的CNVRs。這些CNVRs共覆蓋了20個與繁殖和生長相關(guān)的基因,其中與繁殖相關(guān)的基因14個,與生長相關(guān)的基因6個(表3)。
表3 8個鴨品種中與生長和繁殖潛在相關(guān)的品種特異性CNVRsTable 3 Breed-specific CNVRs potentially related to growth and reproduction in 8 duck breeds
本研究利用CNVnator和CNVcaller軟件,對從公共數(shù)據(jù)庫下載的8個鴨品種共78個個體的全基因組重測序數(shù)據(jù)進(jìn)行了全基因組拷貝數(shù)變異檢測,并只保留兩個軟件檢測結(jié)果中存在至少1 bp重疊的同類型CNVRs,旨在消除假陽性結(jié)果對試驗(yàn)的影響。此外,為了后續(xù)研究方便,本研究在合并CNVRs時,僅考慮只包含重復(fù)或缺失片段的CNVRs,未將同時包含重復(fù)和缺失片段的混合型CNVRs進(jìn)行分析。共發(fā)現(xiàn)了7 550個CNVRs,總長度為16 111.2 kb,平均長度為2 134 bp,約占鴨基因組的1.51%,處于馬、豬、牛和雞中報(bào)道的0.8%~5.12%范圍內(nèi)[36-39]。這7 750個CNVRs長度主要集中分布在1.6~2.1 kb之間,且在常染色體上呈不均勻分布,該結(jié)果與羊、雞等物種基因組的檢測結(jié)果一致[40-41]。
家鴨馴化后經(jīng)過長期的人工選擇,各項(xiàng)生產(chǎn)性能較其祖先綠頭野鴨均有明顯的改變,并根據(jù)人類需要選育出了不同經(jīng)濟(jì)類型的品種,如櫻桃谷鴨、北京鴨、楓葉鴨等肉用型品種,金定鴨、山麻鴨、紹興鴨等蛋用型品種以及蛋肉兼用型品種高郵鴨[42]。在家鴨的選育過程中,與經(jīng)濟(jì)性狀相關(guān)的遺傳變異逐漸在基因組中積累,其中包括拷貝數(shù)變異。因此,本研究為探究品種特異性CNVRs是否與其經(jīng)濟(jì)性狀有關(guān),對各品種特異性CNVRs進(jìn)行了篩選。結(jié)果在8個鴨品種中共篩選到4 304個品種特異性CNVRs,包含了1 230個注釋基因。通過對注釋基因開展GO富集分析,發(fā)現(xiàn)多數(shù)基因富集在細(xì)胞過程、發(fā)育過程、免疫系統(tǒng)過程等生物學(xué)功能上。此外,還有少數(shù)基因與生長和繁殖相關(guān),而這些基因區(qū)域的拷貝數(shù)變異很可能與不同品種鴨特有的經(jīng)濟(jì)性狀有關(guān)。
為鑒別與不同品種鴨的重要經(jīng)濟(jì)性狀潛在相關(guān)的CNVRs,分別對各品種特異性CNVRs進(jìn)行GO富集分析,結(jié)果在8個鴨品種中發(fā)現(xiàn)了13個與生長相關(guān)的品種特異性CNVRs,其中重復(fù)型12個,缺失型1個,這些CNVRs包含了SEMA3E、ANXA6、SEMA3C、ULK1、SLIT2、WWC2等6個基因,其中SLIT2基因已有多個研究報(bào)道與肉牛初生重、斷奶重、骨重以及肉雞體重、骨骼生長發(fā)育等生長性狀相關(guān),但未報(bào)道與該基因的CNV有關(guān)[43-47]。此外,還在8個鴨品種中發(fā)現(xiàn)25個與繁殖相關(guān)的重復(fù)型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個與繁殖相關(guān)的基因。大量研究表明,磷酸二酯酶3A(phosphodiesterase 3A,PDE3A)基因在哺乳動物卵母細(xì)胞發(fā)育成熟中起關(guān)鍵作用,該基因可通過調(diào)控環(huán)磷酸腺苷(cyclic adenosine monophosphate, cAMP)的降解從而促進(jìn)卵母細(xì)胞的減數(shù)分裂并成熟[48-50]。而家禽的雙黃蛋主要是由2個卵泡同時成熟并排卵進(jìn)入輸卵管而形成的,該過程與卵泡發(fā)育密切相關(guān)[51]。因此,推測高郵鴨PDE3A基因區(qū)域的拷貝數(shù)變異很可能與該品種產(chǎn)雙黃蛋的特有經(jīng)濟(jì)性狀有關(guān)。然而,目前關(guān)于PDE3A基因在家禽卵泡發(fā)育中的作用還未見文獻(xiàn)報(bào)道,有待進(jìn)一步的試驗(yàn)進(jìn)行驗(yàn)證。
本研究對NCBI數(shù)據(jù)庫中下載的8個鴨品種共78個個體的重測序數(shù)據(jù)進(jìn)行全基因組拷貝數(shù)變異檢測,共發(fā)現(xiàn)了7 550個CNVRs,總長度為16 111.2 kb。這些CNVRs在鴨基因組呈不均勻分布,覆蓋了鴨基因組的1.51%。在8個鴨品種全基因組中共篩選4 304個潛在的品種特異性CNVRs,覆蓋1 230個注釋基因。通過基因功能GO富集分析,鑒別到38個可能與鴨生長和繁殖相關(guān)的CNVRs。這些結(jié)果對于進(jìn)一步研究CNV與鴨重要經(jīng)濟(jì)性狀的關(guān)聯(lián)性有重要參考價(jià)值。