不同剩余采食量水平的安格斯牛下丘腦單核苷酸多態(tài)性和可變剪接分析

周小南， 馬應(yīng)， 楊朝云， 丁燕玲， 張巖峰， 趙志艷， 康曉龍

(寧夏大學(xué)農(nóng)學(xué)院，寧夏 銀川 750021)

剩余采食量(residual feed intake, RFI)是指實際采食量和預(yù)期的飼料需要量的差值。RFI水平與飼料轉(zhuǎn)化率和干物質(zhì)采食量顯著相關(guān)[1]。當(dāng)RFI的數(shù)值為負(fù)數(shù)時，表明畜禽實際采食量低于預(yù)期采食量，該畜禽為高飼料利用率動物，反之則為低飼料利用率動物。RFI水平不僅反映了個體由遺傳背景所決定的代謝差異，也校正了試驗個體的代謝體質(zhì)量。楊朝云等[2]對2個RFI水平(高/低)的安格斯牛與生長性狀的相關(guān)性分析發(fā)現(xiàn)，低RFI水平組的安格斯牛比高RFI水平組節(jié)約8.8%的飼料。CONNOR等[3]研究發(fā)現(xiàn)，當(dāng)?shù)蚏FI水平組肉牛的飼料攝入量減少10%～12%，此時溫室氣體排放量減少25%～30%，糞便中的養(yǎng)分損失減少15%～17%。單核苷酸多態(tài)性(single nucleotide polymorphism site, SNPs)是指單核苷酸堿基變異，由轉(zhuǎn)換(C/T或G/A)或顛換(C/G、C/A或T/A、T/G)引起[4]。近年來，許多研究利用表達(dá)數(shù)據(jù)進(jìn)行表型相關(guān)表達(dá)基因的SNP分析[5]，鑒定差異表達(dá)基因單核苷酸變異的潛在作用，為進(jìn)一步整合表達(dá)數(shù)據(jù)與基因組信息以及分析其與表型相關(guān)性提供了參考。SAATCHI等[6]通過對4個獨立肉牛群體進(jìn)行全基因組關(guān)聯(lián)研究和功能分析，鑒定出影響肉牛生長的主效數(shù)量性狀基因座(quantitative trait loci，QTL)。在植物方面，王賽等[7]基于SNP遺傳圖譜定位玉米雄穗分枝數(shù)和主軸長QTLs。研究表明，基因的不同剪接方式也會影響機體發(fā)育和生物功能，可變剪接(alternative splicing, AS)是Pre-mRNA通過不同的剪接方式形成序列不同的mRNA剪接異構(gòu)體的過程，其中單個基因通過AS產(chǎn)生多個mRNA異構(gòu)體和蛋白質(zhì)亞型，從而增加蛋白質(zhì)多樣性和復(fù)雜性[8]。周揚[9]對不同發(fā)育階段和不同性別秦川牛的皮下脂肪組織內(nèi)AS相關(guān)基因表達(dá)特征進(jìn)行分析發(fā)現(xiàn)，有28%的基因在牛脂肪組織中出現(xiàn)AS事件。在植物中，AS相關(guān)基因參與調(diào)控植物的抗性特征[10]。下丘腦作為動物采食量調(diào)節(jié)中心，存在復(fù)雜的調(diào)節(jié)食欲和能量平衡的神經(jīng)元網(wǎng)絡(luò)[11]。因此，本研究對安格斯牛下丘腦組織進(jìn)行高通量測序，分析不同RFI水平安格斯牛下丘腦組織的SNPs和表達(dá)基因?qū)?yīng)的AS事件，梳理下丘腦中由于RFI水平不同而發(fā)生的AS事件及AS相關(guān)基因的功能和參與的通路，以期挖掘AS在肉牛剩余采食量表型調(diào)控中的潛在生物學(xué)作用，為后續(xù)解析肉牛飼料利用效率相關(guān)的遺傳學(xué)基礎(chǔ)提供理論支撐。

