基于轉(zhuǎn)錄組方法鑒定貴州白山羊肌肉生長發(fā)育相關(guān)上調(diào)基因

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

（銅仁學(xué)院，貴州省梵凈山地區(qū)生物多樣性保護(hù)與利用重點(diǎn)實(shí)驗(yàn)室，貴州銅仁 554300）

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

地方品種山羊相較于商品羊生產(chǎn)性能較低，但對當(dāng)?shù)氐沫h(huán)境適應(yīng)性好，尤其對地方特定的疾病具有一定耐受性，是發(fā)展地方畜牧業(yè)的重要種質(zhì)資源。貴州白山羊是在云貴高原特殊的生態(tài)環(huán)境下，經(jīng)過長期自然選擇和人工選擇形成的特有品種，具有耐粗飼、性成熟早、繁殖力強(qiáng)、適應(yīng)性強(qiáng)、肉質(zhì)細(xì)嫩、無膻味等優(yōu)點(diǎn)，但同時(shí)存在個(gè)體較小、生長速度慢、產(chǎn)肉量不高等缺點(diǎn)。因此，本研究選用貴州白山羊?yàn)檠芯繉ο?，通過屠宰實(shí)驗(yàn)分析不同體重山羊的屠宰情況，結(jié)合轉(zhuǎn)錄組測序技術(shù)挖掘影響山羊生長性狀的重要基因，探索影響貴州白山羊生長發(fā)育的基因調(diào)控網(wǎng)絡(luò)，以期為肉用山羊生長發(fā)育提供重要的參考數(shù)據(jù)。

1 材料與方法

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

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

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

1.4 引物合成與qRT-PCR 驗(yàn)證 以NCBI 數(shù)據(jù)庫中山羊FGF11、ADIPOQ、NOTCH1、HBEGF基因的mRNA序列為參考序列，以GAPDH基因?yàn)閮?nèi)參基因，應(yīng)用Primer Premier 5.0 軟件設(shè)計(jì)相應(yīng)引物，由生工生物工程（上海）股份有限公司合成，引物序列見表1。

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

qRT-PCR 驗(yàn)證差異表達(dá)mRNA 相對表達(dá)水平，采用3 個(gè)樣本重復(fù)，應(yīng)用北京全式金生物技術(shù)有限公司生產(chǎn)的熒光定量試劑盒進(jìn)行qPCR 擴(kuò)增分析，擴(kuò)增體系為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 反應(yīng)程序?yàn)椋?4℃預(yù)變性1 min；94℃10 s，59℃30 s，72℃10 s，46 個(gè)循環(huán)。使用2-ΔΔCT法計(jì)算各基因間的相對表達(dá)水平。

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

2 結(jié)果與分析

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

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

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

表3 貴州白山羊背最長肌總RNA 轉(zhuǎn)錄組測序數(shù)據(jù)情況

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

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

表4 優(yōu)勢體重群體背最長肌中表達(dá)量最高的10 個(gè)基因

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

表5 貴州白山羊背最長肌中上調(diào)表達(dá)基因GO 分析

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

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

2.5 差異表達(dá)基因qRT-PCR 驗(yàn)證 通過功能分析，隨機(jī)選取4 個(gè)可能與貴州白山羊肌肉生長發(fā)育呈上調(diào)趨勢的調(diào)控基因FGF11、NOTCH1、HBEGF和ADIPOQ進(jìn)行qRT-PCR 驗(yàn)證，具體基因信息見表6。通過qRTPCR 驗(yàn)證，選取的4 個(gè)候選基因相對表達(dá)結(jié)果與轉(zhuǎn)錄組分析結(jié)果一致，說明測序數(shù)據(jù)準(zhǔn)確可靠，結(jié)果見圖3。

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

表6 可能與貴州白山羊背最長肌生長發(fā)育相關(guān)的候選上調(diào)基因

3 討 論

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

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

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

4 結(jié) 論

本研究基于轉(zhuǎn)錄組測序技術(shù)對不同體重的貴州白山羊背最長肌進(jìn)行研究，在體重優(yōu)勢群體中發(fā)現(xiàn)563 個(gè)顯著上調(diào)基因，其中94 個(gè)為新轉(zhuǎn)錄本，主要參與PI3KAkt 等信號通路的調(diào)控，篩選出的顯著上調(diào)基因與轉(zhuǎn)錄組測序預(yù)測基因差異表達(dá)結(jié)果一致，表明差異基因結(jié)果可靠。因此，篩選出的上調(diào)基因可作為提升山羊群體肌肉生長性能的重要候選基因，值得進(jìn)一步探討其對貴州白山羊生長發(fā)育的調(diào)控機(jī)制。