雞繁殖性能近交衰退相關(guān)CpG島差異甲基化基因的篩選

薛 倩，李國輝，殷建玫，張會(huì)永，周成浩，朱云芬，邢偉杰，蘇一軍，鄒劍敏，韓 威*

(1.江蘇省家禽科學(xué)研究所，揚(yáng)州 225125; 2.江蘇省家禽科學(xué)研究所 科技創(chuàng)新中心，揚(yáng)州 225125)

由于近交而引起的后代群體適合度下降的現(xiàn)象被稱為近交衰退[1]，近交衰退是物種資源保護(hù)尤其是瀕危物種保護(hù)過程中備受關(guān)注的問題之一。群體遺傳學(xué)提出了兩個(gè)解釋近交衰退遺傳基礎(chǔ)的經(jīng)典假說: 顯性假說和超顯性假說。顯性假說認(rèn)為，近交過程中，基因組純合性的增加導(dǎo)致了有害隱性等位基因的暴露；超顯性假說提出，在整個(gè)基因組中，雜合子一般優(yōu)于純合子，近交增加了純合子出現(xiàn)的概率，因此會(huì)表現(xiàn)出衰退[2]。近些年研究發(fā)現(xiàn)，這兩種假說對(duì)于某些近交衰退現(xiàn)象卻無法完全解釋[2-6]，比如近交衰退程度會(huì)隨著環(huán)境壓力的變化而變化的不確定現(xiàn)象等[7-8]。DNA表觀修飾也會(huì)受到環(huán)境的影響，進(jìn)而影響到基因的表達(dá)和表型的可塑性[9-11]。并且已有研究表明，近交衰退與DNA表觀遺傳修飾密切相關(guān)[12-14]。表觀遺傳變化是發(fā)生在基因組DNA中的重要變異，它在不直接改變DNA編碼序列的情況下影響著基因的表達(dá)[15-16]。

DNA甲基化是目前研究最廣泛的表觀遺傳修飾方式之一[17]，其參與多種生物學(xué)過程[18-19]，在調(diào)控基因表達(dá)和維持基因組穩(wěn)定性方面起著關(guān)鍵的作用[20]。脊椎動(dòng)物基因組DNA甲基化主要發(fā)生在 CpG胞嘧啶上，并且它的調(diào)控多為通過 CpG 島(CpG island)多個(gè) CpG 位點(diǎn)的甲基化共同作用來完成，而不是通過對(duì)某單個(gè) CpG 位點(diǎn)的甲基化或去甲基化來激活或失活基因[21]。與CpG 島鄰近(～2 kb) 的低CpG 密度區(qū)域被定義為CpG 島岸(CpG island shore)，這些區(qū)域的甲基化與轉(zhuǎn)錄的激活也密切相關(guān)?？锓f等[22]研究表明，近交系小鼠生長發(fā)育緩慢可能與H19基因上游區(qū)域被甲基化進(jìn)而引起該基因的異常表達(dá)有關(guān)；蔣曹德[23]以梅山豬、大白豬群體為親本構(gòu)建了兩個(gè)自繁群體與兩個(gè)雜交群體，發(fā)現(xiàn)群體內(nèi)自繁的子代甲基化水平略高于雜交子代群體，且兩個(gè)群體間某些性狀差異可能與某特定位點(diǎn)甲基化差異相關(guān)[23]；Venney等[24]研究表明，三文魚的近交衰退與其CK-1、GTIIBS和hsp70三個(gè)特定基因的甲基化變化有關(guān)；楊晶淼等[25]研究報(bào)告了馬氏珠母貝的近交與雜交家系的生長性狀和甲基化水平均存在顯著性差異。植物上，關(guān)于山蘿卜[14]和茄科植物[26]的兩項(xiàng)研究均表明，DNA甲基化與近交衰退之間存在著密切的關(guān)系，并且近交個(gè)體比遠(yuǎn)交個(gè)體存在更為廣泛的全基因組DNA甲基化現(xiàn)象。盡管有這些研究，但是人們對(duì)DNA甲基化，尤其是脊椎動(dòng)物CpG 島及其附近區(qū)域甲基化在近交衰退中的調(diào)控作用還沒有完全深刻地理解。