1 材料與方法

1.1 試驗動物

試驗動物選自寧夏地區(qū)某養(yǎng)殖場的安格斯牛，初始體質(zhì)量為(266.7±35.7)kg，所有安格斯牛(n=30)均按照標(biāo)準(zhǔn)程序飼養(yǎng)，飼喂時期為81 d。飼養(yǎng)結(jié)束后，依據(jù)標(biāo)準(zhǔn)方法計算RFI水平，選擇RFI水平最高(high residual feed intake, HRFI)和水平最低(low residual feed intake, LRFI)個體各5頭，分別編號為H1、H2、H3、H4、H5和L1、L2、L3、L4、L5，試驗牛只頸部放血屠宰開顱后，立即從下丘腦組織收集1 cm3的組織，并用生理鹽水沖洗干凈，剪碎后裝入無酶的凍存管，于-80 ℃保存?zhèn)溆?。試驗遵循寧夏大學(xué)動物保護(hù)委員會的指導(dǎo)原則，所有組織均按照無害化處理方法和原則進(jìn)行科學(xué)處理。

1.2 RNA的提取和測序

取適量下丘腦組織樣品，采用TRIzol法分別提取HRFI組和LRFI組樣本的總RNA，將提取到的RNA以質(zhì)量分?jǐn)?shù)1%的瓊脂糖凝膠進(jìn)行檢測，利用 Nano Photometer?分光光度計測定OD260/280及OD260/230比值，將OD260/280介于1.9和2.0之間的RNA樣品用于RNA測序。構(gòu)建cDNA文庫，將構(gòu)建合格的不同文庫按照有效濃度及目標(biāo)下機數(shù)據(jù)量的需求富集后進(jìn)行Illumina測序，RNA測序由北京諾禾致源科技股份有限公司完成。

1.3 測序數(shù)據(jù)分析

對測序得到的原始數(shù)據(jù)(raw reads)進(jìn)行質(zhì)控，過濾去除reads中的接頭序列、帶N堿基以及低質(zhì)量reads序列，得到clean reads。

1.4 SNP及可變剪接分析

將比對結(jié)果進(jìn)行染色體坐標(biāo)排序、去掉重復(fù)的reads等處理后，通過變異檢測軟件GATK2分析各個樣品的SNP位點，同時對原始結(jié)果進(jìn)行過濾。使用rMATS軟件對差異表達(dá)的基因進(jìn)行AS分析，以錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)<0.05作為差異表達(dá)AS事件的篩選標(biāo)準(zhǔn)。經(jīng)rMATS軟件分析包括5種可變剪接類型(圖1)，分別為外顯子跳躍(skipped exon，SE)、5’端可變剪接(alternative 5’ splice site，A5SS)、3’端可變剪接(alternative 3’ splice site，A3SS)、外顯子互斥(mutually exclusive exon，MXE)和內(nèi)含子保留(retained intron，RI)。

注：A，外顯子跳躍；B，外顯子互斥；C，5′端可變剪接；D，3′端可變剪接；E，內(nèi)含子保留。

1.5 差異表達(dá)AS相關(guān)基因功能富集分析

