基于轉錄組方法鑒定貴州白山羊肌肉生長發育相關上調基因

安清明，王 星，宋興超，孟金柱，趙園園，吳震洋

（銅仁學院，貴州省梵凈山地區生物多樣性保護與利用重點實驗室，貴州銅仁 554300）

羊肉因其獨特的風味和較高的營養價值受到消費者喜愛。近年來，我國羊肉需求量逐年增長。國家統計局統計數據顯示，2020 年我國羊肉產量為492.31 萬t，較2012 年的404.50 萬t 增加了近90 萬t，但人均年占有羊肉量只有3.52 kg，較2012 年的人均3.32 kg 幾乎沒有增長，表明我國羊肉產量仍存在巨大的缺口。如何提高肉羊產肉量仍是育種和養殖工作的重點關注方向。研究發現，動物肌肉組織的生長發育狀況直接影響動物總體重，通常認為肌肉的重量可占哺乳動物胴體重的50% 以上[1]。研究認為，動物大多數與肌肉性狀相關的因素都具有顯著的遺傳效應，如骨骼肌生長發育、肌內脂肪含量、嫩度等指標[2-4]，一定程度上均受相應基因的調控。如何高效準確的篩選出與動物肌肉生長發育相關的基因是分子標記育種的一個難點，但隨著測序技術的不斷發展，轉錄組測序技術被廣泛應用到畜禽重要性狀基因篩選和分子機制研究中。如劉瑞莉等[5]基于轉錄組測序技術篩選布萊凱特黑牛和魯西黃牛背最長肌差異基因，篩選出296 個上調基因和258 個下調基因，最終推測MYL2、MYL3、MYLPF基因可能在肉牛骨骼肌生長發育過程中發揮重要作用。徐盼等[6]對蘇姜豬不同部位、不同生長時期骨骼肌進行轉錄組測序，共發現1 655 個差異表達基因，最終篩選鑒定FOXO1、FOXO3、FOXO4、MYF6、A-FABP、H-FABP基因可作為改善豬肉質性狀的候選基因。在綿羊育種中，也有研究通過轉錄組測序技術篩選了部分生產性狀相關的候選基因，如史金平等[7]對雜交綿羊群體進行轉錄組測序研究，挖掘出397 個差異表達基因，篩選出ADIPOQ、IGFBP7、ACTC1基因作為肉質性狀相關的候選基因。在山羊生長性能研究中，曹艷紅等[8]通過對努比亞山羊和隆林山羊進行轉錄組測序發現ACACA、MFGE8、MYH3、ACTN2基因是山羊產肉性狀的關鍵基因，可能在肉質脂肪含量、肌肉生長發育等方面發揮重要作用。孟憲然等[9]篩選發現ADIPOQ、PDK4、CD36基因可作為改良絨山羊肉品質的候選基因。上述研究表明，從遺傳學角度分析畜禽生產相關性狀已十分普遍，但目前針對貴州白山羊肉用性能的研究仍然較少，仍需進一步深入研究。

地方品種山羊相較于商品羊生產性能較低，但對當地的環境適應性好，尤其對地方特定的疾病具有一定耐受性，是發展地方畜牧業的重要種質資源。貴州白山羊是在云貴高原特殊的生態環境下，經過長期自然選擇和人工選擇形成的特有品種，具有耐粗飼、性成熟早、繁殖力強、適應性強、肉質細嫩、無膻味等優點，但同時存在個體較小、生長速度慢、產肉量不高等缺點。因此，本研究選用貴州白山羊為研究對象，通過屠宰實驗分析不同體重山羊的屠宰情況，結合轉錄組測序技術挖掘影響山羊生長性狀的重要基因，探索影響貴州白山羊生長發育的基因調控網絡，以期為肉用山羊生長發育提供重要的參考數據。

1 材料與方法