家禽繁殖性能在家禽生產(chǎn)、育種和品種保護(hù)中具有極其重要的地位，尤其在地方雞遺傳資源活體保護(hù)過程中，群體繁殖性能近交衰退現(xiàn)象時(shí)有發(fā)生，嚴(yán)重影響著保種群的規(guī)模和世代的延續(xù)。本課題組在前期的研究中，已成功組建了狼山雞高近交組(分子近交系數(shù)FIS> 0.15)和低近交組(分子近交系數(shù)FIS< 0.04)，且觀察到狼山雞高近交組繁殖性能發(fā)生了顯著衰退現(xiàn)象[27]，然而，關(guān)于繁殖性能近交衰退個(gè)體DNA甲基化變化，尤其是CpG 島及其附近重要調(diào)控區(qū)域甲基化變化情況尚未進(jìn)行深入探索。基于此，本研究利用新一代測序技術(shù)—全基因組重亞硫酸鹽測序(whole-genome bisulfite sequencing，WGBS)技術(shù)檢測和分析狼山雞高、低近交組性腺軸組織(包括卵巢和下丘腦)全基因組甲基化水平變化情況，著重篩選CpG 島及其附近區(qū)域內(nèi)存在的差異甲基化區(qū)域(differential methylation region，DMR)，并探討其在狼山雞繁殖性能近交衰退中可能發(fā)揮的調(diào)控作用，研究結(jié)果將為今后家禽育種和物種資源保護(hù)工作提供理論參考依據(jù)。

1 材料與方法

1.1 試驗(yàn)動(dòng)物和樣品采集

選取江蘇省家禽科學(xué)研究所國家級(jí)地方雞種基因庫保存的地方雞品種—狼山雞為試驗(yàn)素材，基于前期狼山雞高、低近交組母雞繁殖性狀(包括開產(chǎn)日齡、300天產(chǎn)蛋數(shù)、開產(chǎn)體重和開產(chǎn)蛋重)的表型值記錄，選取高近交組中繁殖性能顯著衰退的個(gè)體3只 和低近交組中繁殖性能處于群體均值附近的個(gè)體3只，快速屠宰并采集其下丘腦和卵巢組織，迅速投入液氮中，然后轉(zhuǎn)移至-80 ℃冰箱中保存?zhèn)溆谩?/p>

1.2 DNA提取和WGBS建庫測序

使用德國凱杰公司提供的組織DNA快速提取試劑盒提取各樣本基因組DNA,并使用安捷倫2100分光光度計(jì)進(jìn)行檢測。DNA樣品經(jīng)超聲片段化后進(jìn)行亞硫酸氫鹽轉(zhuǎn)化，使用DNA甲基化測序文庫構(gòu)建試劑盒(斯威夫特生物科技，美國)，將單鏈DNA片段進(jìn)行接頭連接，對(duì)DNA片段進(jìn)行8個(gè)循環(huán)PCR擴(kuò)增：95 ℃預(yù)變性2 min；8個(gè)循環(huán)(95 ℃變性30 s，65 ℃退火20 s，72 ℃延伸45 s)；最后72 ℃ 延伸性7 min。擴(kuò)增反應(yīng)體系包括預(yù)混反應(yīng)體系25 μL (低濃度EDTA TE 10 μL, Buffer R1 10 μL, 反應(yīng)物R2 4 μL，R3酶1 μL)，DNA純化產(chǎn)物和探針引物25 μL，共計(jì)50 μL。對(duì)擴(kuò)增產(chǎn)物進(jìn)行純化和完整性檢測，最后使用Illumina HiSeq 4000 測序平臺(tái)對(duì)測序文庫進(jìn)行雙末端測序，每個(gè)文庫測序深度達(dá)30×。

1.3 數(shù)據(jù)過濾質(zhì)控和甲基化胞嘧啶(mC)位點(diǎn)識(shí)別