為確定差異表達(dá)AS相關(guān)基因的功能，采用在線分析軟件g:profiler (https://biit.cs.ut.ee/gprofiler/)進(jìn)行基因本體(gene ontology, GO)通路富集分析，同時利用DAVID (https://david.ncifcrf.gov/)功能注釋進(jìn)行京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes，KEGG)通路富集分析，顯著性水平設(shè)置為P<0.05。

1.6 實時熒光定量PCR驗證試驗結(jié)果

為驗證轉(zhuǎn)錄測序的可靠性，隨機挑選7個差異表達(dá)AS相關(guān)基因(PARP2、ITGB1、BPIFBI、LRP11、USP34、MGAT1和KCNN2)進(jìn)行實時熒光定量PCR(quantitative real-time PCR，qRT-PCR)，以GAPDH基因作為內(nèi)參基因，引物委托上海生工生物工程有限公司合成(表1)。配置20 μL的反應(yīng)體系: TB Green Taq II 10 μL，Rox Dye II 0.4 μL，cDNA 2 μL，上下游引物各0.8 μL，dd H2O 6 μL。反應(yīng)程序：95 ℃ 預(yù)變性 4 min；95 ℃變性 10 s，最佳退火溫度30 s，72 ℃延伸10 s，共39個循環(huán)；65 ℃ 5 s，95 ℃ 5 s。使用 2-△△Ct法對基因的表達(dá)量進(jìn)行標(biāo)準(zhǔn)化計算。差異性分析使用SPSS 2.2軟件，以P<0.05 為差異顯著判斷標(biāo)準(zhǔn)。

表1 熒光定量PCR驗證引物序列Table 1 Primer sequences verified by fluorescence quantitative PCR

2 結(jié)果與分析

2.1 高/低RFI水平安格斯牛下丘腦組織轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量分析

對HRFI組和LRFI組的安格斯牛下丘腦組織RNA樣品進(jìn)行測序，結(jié)果如表2所示。所得的raw reads在145 924 188～159 000 376之間，所有樣本平均clean reads比例高于88.48%，Q30堿基比例介于89.48%～93%之間，表明文庫構(gòu)建成功且轉(zhuǎn)錄組測序得到的結(jié)果質(zhì)量良好，測序數(shù)據(jù)可靠，可用于后續(xù)分析。

表2 高/低RFI水平安格斯牛下丘腦組織轉(zhuǎn)錄組測序數(shù)據(jù)評估統(tǒng)計Table 2 Statistics of sequencing data evaluation of Angus cattle with high/low RFI levels

2.2 高/低RFI水平安格斯牛SNP位點統(tǒng)計分析

使用變異檢測軟件GATK2檢測各個樣品潛在的SNP位點，同時進(jìn)行SNP calling和InDel calling分析，結(jié)果如表3所示。根據(jù)測到的SNP位點，統(tǒng)計每組的突變類型的頻率，對2個RFI水平(高/低)的安格斯牛10個樣本進(jìn)行SNP突變類型統(tǒng)計，結(jié)果如圖2—圖4所示。10個樣本中，C/T和A/G的突變頻率最高。所有樣品中檢測到的SNP數(shù)量為48 227～1 000 092個。分析2個RFI水平(高/低)的樣本，發(fā)生頻率最高的是轉(zhuǎn)換型SNP，比例約為71%，而顛換型SNP位點比例約為29%，轉(zhuǎn)換和顛換的比例為2.45，所有樣品中發(fā)生InDel突變事件的數(shù)目為48 276～87 677個。

表3 高/低RFI水平安格斯牛SNP位點分析結(jié)果Table 3 Analysis results of SNP locus of Angus cattle with high/low RFI levels

圖2 高RFI水平安格斯牛SNP突變頻率分布

圖3 低RFI水平安格斯牛SNP突變頻率分布

圖4 高/低RFI水平安格斯牛轉(zhuǎn)換和顛換SNP類型頻率比例

2.3 高/低RFI水平安格斯牛AS差異分析

采用rMATS軟件對HFRI組和LRFI組進(jìn)行AS差異分析，共有12 385個基因?qū)?yīng)35 463個AS事件，鑒定出差異表達(dá)AS事件498個(表4)。2個RFI水平(高/低)發(fā)生的AS事件中，SE類型最普遍，共30 262個，占85.33%；其次是MXE類型共5 057個，占14.26%；A5SS、A3SS和RI類型較少，共占0.41%，分別為33、78和33個。進(jìn)一步對AS事件進(jìn)行顯著性分析發(fā)現(xiàn)，2個RFI水平(高/低)之間差異表達(dá)SE類型有416個，包括199個上調(diào)和217個下調(diào)。此外，差異表達(dá)MXE類型有79個，包括43個上調(diào)和36個下調(diào)；差異表達(dá)A3SS類型有1個下調(diào)；差異表達(dá)RI事件有2個，包括1個上調(diào)和1個下調(diào)；沒有差異表達(dá)A5SS類型。對發(fā)生AS事件的外顯子長度進(jìn)行分析發(fā)現(xiàn)(表5)，RI類型的外顯子長度最長，平均數(shù)為887.9 bp，中位數(shù)為512.0 bp；另外4種類型的AS事件外顯子長度較短，平均數(shù)在136.9～186.0 bp之間。

