基于RNA-Seq技術(shù)挖掘貴州白山羊公羊背最長肌生長發(fā)育下調(diào)基因

安清明，王 星，吳震洋，孟金柱，趙園園 (．銅仁學(xué)院 農(nóng)林工程與規(guī)劃學(xué)院 貴州省梵凈山地區(qū)生物多樣性保護(hù)與利用重點實驗室，貴州 銅仁 554300；2．甘肅農(nóng)業(yè)大學(xué) 草食動物生物技術(shù)重點實驗室，甘肅 蘭州 730070)

目前雖然已有研究對家畜動物的肌肉生長進(jìn)行了篩選，獲得了部分相關(guān)的候選基因，但關(guān)于地方山羊品種肌肉生長發(fā)育的研究較少，尤其對不同性別的山羊肌肉生長發(fā)育還有待深入研究。因此通過研究不同性別貴州白山羊背最長肌組織中差異表達(dá)基因，挖掘和鑒定控制白山羊肌肉組織生長發(fā)育的相關(guān)基因，對改良貴州白山羊胴體肌肉生長發(fā)育具有重要意義。本研究采用RNA-Seq技術(shù)對不同性別貴州白山羊背最長肌進(jìn)行分析，篩選與公羊肌肉生長相關(guān)的下調(diào)基因，并對篩選的下調(diào)基因進(jìn)行GO、KEGG等分析，并通過RT-PCR進(jìn)行驗證，挖掘出可調(diào)節(jié)貴州白山羊肌肉生長發(fā)育的關(guān)鍵候選基因，以期為貴州白山羊在肉質(zhì)和肌肉生長性狀方面的分子機(jī)理研究提供理論依據(jù)。

1 材料與方法

1.1 樣品采集樣品采自貴州省銅仁市沿河土家族自治縣華珍牧業(yè)有限責(zé)任公司，選取同等飼養(yǎng)條件下2周歲齡、健康的白山羊6只(公、母羊各3只)，屠宰后采集背最長肌置于液氮中迅速冷凍，帶回實驗室后置于-80℃保存?zhèn)溆谩?/p>

1.2 背最長肌RNA提取、文庫構(gòu)建及轉(zhuǎn)錄組測序按照TRIzol法提取每只白山羊總RNA，用RNAprep pure Tissue Kit(TianGen，天根生化科技(北京)有限公司)純化，采用Agilent Bioanalyzer 2100檢測完整性，Qubit 2.0 Flurometer測量濃度，待樣品合格后送北京諾和致源生物信息科技有限公司構(gòu)建文庫，應(yīng)用Illumina Hiseq 2500平臺進(jìn)行轉(zhuǎn)錄組測序。