1.1 實驗動物 受試貴州白山羊飼養于貴州省銅仁市沿河縣華珍牧業有限責任公司養殖試驗基地，所選山羊為同等飼養水平下采用半舍飼加放牧方式進行飼養，選取同一種群中6 只2 周歲齡健康、體重不同的貴州白山羊（優勢群體3 只，體重較高；劣勢群體3 只，體重較輕；組內活體重相近），按照國家實驗動物屠宰標準進行屠宰并測定屠宰性能（宰前活重、胴體重、屠宰率、凈肉重和凈肉率），采集背最長肌組織樣，迅速放入液氮中，帶回實驗室轉入-80℃冰箱保存，用于轉錄組測序和qRT-PCR。

1.2 總RNA 提取cDNA 文庫建立 背最長肌中提取總RNA，利用Agilent Bioanalyzer 2100 檢測完整性，Qubit 2.0 Flurometer 測量RNA 濃度，RNA 樣品檢查合格后通過Illumina Hiseq 2500 平臺完成測序，測序分析由北京諾禾致源生物信息科技有限公司完成。

1.3 測序數據處理 將測序得到的原始數據進行評估，去除含有帶接頭和低質量的片段，利用TopHat 2 軟件與已知山羊參考基因組（GCF_001704415.1）進行比對，獲取測序拼裝序列相應特征信息，利用Cufflinks 軟件對轉錄本數據組裝、注釋和合并，利用DESeq2 軟件對mRNA 進行差異表達分析，利用基因本體GO 數據庫對篩選出的差異表達基因進行功能注釋，利用KEGG 數據庫進行DEGs 的信號通路富集分析。

1.4 引物合成與qRT-PCR 驗證 以NCBI 數據庫中山羊FGF11、ADIPOQ、NOTCH1、HBEGF基因的mRNA序列為參考序列，以GAPDH基因為內參基因，應用Primer Premier 5.0 軟件設計相應引物，由生工生物工程（上海）股份有限公司合成，引物序列見表1。

表1 熒光定量引物基因列表

qRT-PCR 驗證差異表達mRNA 相對表達水平，采用3 個樣本重復，應用北京全式金生物技術有限公司生產的熒光定量試劑盒進行qPCR 擴增分析，擴增體系為20 μL，其 中cDNA模板1 μL（1 μg/L），2×qPCR super mix 10 μL，上、下游引物各0.5 μL（10 μmol/L），dd H2O 8 μL。qPCR 反應程序為：94℃預變性1 min；94℃10 s，59℃30 s，72℃10 s，46 個循環。使用2-ΔΔCT法計算各基因間的相對表達水平。

1.5 數據統計與分析 屠宰數據經Excel 整理后，采用SPSS 22.0 軟件對6 只2 周歲齡健康、體重不同的貴州白山羊宰前活重、胴體重、屠宰率、凈肉重和凈肉率進行方差分析。表達量分析ΔΔCt=[（Ct目的基因-CtGAPDH）]實驗組-[（Ct目的基因-CtGAPDH）]對照組。結果應用SPSS 22.0軟件One-way ANOVA 進行顯著性檢驗分析，P＜0.05 時，差異顯著。

2 結果與分析

2.1 貴州白山羊屠宰性能差異分析 實驗測定不同體重貴州白山羊2 歲齡屠宰性能（AG 為體重優勢群體，DG為體重劣勢群體），優勢群體的活體重、胴體重和凈肉重顯著高于劣勢群體，表明優勢群體屠宰性能高于劣勢群體（表2），也說明后續測序試驗樣品具有可靠性。

表2 貴州白山羊屠宰性能差異分析

2.2 RNA-seq 數據比對分析 由表3 可知，最終6 個樣本共獲得78.99 Gb 的有效數據，單個樣品的有效數據均在12Gb 以上，Q20（有效數據質量不小于20 的堿基所占百分比）均在95%以上，Q30（有效數據質量不小于30 的堿基所占百分比）均在89%以上，該結果表明測序數據質量可靠，可用于進一步分析。

表3 貴州白山羊背最長肌總RNA 轉錄組測序數據情況