使用Cutadapt[28]軟件去除測序Reads中含有測序接頭、污染、低質(zhì)量堿基和未知堿基的Reads，序列質(zhì)量進(jìn)一步使用FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)進(jìn)行驗(yàn)證。通過質(zhì)控的Reads被比對(duì)到雞的參考基因組中，使用Samtool[29]軟件對(duì)比對(duì)上的Reads進(jìn)一步去重復(fù)。對(duì)于每個(gè)胞嘧啶位點(diǎn)，其甲基化水平取決于覆蓋到mC位點(diǎn)的Reads數(shù)與總Reads數(shù)的比值。使用二項(xiàng)分布檢驗(yàn)鑒別錯(cuò)誤檢出率(false-discovery rate，F(xiàn)DR)≤0.01的mC位點(diǎn)，所有mC位點(diǎn)測序覆蓋度需≥3。根據(jù)mC位點(diǎn)上、下游堿基種類，將mC位點(diǎn)分為3種類型：mCG、mCHG和mCHH，其中H代表腺苷(A)、胞嘧啶(C)或胸腺嘧啶(T)。

1.4 DMR檢測與CpG島區(qū)DMR相關(guān)基因功能富集分析

使用R語言包MethylKit[30]進(jìn)行DMR檢測，采用滑窗方法(1 000 bp窗口，500 bp步長)掃描基因組，檢測高、低近交組中存在差異甲基化的區(qū)域，對(duì)每個(gè)窗口進(jìn)行Fisher精確檢驗(yàn)，并使用FDR對(duì)多次檢驗(yàn)的P值進(jìn)行校正。FDR<0.05且甲基化水平差異倍數(shù)至少為2倍的窗口被鑒定作為DMR。采用emboss[31-32]進(jìn)行CpG 島鑒定，將CpG島定義為長度至少100 bp、GC含量大于50%，且CpG二核苷酸的出現(xiàn)率(觀測值/期望值)大于60%的基因組區(qū)域。對(duì)CpG 島區(qū)(CpG 島及其上、下游2 kb區(qū)域)甲基化信息進(jìn)行統(tǒng)計(jì)分析，篩選位于CpG 島區(qū)的DMR。為探索這些DMR的功能，對(duì)包含這些DMR的基因(稱為DMR相關(guān)基因或差異甲基化基因)進(jìn)行GO和KEGG數(shù)據(jù)庫功能注釋和富集分析，P<0.05的GO條目和Pathway信號(hào)通路為顯著富集的條目和通路。

1.5 亞硫酸氫鹽測序PCR(bisulfite sequencing PCR,BSP)方法驗(yàn)證WGBS分析結(jié)果

為驗(yàn)證WGBS測序分析結(jié)果的可靠性，從卵巢中隨機(jī)選取2個(gè)DMR相關(guān)基因(CDC27、SRD5A1)，對(duì)其CpG島區(qū)域甲基化水平進(jìn)行驗(yàn)證。利用Methylation Primer Express V1.0軟件設(shè)計(jì)擴(kuò)增引物，引物信息見表1。選用WGBS相同的DNA樣品，每個(gè)樣品取1 mg，用EZ DNA Methylation-GoldTM試劑盒進(jìn)行亞硫酸氫鹽轉(zhuǎn)化。利用設(shè)計(jì)的引物擴(kuò)增轉(zhuǎn)化后的DNA目標(biāo)區(qū)域，PCR產(chǎn)物純化后進(jìn)行TA載體克隆，每個(gè)基因挑選10～15個(gè)克隆子進(jìn)行Sanger測序，對(duì)測序結(jié)果采用BiQ Analyzer軟件分析。

表1 BSP驗(yàn)證試驗(yàn)引物信息

1.6 數(shù)據(jù)處理分析

采用IBM SPSS Statistics 20.0軟件統(tǒng)計(jì)分析高、低近交組各繁殖性狀，通過獨(dú)立樣本t(two-sample t-test)檢驗(yàn)方法估計(jì)兩組間各性狀差異顯著性，當(dāng)P<0.05時(shí)差異顯著，當(dāng)P<0.01時(shí)差異極顯著；根據(jù)測序Reads中覆蓋到mC位點(diǎn)的Reads數(shù)與總Reads數(shù)的比值，計(jì)算全基因組甲基化水平；根據(jù)clean reads中唯一比對(duì)的reads，運(yùn)用R語言psych包對(duì)各樣本進(jìn)行主成分(principle components analysis, PCA)分析。

2 結(jié) 果

2.1 狼山雞高、低近交組繁殖性能差異

比較分析狼山雞高、低近交組間各繁殖性狀的差異，發(fā)現(xiàn)高、低近交組在開產(chǎn)日齡、300天產(chǎn)蛋數(shù)和開產(chǎn)體重3個(gè)性狀上均存在顯著或極顯著差異(P<0.05或P<0.01)(表2)，而在開產(chǎn)蛋重上差異不顯著(P>0.05)。