1.3 數(shù)據(jù)處理與分析將獲得的原始Reads應(yīng)用FastQC軟件(http://www.bioinformatics.babraha- m.ac.uk/pr-ojects/fastqc/)進(jìn)行質(zhì)量評估，去除數(shù)據(jù)中帶有接頭、含N和低質(zhì)量的原始片段，保證獲得的數(shù)據(jù)具有品質(zhì)較高的Clean Reads。應(yīng)用Trinity(http://trinityrnaseq.sourceforge.net/)對Clean Reads進(jìn)行組裝，獲得unigenes。利用HISAT2軟件將獲得數(shù)據(jù)與已知山羊參考基因組進(jìn)行比對，獲得序列在參考基因組上的相關(guān)信息。應(yīng)用DESeq2軟件對獲得mRNA進(jìn)行統(tǒng)計分析，篩選差異表達(dá)的下調(diào)基因mRNA。利用GO seq軟件對篩選的下調(diào)表達(dá)基因進(jìn)行GO分析，應(yīng)用KOBAS軟件進(jìn)行KEGG通路分析。

1.4 cDNA及引物合成將轉(zhuǎn)錄組測序剩余總RNA反轉(zhuǎn)錄為cDNA：按照北京全式金生物技術(shù)有限公司反轉(zhuǎn)錄試劑盒說明進(jìn)行，42℃孵化15 min，再85℃加熱5 s。從轉(zhuǎn)錄組測序篩選出的差異下調(diào)表達(dá)基因中隨機(jī)挑選4個基因進(jìn)行驗證(KLF11、TNS1、SIK2和IGF1R)，以GAPDH基因為參考基因，引物由生工生物工程(上海)股份有限公司合成，引物序列見表1。

表1 熒光定量引物信息

1.5 熒光定量PCR驗證通過熒光定量驗證轉(zhuǎn)錄組測序中貴州白山羊公羊(MF)與母羊(FL)背最長肌中差異表達(dá)下調(diào)mRNA的相對表達(dá)量。每個樣品設(shè)3個重復(fù)，應(yīng)用北京全式金生物技術(shù)有限公司試劑盒進(jìn)行相對定量檢測分析，反應(yīng)體系：cDNA模板1 μL，2×qPCR super mix 10 μL，上、下游引物各0.5 μL，RNA-free H2O補(bǔ)齊至20 μL。反應(yīng)條件：94℃ 1 min，94℃ 10 s，59℃ 30 s，72℃ 10 s，43個循環(huán)。使用2-△△Ct法計算各基因間的相對表達(dá)水平。

2 結(jié)果

2.1 RNA-seq數(shù)據(jù)比對分析對采集6個背最長肌樣品的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行分析，對原始獲得的Raw data將帶接頭、含N和低質(zhì)量的Reads片段過濾篩選后，最終在MF中獲得95739024條Clean reads，結(jié)果中堿基質(zhì)量在Q30以上的占92.86%；在FL中獲得89 806 996條Clean reads，結(jié)果中堿基質(zhì)量在Q30以上的占89.86%(圖1)。將獲得的Clean reads與山羊參考基因組序列進(jìn)行比對，結(jié)果表明MF中比對效率為94.18%，F(xiàn)L比對效率為94.41%，其中在參考序列中能夠唯一比對到位置的Reads數(shù)量分別占86.24%和86.07%，能夠多點定位的Reads分別占7.94%和8.34%，表明測序質(zhì)量高，數(shù)據(jù)可靠，可用于后續(xù)分析(表2)。

表2 測序數(shù)據(jù)與參考山羊基因組的序列比對結(jié)果

圖1 MF和FL背最長肌中原始數(shù)據(jù)的分類

2.2 差異表達(dá)基因的篩選、功能注釋和富集分析

2.2.1下調(diào)差異表達(dá)基因的篩選 以FPKM值作為衡量轉(zhuǎn)錄基因下調(diào)差異表達(dá)水平的指標(biāo)，分析6個樣本共得到33 896個轉(zhuǎn)錄本，對MF和FL的FPKM進(jìn)行標(biāo)準(zhǔn)化后，篩選出FPKM≥1的轉(zhuǎn)錄本有25 089個。應(yīng)用DESeq2軟件對獲得的25 089個mRNA進(jìn)行差異表達(dá)篩選，設(shè)定參數(shù)FPKM≥1，且P<0.05時，共獲得514差異表達(dá)下調(diào)基因，其中100個為新發(fā)現(xiàn)的基因(轉(zhuǎn)錄本)，表3列出MF背最長肌中FPKM數(shù)值最高的前20個基因，表4列出MF相較于FL背最長肌中差異表達(dá)最高的10個下調(diào)基因及其功能。

表3 MF背最長肌中FPKM數(shù)值最高的20個基因

表4 可能與貴州白山羊背最長肌發(fā)育相關(guān)的候選基因

2.2.2差異表達(dá)基因的GO功能富集 應(yīng)用Goseq軟件對篩選出的514個差異表達(dá)下調(diào)基因進(jìn)行GO功能富集分析，結(jié)果共分為2大類21個功能亞類，主要參與生物學(xué)過程(biological process，BP)和細(xì)胞組分(cellular component，CC)，而分子功能(molecular function，MF)的功能富集在本研究中未檢測到，結(jié)果見表5。

表5 貴州白山羊MF和FL中上調(diào)表達(dá)基因GO分析

2.2.3下調(diào)差異表達(dá)基因KEGG信號通路分析 為了篩選出與貴州白山羊背最長肌生長發(fā)育相關(guān)的信號通路，應(yīng)用Kobas軟件對下調(diào)差異表達(dá)基因進(jìn)行KEGG信號通路分析，F(xiàn)DR≤0.05的通路為差異表達(dá)基因中顯著富集的信號通路，結(jié)果表明被注釋的差異基因共參與145個信號通路，其中參與核糖體信號通路的基因最為富集，共有27個被注釋基因，圖2為下調(diào)表達(dá)基因中KEGG信號通路結(jié)果中顯著性最高的20個通路。

圖2 貴州白山羊MF和FL背最長肌中下調(diào)表達(dá)基因KEGG通路富集散點圖

富集在核糖體信號通路中的差異蛋白位點發(fā)現(xiàn)共有25個(包括27個注釋基因)，核糖體信號通路圖見圖3。25種差異蛋白分別分布在核糖體的60S大亞基和40S小亞基中，其中60S大亞基核糖體蛋白15種(包括17個基因)，40S小亞基核糖體蛋白9種(包括10個基因)，其中S38e小亞基核糖體蛋白未注釋到對應(yīng)基因(表6)。

圖3 核糖體信號通路

表6 核糖體信號通路中的差異蛋白位點及基因

2.2.4下調(diào)差異表達(dá)基因RT-PCR驗證分析 通過功能分析，在差異表達(dá)的下調(diào)基因中，隨機(jī)選擇4個可能與貴州白山羊公羊背最長肌生長發(fā)育呈下調(diào)趨勢的調(diào)控基因KLF11、TNS1、SIK2和IGF1R進(jìn)行驗證，具體信息見表7。采用RT-PCR進(jìn)行驗證這4個基因的相對表達(dá)量(圖4)。

**表示P<0.05圖4 候選基因在不同性別白山羊MF和FL背最長肌中的相對表達(dá)量

表7 貴州白山羊背最長肌負(fù)調(diào)控候選基因

3 討論

肌肉的生長發(fā)育受到多種基因和信號通路的調(diào)節(jié)，同時還受品種、性別、飼養(yǎng)環(huán)境、營養(yǎng)等諸多因素的影響，同時肌肉的生長發(fā)育直接影響動物的體質(zhì)量和肉品質(zhì)。因此，本研究以同等飼養(yǎng)條件下、同年齡段的貴州白山羊公、母羊背最長肌為研究材料，基于RNA-Seq測序技術(shù)篩選與肌肉生長發(fā)育相關(guān)的基因，結(jié)合相關(guān)生物信息學(xué)功能分析，篩選出公羊背最長肌對于母羊背最長肌存在差異的514個下調(diào)基因，其中100個為新發(fā)現(xiàn)的關(guān)聯(lián)基因，這些差異表達(dá)基因可能對不同性別貴州白山羊骨骼肌的生長發(fā)育具有一定影響，值得進(jìn)一步深入研究。

本研究通過功能富集分析發(fā)現(xiàn)，差異最為富集的蛋白均屬于核糖體信號通路。核糖體是機(jī)體細(xì)胞中重要的細(xì)胞器，是合成蛋白質(zhì)的重要場所，核糖體中最主要的成分為核糖體蛋白，其廣泛分布于各種組織中，在核糖體的翻譯和蛋白質(zhì)的合成中發(fā)揮著重要作用。研究還發(fā)現(xiàn)，核糖體蛋白同時還具有調(diào)控基因轉(zhuǎn)錄、調(diào)控mRNA的翻譯、同時還能影響細(xì)胞分化、增殖和凋亡等作用[11]。本研究在核糖體通路中共發(fā)現(xiàn)了25種差異表達(dá)的蛋白，分布于核糖體的兩種不同亞基中，這些差異蛋白是否會影響背最長肌中核糖體的蛋白質(zhì)生物合成過程，進(jìn)而影響肌肉的生長發(fā)育尚未可知，有待進(jìn)一步研究。

在篩選出的514個差異下調(diào)基因中，隨機(jī)選擇4個可能與貴州白山羊公羊背最長肌肌肉生長發(fā)育呈下調(diào)趨勢的調(diào)控基因KLF11、TNS1、SIK2和IGF1R進(jìn)行RT-PCR驗證，結(jié)果與RNA-Seq分析結(jié)果一致，表明RNA-Seq分析結(jié)果可靠，因此篩選出的514個關(guān)聯(lián)基因可用于進(jìn)一步分析。本研究隨機(jī)選擇的4個下調(diào)基因中，KLF11基因編碼的轉(zhuǎn)錄因子KLF11屬于Kruppel轉(zhuǎn)錄因子家族(Kruppel like factorsfamily，KLF)成員，該家族轉(zhuǎn)錄因子成員眾多，目前共發(fā)現(xiàn)18個成員，且功能復(fù)雜多樣，在不同組織和不同生理條件下參與細(xì)胞分化、凋亡、衰老、細(xì)胞間黏附、糖脂代謝和脂肪細(xì)胞分化等生命活動[12-15]。KLF11基因有4個外顯子和3個內(nèi)含子組成，共編碼512個氨基酸，編碼蛋白相對分子質(zhì)量為74 kDa。已有研究表明，KLF11基因編碼的轉(zhuǎn)錄因子可以與不同作用因子結(jié)合，發(fā)揮不同的轉(zhuǎn)錄調(diào)控作用。有研究表明，KLF11可以調(diào)節(jié)TGH-β信號通路的下游信號，從而參與由TGH-β信號通路調(diào)控的生理功能，如細(xì)胞分化凋亡、抑制細(xì)胞增殖等[16-17]。LOMBERK等[18]研究發(fā)現(xiàn)KLF11還能參與調(diào)控葡萄糖攝取、糖酵解和脂質(zhì)代謝。王江林等[19]研究發(fā)現(xiàn)，通過干擾KLF11基因可促進(jìn)山羊肌內(nèi)前體脂肪細(xì)胞分化。綜上，KLF11基因可能是與白山羊肌肉生長發(fā)育的一個重要候選基因。

TNS1基因(tensin 1)編碼蛋白稱張力蛋白1，屬于黏著蛋白家族，該家族蛋白是一種多結(jié)構(gòu)域蛋白，其主要參與細(xì)胞膜交互、構(gòu)建肌動蛋白細(xì)胞骨架，因此又被稱作支架蛋白，同時還能參與細(xì)胞黏著、運動和生長等生物學(xué)活動[20]。已有的研究報道主要集中在TNS1基因與癌癥的關(guān)聯(lián)分析中，如ZHAN等[21]研究發(fā)現(xiàn)TNS1基因高表達(dá)狀態(tài)與乳腺癌細(xì)胞的轉(zhuǎn)移呈負(fù)相關(guān)。曹一鑫等[22]研究發(fā)現(xiàn)TNS1基因高表達(dá)狀態(tài)與胃癌腫瘤細(xì)胞轉(zhuǎn)移呈負(fù)相關(guān)。任琎等[23]研究發(fā)現(xiàn)TNS1基因可以抑制非小細(xì)胞肺癌細(xì)胞增殖和侵襲。在家畜動物中，僅靳聰飛等[24]在牛肌肉組織中克隆獲得了TNS1基因全部序列，通過生物信息學(xué)分析預(yù)測TNS1基因編碼蛋白不屬于分泌蛋白，磷酸化位點均在絲氨酸、蘇氨酸和酪氨酸殘基上。目前尚未見TNS1基因影響肌肉生長發(fā)育的相關(guān)報道，因此，TNS1基因是否可以作為影響動物肌肉組織的生長發(fā)育的候選基因，仍需進(jìn)一步深入研究。

SIK2基因(salt-inducible kinase 2)編碼蛋白稱為鹽誘導(dǎo)激酶2，是一種絲氨酸/蘇氨酸蛋白激酶，屬于AMP活化蛋白激酶(adenosine monophosphate-activated protein kinase， AMPK)亞家族。研究發(fā)現(xiàn)SIK2基因在脂肪組織中的表達(dá)最高，且在脂肪細(xì)胞中過表達(dá)SIK2時，可以顯著降低細(xì)胞核類固醇調(diào)節(jié)元件結(jié)合蛋白1(sterol regulatory element binding protein-1，SREBP-1)水平，進(jìn)而抑制多個參與脂肪合成代謝調(diào)節(jié)基因的表達(dá)，如乙酰輔酶A羧化酶2基因(ACC2)和脂肪酸合成酶基因(FAS)[25-26]。也有研究表明，SIK2可以通過磷酸化過程參與胰島素和胰高血糖素對血糖濃度的調(diào)節(jié)過程中[27]。綜上，SIK2基因在脂質(zhì)代謝和糖代謝中發(fā)揮重要的作用，該基因可能是通過參與到調(diào)節(jié)肌肉中脂肪的生長，進(jìn)而影響肌肉的嫩度，可以作為調(diào)節(jié)白山羊肌肉生長發(fā)育的候選基因。

胰島素樣生長因子1受體(insulin-like growth factor 1 receptor，IGF1R)是一種表面受體蛋白，研究表明，IGF1R主要通過介導(dǎo)IGF1和IGF2信號參與調(diào)節(jié)多種組織及器官的生長發(fā)育，如IGF1R基因缺失會導(dǎo)致小鼠胚胎死亡或肌肉發(fā)育遲緩，而過表達(dá)則會導(dǎo)致肌肉組織中肌肉發(fā)育肥大[28]。有研究通過抑制IGF1R表達(dá)可以減弱肌肉的再生能力及細(xì)胞的增殖[29-30]。FüRSTENBERGER等[31]研究發(fā)現(xiàn)IGF1基因具有可以調(diào)節(jié)細(xì)胞的增殖分化和細(xì)胞凋亡的作用，而IGF1R是IGF1發(fā)揮具體功能的關(guān)鍵受體。而KLOTING等[32]通過降低小鼠脂肪組織中IGF1R的活性可以增加脂肪質(zhì)量。YANG等[33]研究發(fā)現(xiàn)，IGF1R可以通過AKT信號通路促進(jìn)肌肉的分化。在畜牧生產(chǎn)中，占思遠(yuǎn)等[34]研究發(fā)現(xiàn)IGF1R基因在山羊肌肉組織和肌肉細(xì)胞中均為高表達(dá)基因，且IGF1R與IGF2R可以協(xié)同促進(jìn)肌肉組織發(fā)育和肌細(xì)胞增殖。馬懿磊等[35]研究表明IGF1R基因在胎牛和成年牛時期的肌肉中的表達(dá)存在差異，且干擾IGF1R基因能夠抑制牛肌細(xì)胞的分化。綜上，IGF1R基因具有調(diào)節(jié)肌肉生長發(fā)育的功能，是研究肌肉生長發(fā)育的候選基因。

綜上，本研究基于RNA-Seq技術(shù)篩選出貴州白山羊公羊背最長肌的差異下調(diào)基因514個，其中100個為新發(fā)掘基因或轉(zhuǎn)錄本，發(fā)現(xiàn)參與核糖體信號通路的基因最為富集，為27個基因。通過功能分析及RT-PCR驗證，初步認(rèn)為RNA-Seq技術(shù)篩選出的相關(guān)下調(diào)基因數(shù)據(jù)真實可靠，可作為貴州白山羊公羊背最長肌生長發(fā)育相關(guān)的重要候選基因，對進(jìn)一步研究白山羊肌肉生長發(fā)育的遺傳機(jī)制提供了依據(jù)。