表4 高/低RFI水平安格斯牛AS差異分析統(tǒng)計Table 4 Differential AS statistics of Angus cattle with high/low RFI levels

表5 高/低RFI水平安格斯牛AS外顯子長度Table 5 AS exon length of Angus cattle with high/low RFI levels

2.4 高/低RFI水平安格斯牛差異表達(dá)AS相關(guān)基因GO和KEGG通路富集分析

對獲得的高/低RFI水平安格斯牛差異表達(dá)AS相關(guān)基因進(jìn)行GO通路富集分析(圖5)，共檢索到34條GO條目，包括分子功能6條，細(xì)胞組成11條，生物學(xué)過程17條。在分子功能中有35個差異表達(dá)AS相關(guān)基因富集在ATP結(jié)合方面，生物學(xué)過程中分別有11個和8個差異表達(dá)AS相關(guān)基因富集在光感受器連接體和突觸纖毛組裝。對顯著差異表達(dá)的AS相關(guān)基因進(jìn)行KEGG通路富集分析發(fā)現(xiàn)(圖6)，共有82個基因顯著富集在13條通路上，包括胰島素信號通路、催乳激素信號通路、膽堿代謝、哺乳動物雷帕霉素靶蛋白(mammalian target of rapamycin，mTOR)信號通路、環(huán)磷酸腺苷(cyclic adenosine monophosphate，cAMP)信號通路、酪氨酸激酶受體(erbB)和晝夜節(jié)律等。

圖6 高/低RFI水平安格斯牛差異表達(dá)基因KEGG通路富集散點圖Fig.6 Enrichment distribution point diagram of differentially expressed gene KEGG pathway of Angus cattle with high/low RFI levels

2.5 差異表達(dá)AS相關(guān)基因測序結(jié)果驗證

通過qRT-PCR驗證PARP2、ITGB1、BPIFBI、LPR11、USP34、MGAT1和KCNN2共7個差異表達(dá)AS相關(guān)基因的表達(dá)情況(圖7)。經(jīng)qRT-PCR得到的基因相對表達(dá)量和RNA-seq測序結(jié)果的差異倍數(shù)基本一致，通過計算2個結(jié)果的皮爾遜相關(guān)性，相關(guān)系數(shù)r=0.95，表明RNA-seq數(shù)據(jù)相對可靠。

圖7 qRT-PCR 對RNA-Seq檢測出的差異表達(dá)基因驗證

3 結(jié)論與討論