2.2 WGBS數(shù)據(jù)概述以及狼山雞全基因組DNA甲基化水平分析

測序樣本總共12個(gè)，包含4個(gè)組，每組3個(gè)生物學(xué)重復(fù)，每個(gè)樣本測序深度≥30×，產(chǎn)生約154.9 Mb clean reads，平均69.14% 的clean reads可以唯一比對(duì)到參考基因組中，這些clean reads用于后續(xù)的分析。為證明樣本采集的準(zhǔn)確性和可靠性，對(duì)其進(jìn)行主成分分析(PCA)，結(jié)果表明，來自同一組的3個(gè)生物學(xué)重復(fù)相距較近(圖1)，表明試驗(yàn)樣本采集以及測序數(shù)據(jù)較為準(zhǔn)確、可靠。

根據(jù)測序獲得reads中覆蓋到mC位點(diǎn)的reads數(shù)以及reads總數(shù)，計(jì)算并分析了狼山雞高近交組和低近交組中卵巢和下丘腦全基因組DNA甲基化水平，結(jié)果發(fā)現(xiàn)，狼山雞高近交組卵巢和下丘腦全基因組甲基化水平分別為3.49%和4.31%，低近交組卵巢和下丘腦全基因組甲基化水平分別為3.48%和4.35%。相同組織，狼山雞高近交組與低近交組全基因組甲基化水平差異不顯著(P>0.05)，然而兩種組織間比較發(fā)現(xiàn)，下丘腦基因組整體甲基化水平顯著高于卵巢組織(P<0.05，圖2)。

表2 高、低近交組狼山雞各繁殖性狀差異

2.3 狼山雞高、低近交組間差異甲基化區(qū)域檢測

差異甲基化區(qū)域(DMR)是表觀遺傳變化的一個(gè)重要標(biāo)志，涉及到許多生物學(xué)過程調(diào)控。為探索與狼山雞繁殖性能近交衰退相關(guān)的差異甲基化區(qū)域, 本研究篩選高近交組和低近交組間甲基化水平差異大于等于2倍，且 FDR<0.05的區(qū)域作為差異甲基化區(qū)域。結(jié)果在兩組狼山雞的下丘腦和卵巢中分別檢測到5 948和4 593個(gè)差異甲基化區(qū)域，其中分別有1 798和995個(gè)差異甲基化區(qū)域位于CpG 島區(qū)。下丘腦CpG 島區(qū)檢測到的差異甲基化區(qū)域中，有809個(gè)上調(diào)，989個(gè)下調(diào)；卵巢CpG 島區(qū)的995個(gè)差異甲基化區(qū)域中，有374個(gè)上調(diào)，621個(gè)下調(diào)(圖3)。

脊椎動(dòng)物CpG 島區(qū)DNA甲基化在基因轉(zhuǎn)錄、沉默等過程中起重要調(diào)控作用，因此，對(duì)下丘腦和卵巢基因組CpG 島區(qū)差異甲基化區(qū)域進(jìn)行注釋分析發(fā)現(xiàn)，下丘腦和卵巢基因組CpG 島區(qū)差異甲基化區(qū)域分別注釋到了1 020 和 552個(gè)基因，這些基因稱為DMR相關(guān)基因或差異甲基化基因。

2.4 下丘腦基因組CpG島區(qū)差異甲基化基因功能注釋分析

對(duì)下丘腦基因組CpG 島區(qū)差異甲基化基因進(jìn)行GO和KEGG注釋和功能富集分析發(fā)現(xiàn)，分別有920和686個(gè)基因被注釋到了GO和KEGG數(shù)據(jù)庫。GO富集分析結(jié)果顯示，有141個(gè)GO條目被這些基因顯著富集(P<0.05)，這些GO條目主要涉及到信號(hào)轉(zhuǎn)導(dǎo)、神經(jīng)系統(tǒng)發(fā)育、轉(zhuǎn)錄調(diào)控、蛋白激酶活性等生物學(xué)過程，另外，卵母細(xì)胞成熟調(diào)控以及生殖系統(tǒng)發(fā)育相關(guān)的GO條目也被顯著富集到。KEGG富集分析表明，這些差異甲基化基因顯著富集于19條KEGG信號(hào)通路中(P<0.05)，包括轉(zhuǎn)化生長因子β信號(hào)通路、乙型肝炎、脂肪酸代謝、胰島素信號(hào)通路、調(diào)控多能性干細(xì)胞信號(hào)通路、糖酵解/糖異生、氨基酸代謝等信號(hào)通路(圖4)，其中Smad7a、SMAD6、BMP2、TGFB2、MAPK6等繁殖相關(guān)基因包含在這些通路中。