2.3 差異表達基因的篩選、功能注釋及富集分析 FPKM是衡量轉錄組測序中基因表達水平的估算指標。在6個樣本中共獲得33 896 個轉錄本，篩選出FPKM ≥1的轉錄本有25 089 個。表4 列出在體重優勢群體背最長肌中表達量前10 的候選基因。應用DESeq2 軟件對篩選出的基因進行差異表達篩選，當設定參數FPKM ≥1 且P值小于0.05 時，在優勢群體背最長肌中共篩選出563 個差異表達上調基因，其中94 個為新發現基因（轉錄本）。皮爾遜相關分析結果顯示，2 組3個不同的重復樣品的R2值基本都在0.85 以上（R2值越高，表明樣品可靠性和可重復性越高），表明實驗樣品具有較高的可靠性（圖1A），聚類分析熱圖顯示2 個不同群體間差異表達基因重復性較好，差異基因分別聚類在一起（圖1B）。

圖1 貴州白山羊不同群體背最長肌樣品及差異基因聚類分析

表4 優勢體重群體背最長肌中表達量最高的10 個基因

2.4 差異表達基因GO 功能注釋及KEGG 富集分析 GO注釋系統可分為3 大類，分別為：生物學過程、細胞組分和分子功能。通過Goseq 軟件對篩選出的差異表達上調基因進行GO 功能富集分析，其中36 個為顯著富集（P＜0.05），這些基因在3 個分支中均有分布，包括9個生物學過程，20 個細胞組分和7 個分子功能（表5）。

表5 貴州白山羊背最長肌中上調表達基因GO 分析

應用Kobas 軟件對差異表達基因進行KEGG 通路分析，結果表明注釋到的差異表達上調基因共參與220條信號通路，其中FDR ≤0.05 的通路為表達基因中顯著富集的信號通路，共有21 條，圖2 列出了前20 個顯著富集到的通路，最為富集的前5 條通路分別為代謝通路途徑（含64 基因）、磷脂酰肌醇3-激酶信號通路（含25 個基因）、阿爾茨海默?。ê?4 個基因）、致癌的途徑（含21 個基因）、帕金森?。ê?0 個基因）和非酒精性脂肪肝（含20 個基因）。

圖2 貴州白山羊背最長肌中上調表達基因KEGG 通路富集散點圖

2.5 差異表達基因qRT-PCR 驗證 通過功能分析，隨機選取4 個可能與貴州白山羊肌肉生長發育呈上調趨勢的調控基因FGF11、NOTCH1、HBEGF和ADIPOQ進行qRT-PCR 驗證，具體基因信息見表6。通過qRTPCR 驗證，選取的4 個候選基因相對表達結果與轉錄組分析結果一致，說明測序數據準確可靠，結果見圖3。

圖3 候選基因在不同群體貴州白山羊背最長肌中的相對表達量

表6 可能與貴州白山羊背最長肌生長發育相關的候選上調基因

3 討 論

貴州白山羊是在云貴高原特有的生態環境下，經過長期自然選擇和人工選擇而形成的地方品種，具有耐粗飼、性成熟早、繁殖力強、適應性強、肉質細嫩、無膻味等優點，但同時存在個體較小、生長速度慢、產肉量不高等缺點。近年來，由于自然和人為因素的影響，貴州白山羊個體生產性能差異較大，品種性能退化嚴重，并呈進一步加速趨勢。因此，研究可用于輔助篩選貴州白山羊肌肉生長發育的候選基因及調控機制，對進一步提升貴州白山羊生產性能具有重要指導意義。肌肉生長發育是影響山羊生產性能的重要指標，本研究對同等飼養條件下、不同體重白山羊的背最長肌進行全轉錄組測序，篩選出優勢體重群體相對于劣勢體重群體的差異上調基因563 個，其中94 個為新發現基因（轉錄本），這些差異基因可能與貴州白山羊肌肉生長發育存在一定的關聯性，值得進一步研究。