下丘腦內(nèi)存在復(fù)雜的調(diào)節(jié)食欲和能量平衡的神經(jīng)元網(wǎng)絡(luò)[12]，通過促進(jìn)或抑制攝食神經(jīng)肽的合成作用調(diào)節(jié)采食量[13]。本研究通過對安格斯牛下丘腦的轉(zhuǎn)錄組測序分析，分析不同RFI水平的AS相關(guān)基因。這些基因均與食欲調(diào)節(jié)相關(guān)，如阿黑皮素原(proopiomelanocortin, POMC)相關(guān)基因和生長激素釋放激素受體(growth hormone-releasing hormone receptor , GHRHR)相關(guān)基因均發(fā)生了AS事件[14-15]。本研究通過對2個RFI水平(高/低)安格斯牛下丘腦組織樣品測序數(shù)據(jù)進(jìn)行SNP檢測分析發(fā)現(xiàn)，A/T和C/G的突變頻率最高，約占總數(shù)的71%;而顛換型SNP位點比例為29%，轉(zhuǎn)換和顛換的比例為2.45。這與其他動物如鴨、雞和牛的研究結(jié)果一致[16-18]。SNP是基因組中最常見的遺傳變異之一，但是轉(zhuǎn)換偏差現(xiàn)象使得每個SNP 位點發(fā)生變異形式由4種變?yōu)橹挥?種，即轉(zhuǎn)換型和顛換型，且比例為1∶2[19-21]。該現(xiàn)象可能是物種長期進(jìn)化的結(jié)果，且這種變異在一定程度上有利于保持編碼蛋白原有的結(jié)構(gòu)，減少不利突變的概率[22]。因此，在SNP數(shù)據(jù)集中觀察到的轉(zhuǎn)換和顛換的比例相對較高。此外，本研究結(jié)果顯示，共有12 385個基因發(fā)生35 463個AS事件，這也說明了單個基因通過AS可以形成不同的mRNA剪接異構(gòu)體。對AS事件的類型進(jìn)行統(tǒng)計，其中SE類型比例最高，達(dá)85.33%；其次是MXE類型，比例為14.30%。王旭平等[23]對番鴨下丘腦轉(zhuǎn)錄組數(shù)據(jù)分析發(fā)現(xiàn)，有5 034個基因發(fā)生7 528個AS事件，其中SE類型占比92.76%，而A5SS類型和RI類型僅占0.027%。鮑晶晶等[24]對綿陽不同發(fā)育時期的背最長肌組織的AS事件進(jìn)行統(tǒng)計，其中SE類型比例最高，為70.69%，這與本試驗的研究結(jié)果一致。本研究統(tǒng)計了各類型AS事件可變外顯子的長度，其中發(fā)生頻率較高的SE和MXE類型的AS事件的外顯子長度平均數(shù)較小，而A5SS和RI類型的可變外顯子平均長度較大，這表明較小的外顯子更容易發(fā)生AS事件，該結(jié)果與STAMM等[25]的研究結(jié)果一致。

在本研究中，從2個RFI水平(高/低)安格斯牛中檢測出35 463個AS事件，將差異表達(dá)的AS相關(guān)基因比對到KEGG數(shù)據(jù)庫中，共篩選到13條AS相關(guān)基因顯著富集的通路。其中，顯著富集基因最多的是胰島素信號通路，其次是催乳激素信號通路，還有較多的基因富集在膽堿代謝、mTOR信號通路、cAMP信號通路、erbB信號通路和晝夜節(jié)律等通路，這些通路廣泛參與調(diào)節(jié)能量代謝[26-28]。其中，mTOR信號通路在細(xì)胞的增殖、生長和分化過程起到重要的調(diào)控作用[29-30]。研究表明，mTOR可監(jiān)測下丘腦中營養(yǎng)和激素傳輸信號，當(dāng)外界刺激增強下丘腦mTOR信號通路活性時能夠抑制動物的采食，但當(dāng)在腦室注射mTOR信號通路的抑制劑則顯著增加禁食動物的短期采食量[31]。同時，本研究結(jié)果顯示，有8個差異表達(dá)AS相關(guān)基因顯著富集于cAMP信號通路。cAMP作為第二信使在脂肪分解中具有重要作用，其表達(dá)受Ca2+水平的影響。在豬脂肪組織中，cAMP磷酸二酯酶作為中間體，影響脂解過程[28]。此外，喂養(yǎng)治療試驗證實了肥胖和cAMP途徑之間的可能聯(lián)系[32]。本研究結(jié)果揭示AS事件及表達(dá)基因的SNPs在安格斯牛采食量調(diào)控中的潛在功能，表明不同剩余采食量水平的差異表達(dá)AS相關(guān)基因可能通過能量代謝相關(guān)的信號通路和影響食欲的通路來調(diào)節(jié)能量代謝，進(jìn)而參與肉牛采食量性狀，但具體生物學(xué)機制有待進(jìn)一步的深入研究與探討。