2.5 卵巢基因組CpG島區(qū)差異甲基化基因功能注釋分析

對(duì)卵巢基因組CpG 島區(qū)差異甲基化基因進(jìn)行GO和KEGG注釋和功能富集分析發(fā)現(xiàn)，分別有500和369個(gè)基因被注釋到了GO和KEGG數(shù)據(jù)庫。GO富集分析結(jié)果表明，有105個(gè)GO條目被這些基因顯著富集(P<0.05)，其中涉及到神經(jīng)系統(tǒng)發(fā)育、細(xì)胞增殖、離子轉(zhuǎn)運(yùn)、轉(zhuǎn)錄調(diào)控等(P<0.05)。KEGG富集分析結(jié)果顯示，12條信號(hào)通路被顯著富集(P<0.05)，包括慢性骨髓白血病、流感A、精氨酸和脯氨酸代謝、黏著連接、調(diào)控多能性干細(xì)胞信號(hào)通路等(圖5)。另外，還發(fā)現(xiàn)一些與卵子發(fā)育和性激素分泌相關(guān)的通路被富集到，如黃體酮介導(dǎo)的卵母細(xì)胞成熟、卵母細(xì)胞減數(shù)分裂、GnRH信號(hào)通路、雌激素信號(hào)通路等，其中包含CDC27、ADCY8、AKT3等10個(gè)差異甲基化基因(表3)。

2.6 差異甲基化區(qū)域BSP驗(yàn)證

為驗(yàn)證全基因組重亞硫酸鹽測序分析結(jié)果的可靠性，隨機(jī)挑選2個(gè)CpG島區(qū)差異甲基化區(qū)域(Chr27:2838112-2838584和Chr2:79746219-79746715)，這兩個(gè)差異甲基化區(qū)域分別位于CDC27和SRD5A1基因的CpG島中。采用BSP—克隆測序法對(duì)這兩個(gè)區(qū)域甲基化情況進(jìn)行檢測，結(jié)果顯示，在SRD5A1和CDC27基因的驗(yàn)證區(qū)域分別檢測到了6和14個(gè)甲基化CpG位點(diǎn)(圖6B)，且這些位點(diǎn)的平均甲基化水平在高近交組與低近交組間的變化趨勢與全基因組重亞硫酸鹽測序結(jié)果(圖6A)相一致，表明了全基因組重亞硫酸鹽測序分析結(jié)果的準(zhǔn)確性和可靠性。

表3 卵巢中富集的繁殖相關(guān)通路和CpG 島區(qū)差異甲基化基因

3 討 論

DNA甲基化尤其是脊椎動(dòng)物CpG島區(qū)的甲基化在調(diào)控基因轉(zhuǎn)錄以及轉(zhuǎn)錄后生物學(xué)過程中發(fā)揮重要作用[21]，大量研究表明，DNA甲基化與動(dòng)植物的近交衰退現(xiàn)象密切相關(guān)[10-12]。本研究基于課題組前期狼山雞高、低近交組的成功組建，以及狼山雞高近交組個(gè)體出現(xiàn)的繁殖性能近交衰退現(xiàn)象，通過全基因組重亞硫酸鹽測序技術(shù)開展狼山雞高近交組與低近交組間性腺軸組織DNA甲基化差異比較分析，試圖揭示DNA甲基化在狼山雞繁殖性能近交衰退中的調(diào)控作用。