KEGG 富集分析得到了21 條富集的信號通路，最為富集的5 條通路中，PI3K-Akt 信號通路在細胞增殖、生長、代謝和生存中起著重要調控作用[10]，其中磷脂酰肌醇3 激酶I 類中的p85α亞基可以介導胰島素的調控進而影響葡萄糖轉運[11]。通路中蛋白激酶B（Akt/Pkb）重要的亞型Akt2 主要在骨骼肌、棕色脂肪中表達，而且研究發現，活化后的Akt 可以激活上百種底物，可以調控細胞生長、增殖、代謝、運動等功能[12-13]。呂欣等[14]也認為通過影響PI3K/Akt 信號通路中的蛋白表達、誘導通路的磷酸化，可以促進肌衛星細胞增殖分化、抑制凋亡，影響骨骼肌的再生修復。另外，篩選出基因最多的代謝途徑信號通路中，通路信號眾多，主要包括多聚糖生物合成代謝通路、脂類代謝通路、碳水化合物代謝通路、能量代謝通路、氨基酸代謝通路、核酸代謝通路、萜類化合物和多酮類化合物代謝通路等，均是參與能量代謝調控的通路。因此，本研究篩選出的通路可能在調解貴州白山羊肌肉生長發育中具有重要作用。

在篩選出的差異表達基因中，已有研究表明部分基因與肌肉的生長發育相關，本研究隨機篩選4 個基因進行了驗證。其中FGF11基因編碼蛋白是一種自分泌型成纖維細胞因子，主要通過電壓門控的鈉離子通道相互作用調控神經系統的興奮性[15]。已有研究表明，FGF11與毛細血管的結構穩定性及形成有關，同時還可能是腫瘤的重要調控因子[16]。Tejedor 等[17]最新研究發現FGF11基因敲除能夠阻礙小鼠前肢的再生發育，表明FGF11基因可能是動物肢體發育的重要調控因子。本研究中發現，FGF11基因在體重優勢群體中的表達量遠高于體重劣勢群體。NOTCH1基因編碼蛋白是NOTCH信號通路中的重要受體，研究發現NOTCH 信號通路可以調解骨骼肌細胞的穩態、脂肪組織的相互轉變、較少肥胖等[18]。亦有研究分析認為NOTCH1可以調解下游的PI3K、Akt 信號通路，而PI3K 和Akt 信號通路是參與調控肌肉細胞增殖、分化和發育的重要通路，同時該通路還參與葡萄糖的攝取與儲存[19-20]。Yamaguchi 等[21]研究發現NOTCH1在小鼠體內的合成不足可以加速體內脂肪的生成。HBEGF基因編碼蛋白是EGF 家族的糖蛋白成員，最早在巨噬細胞培養后的產物中獲得，研究發現HBEGF可以促進平滑肌細胞、成纖維細胞和上皮細胞的增殖和遷移[22-23]。亦有研究表明，HBEGF基因是影響豬繁殖性狀的重要候選基因[24-25]。本研究中發現HBEGF基因在體重優勢群體中的表達量雖高于體重劣勢群體，但差異不明顯。因此對該基因的具體功能需進一步研究。ADIPOQ基因編碼蛋白是一種對肥胖呈負相關調控的蛋白類激素，其主要參與動物機體內葡萄糖代謝和脂肪代謝。研究表明，ADIPOQ基因在脂肪組織中的表達量和在血液中的水平會顯著升高，而當脂肪過度沉積時，在表達量和在血液中的水平均降低[26]。研究還發現ADIPOQ基因表達量與豬的胴體性狀、肉質[27]，牛的眼肌面積和背膘厚度顯著相關[28]，本研究團隊也曾發現ADIPOQ基因與綿羊的生長速度和胴體性狀相關[29]。羅宗剛等[30]研究發現ADIPOQ基因的表達在動物個體背最長肌中存在性別差異，表明ADIPOQ基因可能參與調解不同性別個體肌肉的生長調控。還有研究發現ADIPOQ基因在山羊脂肪組織中的表達量具有發育階段的特異性[31]。在本研究中，不同群體間ADIPOQ基因的表達也存在顯著差異，因此，ADIPOQ基因可能是肌肉生長發育的重要候選基因，其功能值得深入研究。

4 結 論

本研究基于轉錄組測序技術對不同體重的貴州白山羊背最長肌進行研究，在體重優勢群體中發現563 個顯著上調基因，其中94 個為新轉錄本，主要參與PI3KAkt 等信號通路的調控，篩選出的顯著上調基因與轉錄組測序預測基因差異表達結果一致，表明差異基因結果可靠。因此，篩選出的上調基因可作為提升山羊群體肌肉生長性能的重要候選基因，值得進一步探討其對貴州白山羊生長發育的調控機制。