通過比較分析狼山雞高、低近交組卵巢和下丘腦基因組整體甲基化水平差異，發(fā)現(xiàn)高近交組和低近交組間同一組織基因組整體甲基化水平差異不顯著(P>0.05)，然而，將基因組劃分成許多個(gè)小的區(qū)域，分析各區(qū)域高、低近交組間的甲基化水平差異時(shí)，卻發(fā)現(xiàn)了大量的上調(diào)或下調(diào)的差異甲基化區(qū)域，可能是由于這些區(qū)域上、下調(diào)趨勢正好相反，導(dǎo)致整個(gè)基因組甲基化水平在高、低近交組間差異不顯著。差異甲基化區(qū)域是表觀遺傳變化的一個(gè)重要標(biāo)志，被看作為參與基因轉(zhuǎn)錄水平調(diào)控的功能性區(qū)域，尤其是脊椎動(dòng)物CpG島區(qū)DNA甲基化的變化一直以來受到研究者們的極大關(guān)注[33-36]。本研究在下丘腦和卵巢基因組CpG島區(qū)分別檢測到了1 798和995個(gè)差異甲基化區(qū)域，這些差異甲基化區(qū)域分別注釋到了1 020 和 552個(gè)基因，為了進(jìn)一步探索這些差異甲基化區(qū)域的功能，對(duì)差異甲基化區(qū)域相關(guān)基因進(jìn)行了功能注釋和富集分析，發(fā)現(xiàn)了大量差異甲基化基因富集在神經(jīng)系統(tǒng)發(fā)育、信號(hào)轉(zhuǎn)導(dǎo)和離子轉(zhuǎn)運(yùn)等神經(jīng)內(nèi)分泌調(diào)節(jié)相關(guān)生物學(xué)過程GO條目，這些生物學(xué)過程與雞的繁殖性能密切相關(guān)。KEGG功能富集分析顯示，卵巢中一些差異甲基化基因富集在了卵子發(fā)育和性激素分泌相關(guān)的通路，如黃體酮介導(dǎo)的卵母細(xì)胞成熟、卵母細(xì)胞減數(shù)分裂、GnRH信號(hào)通路、雌激素信號(hào)通路等，包含CDC27、ADCY8、AKT3等10個(gè)差異甲基化基因，推測這些基因CpG島區(qū)的差異甲基化在狼山雞繁殖性能近交衰退中可能發(fā)揮重要調(diào)控作用[37]。其中，CDC27基因CpG島區(qū)甲基化差異情況進(jìn)一步得到了BSP-克隆測序方法的驗(yàn)證，CDC27是編碼細(xì)胞周期進(jìn)程所必需的一種蛋白質(zhì)[38-39]，被認(rèn)為與細(xì)胞增殖和細(xì)胞分裂相關(guān)[40]，本研究中，該基因CpG島區(qū)甲基化水平在狼山雞高近交組中顯著上調(diào)，推測該基因CpG島區(qū)甲基化上調(diào)影響了其轉(zhuǎn)錄表達(dá)，進(jìn)而影響了卵母細(xì)胞分裂和成熟，導(dǎo)致狼山雞高近交組繁殖性能衰退。

另外，在下丘腦和卵巢差異甲基化基因顯著富集的通路中，發(fā)現(xiàn)多條通路與氨基酸、脂肪和糖代謝相關(guān)，推測狼山雞的高度近交可能影響了其基礎(chǔ)代謝，進(jìn)而可能間接影響到其生殖系統(tǒng)發(fā)育，導(dǎo)致其繁殖性能下降。一些與疾病相關(guān)的信號(hào)通路也被下丘腦和卵巢差異甲基化基因顯著富集到，如乙型肝炎、慢性骨髓白血病、流感A等，表明狼山雞高度近交后，其抗病力或免疫力方面也受到了影響，可能需要后續(xù)進(jìn)一步關(guān)注和研究。

4 結(jié) 論

狼山雞高、低近交組間基因組整體DNA甲基化水平差異不顯著(P>0.05)，但全基因組范圍內(nèi)仍檢測到大量差異甲基化區(qū)域，且發(fā)現(xiàn)大量差異甲基化基因富集到繁殖相關(guān)的生物學(xué)過程和通路中，如CDC27、ADCY8、AKT3等差異甲基化基因，推測這些基因CpG島區(qū)甲基化在狼山雞繁殖性能近交衰退過程中發(fā)揮重要作用。本研究結(jié)果為今后深入研究雞繁殖性能近交衰退調(diào)控機(jī)制提供了有力線索，為物種資源保護(hù)工作提供了科學(xué)依據(jù)。