李 潔，趙瑞鵬，陳楚雯，楊朝武，吳錦波，李志雄*

(1.西南民族大學(xué) 青藏高原動物遺傳資源保護與利用教育部重點實驗室，成都 610041；2.西南民族大學(xué)畜牧獸醫(yī)學(xué)院，成都 610041； 3.西南民族大學(xué)青藏高原研究院，成都 610041；4.四川省畜牧科學(xué)研究院，成都 610066； 5.阿壩藏族羌族自治州畜牧科學(xué)技術(shù)研究所，阿壩 624402)

藏雞是我國特有的雞種，于2020年5月29日被列入《國家畜禽遺傳資源品種名錄》，主要分布在西藏、青海等海拔較高的地區(qū)，是當(dāng)?shù)胤浅Ｖ匾男笄萜贩N之一。該品種具有抗病力強且耐粗飼的優(yōu)點，其肉質(zhì)鮮美，肌肉中的必需氨基酸和鮮味氨基酸含量豐富[1]。經(jīng)過長期的自然選擇，藏雞表現(xiàn)出對高海拔環(huán)境穩(wěn)定的遺傳適應(yīng)性[2]。因此，藏雞是發(fā)展高原養(yǎng)殖業(yè)不可替代的品種資源。雖然藏雞肉品質(zhì)好，但生長速度慢，與普通的肉雞相比體重和體尺較小，產(chǎn)肉率和屠宰性能不佳。大恒肉雞是由四川省畜牧科學(xué)研究院培育的肉雞新品種，該品種不僅保持了地方雞種的優(yōu)質(zhì)風(fēng)味和外觀性狀，且生長速度快，生產(chǎn)性能和屠宰率高[3-4]。根據(jù)2018年的統(tǒng)計顯示，商品代大恒肉雞70日齡公雞和母雞的體重分別為2.32和1.79 kg，而成年藏雞公雞和母雞的體重分別僅為1.44和1.15 kg[5]。禽體的肌肉約占體重的30%～40%，肌纖維是構(gòu)成肌肉的基本單位，肌纖維的數(shù)量在雛雞孵化前就已基本固定。在雞胚胎期肌肉發(fā)育的兩個主要階段中，對肌肉生長發(fā)育有決定作用的肌細胞在胚胎期間不可逆的退出細胞周期并不再增殖[6]。除此之外，在胚胎末期，成熟的肌纖維周圍會出現(xiàn)參與肌纖維肥大和成熟的肌衛(wèi)星細胞[7]，因此，比較分析藏雞和大恒肉雞這兩個品種在胚胎期骨骼肌生長發(fā)育上的差異對解決藏雞生長速度慢，產(chǎn)肉率低等問題具有重要的意義。

影響動物肌肉生長發(fā)育的因素有很多，例如：基因、營養(yǎng)、飼養(yǎng)管理、環(huán)境、品種以及年齡等。其中基因的影響是最根本、最具有決定性的。近年來，隨著分子生物學(xué)的快速發(fā)展以及生物信息技術(shù)的不斷進步，基因在調(diào)控畜禽經(jīng)濟性狀中的功能以及機制不斷被挖掘和破解。例如：肌分化因子MyoD1[8]通過激活肌肉特異性基因來調(diào)控肌肉的生長發(fā)育、生肌調(diào)節(jié)因子(myogenic regulatory factors，MRFs)家族，如肌細胞生成素MyoG參與調(diào)控骨骼肌細胞的增殖、分化以及融合[9]，MyoG、MyoD和MRF4是成肌細胞分化所必需的[10]、肌肉骨骼胚胎核蛋白1(MUSTN1)基因[11]與肌肉的分化和再生有關(guān)等，這些都是在肌肉生長發(fā)育過程中起重要作用的基因。

lncRNA是長度大于200 nt的非編碼RNA[12]。迄今為止，已有大量的研究證明lncRNA可以通過干擾基因的表達、與蛋白質(zhì)相結(jié)合、作為競爭性內(nèi)源RNA(competitive endogenous RNA，ceRNA)以及編碼肽類等多種作用方式直接或間接的參與家畜生命活動的幾乎每一個過程[13-16]。在畜禽肌肉的生長發(fā)育的過程中，lncRNA常作為ceRNA與肌肉生長發(fā)育相關(guān)基因競爭結(jié)合同一個miRNA從而間接調(diào)控肌肉生長發(fā)育相關(guān)基因的表達，最終影響畜禽肌肉的生長發(fā)育，參與調(diào)控。也有越來越多特定的lncRNA被發(fā)現(xiàn)參與調(diào)控畜禽肌肉的生長發(fā)育，例如有研究表明lncRNA MEG3參與豬骨骼肌分化的調(diào)控，并且作為提高豬肉產(chǎn)量的候選基因[17]，lncRNA-Six1能夠通過編碼微肽間接的促進雞骨骼肌細胞的增殖[18]，lncMD1可作為ceRNA與miR-133和miR-135相結(jié)合，進而調(diào)節(jié)與晚期肌細胞增殖分化相關(guān)的轉(zhuǎn)錄因子表達，從而調(diào)控牛肌肉的生長發(fā)育[15]，以及l(fā)nc-002783能夠抑制山羊骨骼肌衛(wèi)星細胞的增殖和分化[19]。大量的研究證實，lncRNA在很多畜禽經(jīng)濟性狀的調(diào)控中普遍存在，并發(fā)揮著重要的作用。但lncRNA在藏雞肌肉生長發(fā)育中的相關(guān)研究卻很少。因此，本研究通過轉(zhuǎn)錄組測序篩選了胚胎期藏雞和大恒肉雞腿肌中差異表達的mRNAs和lncRNAs，挑選與肌肉生長發(fā)育相關(guān)的mRNAs和lncRNAs構(gòu)建競爭性調(diào)控網(wǎng)絡(luò)，為后續(xù)深入分析兩雞種在肌肉生長發(fā)育方面呈現(xiàn)差異的分子機制奠定基礎(chǔ)。

1 材料與方法

1.1 試驗動物及樣品采集

本試驗所用藏雞和大恒肉雞的種蛋分別由阿壩州九頂原生態(tài)畜禽養(yǎng)殖有限責(zé)任公司和四川大恒家禽育種有限公司(國家肉雞核心育種場)提供。將所有種蛋按照常規(guī)程序在相同環(huán)境條件下進行孵化，胚胎期第18天，在兩個品種的種蛋中分別挑選3個發(fā)育正常的雄性雞胚，快速采集腿肌組織并裝入無RNA酶的凍存管中，藏雞的腿肌組織標(biāo)記為TC_E18，大恒肉雞標(biāo)記為DH_E18，之后迅速置于液氮中保存?zhèn)溆谩?/p>

1.2 文庫構(gòu)建及轉(zhuǎn)錄組測序

按照TRIzol試劑盒(Invitrogen)說明書提取藏雞和大恒肉雞腿肌組織的總RNA。之后通過1%瓊脂糖凝膠電泳檢測RNA的質(zhì)量，利用NanoDrop ND-2000(Agilent)檢測RNA的濃度和純度。將檢測合格后的總RNA使用PrimeScriptTMRT reagent Kit with gDNA Eraser(Perfect Real Time)反轉(zhuǎn)錄合成cDNA鏈，置于-80 ℃保存?zhèn)溆?。兩個品種采用去除核糖體RNA的方法共構(gòu)建6個特異性文庫并進行轉(zhuǎn)錄組測序(轉(zhuǎn)錄組文庫的構(gòu)建和測序由上海美吉生物醫(yī)藥科技有限公司完成)。

1.3 測序數(shù)據(jù)質(zhì)控及序列比對

通過Illumina二代測序技術(shù)進行測序后產(chǎn)生大量以fastq格式儲存起來的原始數(shù)據(jù)(raw data)，對每個樣本的原始測序數(shù)據(jù)進行相關(guān)質(zhì)量評估。通過使用fastp軟件對堿基質(zhì)量分布、堿基錯誤率分布以及堿基含量分布進行統(tǒng)計以保證之后生物信息分析的準確可靠性。由于原始數(shù)據(jù)中會包含測序接頭序列、低質(zhì)量讀段、不確定堿基信息(N率)較高的序列以及長度過短序列等，這些都會嚴重影響后續(xù)分析的準確性。因此統(tǒng)計結(jié)果符合要求之后，使用SeqPrep(https://github.com/jstjohn/SeqPrep)和Sickle(https://github.com/najoshi/sickle)軟件對原始序列進行過濾。去除接頭序列、由于接頭自連等原因?qū)е碌臎]有插入片段的reads、含有不確定堿基的reads以及修剪后長度較小(< 30 bp)的序列，剔除質(zhì)量較差序列(將序列3′末端質(zhì)量值<20的堿基修減后，若剩余序列中仍然有質(zhì)量值< 10的堿基，則將整條序列剔除)。由此得到質(zhì)量較高的測序數(shù)據(jù)(clean data/clean reads)。

將質(zhì)控后得到的Clean reads利用Hisat(http://ccb.jhu.edu/software/hista2/index.shtml)軟件比對到雞參考基因組上(http://asia.ensembl.org/Gallus_gallus/Info/Index)[20]，參考基因組版本為GRCg6a。之后將比對結(jié)果利用StringTie(http://ccb.jhu.edu/software/stringtie/)軟件進行轉(zhuǎn)錄本組裝獲得每個樣本的轉(zhuǎn)錄本信息以及發(fā)現(xiàn)新的轉(zhuǎn)錄本或新基因。

1.4 差異表達mRNA和lncRNA的篩選及功能分析

以大恒肉雞為對照組，藏雞為處理組，使用RSEM軟件(http://deweylab.github.io/RSEM/)[21]通過單端或雙端測序數(shù)據(jù)來計算基因或者轉(zhuǎn)錄本的表達量，并以TPM(Transcripts Per Million reads)作為表達衡量指標(biāo)。使用DESeq2、DEGseq軟件以及備選軟件edgeR檢測差異表達的mRNAs和lncRNAs，分別以Padjust<0.05、Padjust<0.001(Padjust為多重檢驗校正后的P值)和|log2FC|≥1作為篩選差異表達mRNAs和lncRNAs的條件。

為進一步分析大恒肉雞和藏雞腿肌中差異表達mRNAs和lncRNAs的功能，將篩選得到的差異表達mRNAs(differentially expressed mRNAs，DEMs)進行GO(http://www.geneontology.org/)富集分析和KEGG(Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)通路分析，以Padjust<0.05為標(biāo)準，篩選顯著富集的GO條目和KEGG通路。對于篩選得到的差異表達lncRNAs(differentially expressed lncRNAs，DELs)，首先根據(jù)lncRNA不同的作用模式進行順式(cis)作用靶基因預(yù)測和反式(trans)作用靶基因預(yù)測，然后將預(yù)測得到的DELs的靶基因進行GO富集分析和KEGG通路分析，最后篩選顯著富集的GO條目和KEGG通路。挑選得到與肌肉生長發(fā)育和脂質(zhì)代謝相關(guān)的GO條目和KEGG通路。

1.5 lncRNA-miRNA-mRNA競爭性調(diào)控網(wǎng)絡(luò)的構(gòu)建

收集富集到與肌肉生長發(fā)育和脂質(zhì)代謝相關(guān)的GO條目和KEGG通路上的基因。通過lncRNA靶基因預(yù)測，發(fā)現(xiàn)有些lncRNA可以同時作用于兩個基因；之后通過miRNA靶基因預(yù)測發(fā)現(xiàn)，有miRNA靶向于與肌肉生長發(fā)育和脂質(zhì)代謝相關(guān)的lncRNA和基因上；還有些lncRNA雖然沒有直接作用于與肌肉生長發(fā)育和脂質(zhì)代謝相關(guān)的基因上，但通過作用于miRNA可間接作用于相關(guān)基因。因此，通過lncRNA和miRNA靶基因預(yù)測，收集與肌肉生長發(fā)育和脂質(zhì)代謝相關(guān)的lncRNAs、miRNAs和基因整理其聯(lián)系。利用Cytoscape V-3.6.0軟件進行l(wèi)ncRNA-miRNA-mRNA競爭性調(diào)控網(wǎng)絡(luò)的構(gòu)建。

1.6 實時熒光定量PCR(qRT-PCR)驗證

為保證測序結(jié)果的準確可靠性，隨機選取基因和lncRNAs各5個，以GAPDH為內(nèi)參基因，使用SYBR Green法進行qRT-PCR。qRT-PCR驗證過程中所用到的lncRNAs、基因以及引物信息見表1。反應(yīng)的體系為：2×Taq Pro Universal SYBR qPCR Master Mix 10 μL，上、下游引物(10 μ mol·L-1)各0.4 μL，cDNA模板2 μL，ddH2O 7.2 μL，共20 μL。擴增程序為：95 ℃ 30 s；95 ℃ 10 s，60 ℃ 30 s，30個循環(huán)；熔解曲線程序：95 ℃ 15 s，60 ℃ 1 min，95 ℃ 15 s。每個樣本設(shè)置3個技術(shù)重復(fù)，采用2-ΔΔCt法計算基因和lncRNAs的相對表達量。

表1 實時熒光定量PCR引物信息

2 結(jié) 果

2.1 雞胚胎腿肌轉(zhuǎn)錄組測序數(shù)據(jù)統(tǒng)計

測序結(jié)果表明，原始數(shù)據(jù)經(jīng)質(zhì)控后，藏雞(TC)和大恒肉雞(DH)組分別獲得359 085 628和344 771 590條clean reads，兩個組中各個樣本Q30的百分比(測序質(zhì)量在99.9%以上的堿基占總堿基的百分數(shù))均在94%以上，測序堿基平均錯誤率(error rate)均在0.03%以下，見表2。表明測序的質(zhì)量高，數(shù)據(jù)可靠準確，可用于后續(xù)分析。

表2 測序數(shù)據(jù)質(zhì)量分析

2.2 雞胚胎腿肌轉(zhuǎn)錄組測序結(jié)果比對分析

將質(zhì)控后得到的高質(zhì)量讀段(clean reads)與雞參考基因組比對，結(jié)果顯示DH_E18和TC_E18組比對到雞參考基因組上的clean reads比對率為90.05%～94.35%，在參考序列上有多個比對位置的clean reads比率為12.69%～17.72%，比對效果較好，滿足后續(xù)分析要求(表3)。

表3 比對統(tǒng)計結(jié)果

2.3 雞胚胎腿肌中l(wèi)ncRNA的長度分布及分類統(tǒng)計

統(tǒng)計文庫中所有l(wèi)ncRNA長度分布，大多數(shù)lncRNA的長度超過2 000 bp(72.83%)，其次主要分布在1 001～1 500 bp之間(45.55%)(圖1A)。根據(jù)lncRNA相對于基因組上蛋白編碼基因的位置，lncRNA有sense exon overlap lncRNA，intergenic lncRNA，antisense lncRNA，sense intron overlap lncRNA以及bidirection lncRNA五種[22]。經(jīng)統(tǒng)計，文庫中大多數(shù)為intergenic lncRNA(72%)，其次是antisense lncRNA(13%)，sense intron overlap lncRNA的數(shù)量最少(1%)(圖1B)。

2.4 藏雞和大恒肉雞胚胎腿肌mRNA和lncRNA的差異表達分析

通過RSEM軟件對mRNA和lncRNA的表達量進行分析，發(fā)現(xiàn)共有10 460個mRNAs在藏雞和大恒肉雞中共同存在，有227個mRNAs只存在于大恒肉雞中，213個mRNAs只存在于藏雞中(圖2A)；在藏雞和大恒肉雞中都存在的lncRNAs有4 406個，只存在于大恒肉雞中的lncRNAs有1 315個，只存在于藏雞中的lncRNAs有1 347個(圖2B)。

圖1 藏雞和大恒肉雞腿肌lncRNA的長度分布 (A) 及分類統(tǒng)計 (B)Fig.1 Length distribution (A) and classification statistics (B) of lncRNA in leg muscle of Tibetan chicken and Daheng broiler

以DH-E18為對照組，TC-E18為處理組，差異表達分析結(jié)果顯示，差異表達的mRNAs有106個，其中上調(diào)的有48個，下調(diào)的有58個(圖3A)；得到差異表達的lncRNAs有28個，其中上調(diào)的有10個，下調(diào)的有18個(圖3B)。之后對DELs的靶基因進行預(yù)測，得到48個trans作用的靶基因。

圖2 藏雞和大恒肉雞腿肌mRNAs (A) 和lncRNAs (B) 表達量分析Fig.2 Expression analysis of mRNAs (A) and lncRNAs (B) in leg muscle of Tibetan chicken and Daheng broiler

圖3 差異表達mRNAs (A) 和lncRNAs (B) 火山圖Fig.3 Volcano chart of differentially experssed mRNAs (A) and lncRNAs (B)

2.5 藏雞和大恒肉雞胚胎腿肌DEMs和DELs的功能分析

DEMs的GO富集結(jié)果顯示，共有8個GO條目上出現(xiàn)顯著富集(圖4)，其中與骨骼肌生長發(fā)育有關(guān)的條目是骨骼肌細胞分化(GO： 0035914)，有較多DEMs富集在該條目。KEGG富集結(jié)果顯示，共有9個KEGG通路顯著富集(圖5)，其中與脂質(zhì)代謝相關(guān)的通路是類固醇生物合成信號通路(map00100)。

對DELs的靶基因進行KEGG富集分析，結(jié)果顯示顯著富集的KEGG通路共有10條(圖6)。富集到靶基因最多的KEGG通路是類固醇生物合成信號通路(map00100)，除此之外，與脂質(zhì)代謝相關(guān)的KEGG通路還有脂肪酸生物合成信號通路(map00061)。

圖4 DEMs的GO富集分析Fig.4 GO enrichment analysis of DEMs

2.6 lncRNA-miRNA-mRNA競爭性調(diào)控網(wǎng)絡(luò)分析

利用Cytoscape構(gòu)建lncRNA-miRNA-mRNA競爭性調(diào)控網(wǎng)絡(luò)，結(jié)果顯示，調(diào)控網(wǎng)絡(luò)共有由9個mRNAs、16個miRNAs以及5個lncRNAs組成的30個節(jié)點，產(chǎn)生了28個相互作用(圖7)。其中，一些lncRNAs既調(diào)控與骨骼肌細胞分化相關(guān)基因，又調(diào)控與脂質(zhì)代謝相關(guān)基因，如MSTRG.21490.5和ENSGALT00000098920；除此之外，一些lncRNAs還作用于調(diào)控雞骨骼肌細胞分化的miRNAs；同樣，與雞骨骼肌細胞分化和脂質(zhì)代謝相關(guān)的基因也同時受到多個lncRNAs和miRNAs的調(diào)控。

2.7 實時熒光定量PCR驗證

為驗證測序結(jié)果的準確性，結(jié)合富集分析以及競爭性調(diào)控網(wǎng)絡(luò)挑選DEMs和DELs各5個進行qRT-PCR。結(jié)果顯示，qRT-PCR中DEMs(圖8A)和DELs(圖8B)的上下調(diào)趨勢與測序結(jié)果一致。表明轉(zhuǎn)錄組測序結(jié)果準確可靠。

3 討 論

藏雞作為我國特有的地方品種，是西藏青海等高海拔地區(qū)發(fā)展高原養(yǎng)殖業(yè)的優(yōu)良雞種。關(guān)于藏雞的研究大多集中在產(chǎn)蛋性能以及對高海拔的適應(yīng)性上，有關(guān)肌肉生長發(fā)育以及肉品質(zhì)方面的研究還較少。本研究對胚胎第18天的藏雞和大恒肉雞腿肌進行全轉(zhuǎn)錄組測序，篩選DEMs和DELs，然后進一步對其進行功能分析，并建立競爭性調(diào)控網(wǎng)絡(luò)進行聯(lián)合分析，為進一步探索雞肌肉生長發(fā)育的分子機制提供一定的理論基礎(chǔ)。

圖5 DEMs的KEGG富集分析Fig.5 KEGG enrichment analysis of DEMs

圖6 DELs靶基因的KEGG富集分析Fig.6 KEGG enrichment analysis of DELs target genes

對篩選到的DEMs進行GO富集分析，顯著富集的條目主要是生物過程中的補體活化、免疫應(yīng)答、骨骼肌細胞分化以及細胞組分中的胞外空間。在KEGG富集分析中也有2條通路與免疫系統(tǒng)相關(guān)。補體活化系統(tǒng)是免疫應(yīng)答中固有免疫的重要組成部分，研究表明，補體活化途徑與骨骼的生長發(fā)育以及骨重塑之間有著密切的聯(lián)系[23]。補體C3是補體活化單位之一，Matsuoka等[24]研究表明，在核因子κB受體活化因子配體(RANKL)注射的高骨轉(zhuǎn)化狀態(tài)下，補體C3的表達量升高，在每天服用抑制補體C3激活的藥物狀態(tài)下，導(dǎo)致骨形成減弱，骨丟失加??；動物的產(chǎn)肉量和體型大小與骨骼的生長發(fā)育密切相關(guān)，成骨細胞和破骨細胞對于骨骼的形成、發(fā)育、生長以及修復(fù)具有重要的作用，而補體活化途徑對成骨細胞和破骨細胞的生長和分化具有重要的調(diào)節(jié)作用[22]，進而會影響骨骼的形成、發(fā)育和生長。提示可以從補體活化的角度進行深入研究，可能會在一定程度上提高藏雞的生長速度，增大其體格。近年來，越來越多的研究結(jié)果都表明，免疫應(yīng)答對肌肉的生長發(fā)育有一定的影響，Li等[25]對瀘寧雞和白羽肉雞腓腸肌的差異表達基因進行GO富集分析發(fā)現(xiàn)主要富集的條目為免疫應(yīng)答，表明免疫應(yīng)答對肌肉的生長發(fā)育具有一定的影響；Dou等[26]以大體型的武定雞和小體型的大圍山雞為樣本進行比較研究，發(fā)現(xiàn)在胸肌組織中，與大圍山雞相比，武定雞的免疫功能相關(guān)基因表達增加；袁茂等[27]對不同發(fā)育階段的藏雞腿肌的差異表達基因進行KEGG富集分析發(fā)現(xiàn)，在KEGG顯著富集的19條通路中有6條與免疫有關(guān)。骨骼肌細胞分化是直接與生長速度和產(chǎn)肉性狀相聯(lián)系的GO條目，富集于該條目的DEMs有生肌因子6(myogenic factor 6,MYF6)、Fos原癌基因(Fos proto-oncogene,FOS)、早期生長應(yīng)答蛋白基因(early growth response 1,EGR1)、錨蛋白重復(fù)序列和細胞信號抑制因子盒蛋白2基因(ankyrin repeat and suppressor of cytokine signalling box containing protein 2,ASB2)以及錨蛋白重復(fù)域1基因(ankyrin repeat domain 1,ANKRD1)。

圖中矩形節(jié)點代表lncRNAs；圓形節(jié)點代表miRNAs；三角形節(jié)點代表mRNAs，上方4個基因是與骨骼肌細胞分化相關(guān)基因，下方5個基因是與脂質(zhì)代謝相關(guān)基因In the diagram, square nodes represent lncRNAs; circle nodes represent miRNAs; triangular nodes represent mRNAs, which the top 4 genes are related to skeletal muscle cell differentiation and the bottom 5 genes are related to lipid metabolism圖7 lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)分析Fig.7 Regulatory network analysis of lncRNA-miRNA-mRNA

圖8 轉(zhuǎn)錄組測序中DEMs (A) 和DELs (B) 的qRT-PCR驗證Fig.8 Validation of DEMs (A) and DELs (B) in RNA-seq by qRT-PCR

其中，MYF6基因是生肌決定因子(mycgenic determination factor, MyoD)家族的成員[28]，對于肌肉的生長發(fā)育具有重要的調(diào)控作用。研究表明，MYF6基因參與肌纖維的成熟，并且維持肌纖維的分化狀態(tài)[29]，與MyoD家族中的其它成員協(xié)同調(diào)控，在肌肉生長發(fā)育的各階段發(fā)揮關(guān)鍵作用[30-32]。湯展毅等[33]通過克隆MYF6基因并構(gòu)建載體，運用脂質(zhì)體轉(zhuǎn)染法轉(zhuǎn)染魯西黃牛成纖維細胞，檢測MYF6對成肌纖維細胞的影響，發(fā)現(xiàn)與肌肉生長發(fā)育相關(guān)基因的表達量顯著升高，且轉(zhuǎn)染后的成纖維細胞發(fā)育為肌管，表明了MYF6具有促進肌細胞分化的功能；Li等[34]通過研究MYF6基因在家鴨胚胎和孵化后胸肌和腿肌中的表達模式證明，MYF6與肌肉生長發(fā)育有關(guān)，是骨骼肌生長發(fā)育的一個組成部分。FOS和EGR1基因?qū)儆诹⒃缁?IEGs)，在細胞的增殖分化中具有重要的作用[35-36]。FOS參與肌肉的生成[37]，李志雄等[38]研究發(fā)現(xiàn)，在注射了使肉雞腓腸肌纖維直徑顯著減小的miR-499-5p的處理組中FOS的表達量顯著高于對照組，提示FOS基因會抑制肌肉纖維的生長發(fā)育；但也有研究表明，F(xiàn)OS下調(diào)，通過調(diào)節(jié)TNF-α的表達進而激活NF-κB以及IGF-1信號通路的Akt酶從而抑制肌管的形成以及肌源性分化[39]，通過抑制酪氨酸磷酸化從而影響肌生成素的表達[40]，提示FOS可能促進肌肉的生長發(fā)育。表明肌肉的生長發(fā)育機制是一個極其復(fù)雜的過程。EGR1主要在肌肉中表達，已有研究表明，EGR1調(diào)控肌肉細胞的生長和分化[41]。Yang等[42]通過抑制靶向于EGR1的lncRNA-COX2，使成纖維細胞中的EGR1顯著增加，纖維生成顯著促進，并且成纖維細胞顯示出比對照組更多的遷移細胞，而在EGR1被抑制后，這些遷移被阻斷，說明EGR1基因具有促進成纖維細胞活化以及肌纖維生成的作用；崔亞鳳[43]利用CRISPRi技術(shù)激活和抑制EGR1，發(fā)現(xiàn)EGR1的表達量與小鼠成肌細胞C2C12的分化呈顯著的正相關(guān)，證明了EGR1能夠促進肌肉的分化。ASB2基因是ASB家族成員之一，有研究發(fā)現(xiàn)該基因可能參與骨骼肌的分化、參與肌肉再生以及肌肉損傷修復(fù)等過程[44]?？低ね45]在研究豬骨骼肌衛(wèi)星細胞時，發(fā)現(xiàn)ANKRD1在其分化過程中的表達量較高且在約克夏豬和監(jiān)利豬中的表達量差異較大，推測ANKRD1可能參與調(diào)控肌衛(wèi)星細胞的分化過程；除此之外，也有研究表明ANKRD1可能也參與調(diào)控成肌纖維細胞的生長[46]。這些基因都能夠促進肌肉的生長發(fā)育，但其表達量均顯著下調(diào)，由此提示相較于大恒肉雞來說，藏雞的生長速度較慢、產(chǎn)肉率低可能與這些基因的表達水平較低有關(guān)?？梢詫@些基因進行深入探索，揭示其具體的調(diào)控機制，為提高藏雞生長速度，改良其產(chǎn)肉性能的研究提供了思路。

lncRNA靶基因的KEGG顯著富集的結(jié)果中最可靠的是脂質(zhì)代謝相關(guān)通路，其中類固醇生物合成是最顯著富集的KEGG通路，該通路同樣出現(xiàn)在DEMs的KEGG顯著富集通路中。富集在該通路的基因分別是DHCR24、SQLE以及LSS，且ENSGALT00000098920(lnc-920)共同靶向于這3個基因。這些基因都與膽固醇的合成有關(guān)，分別是催化膽固醇合成的起始、最終形成以及膽固醇前體的重要因子[47-49]。Mirza等[50]敲除小鼠肌肉組織中的DHCR24后，敲除小鼠的體積較同窩中未敲除小鼠的體積小且在幾小時后死亡，表示DHCR24可能對肌肉的生長發(fā)育也有一定的影響。查詢差異表達量，發(fā)現(xiàn)這3個基因均顯著上調(diào)，而lnc-920顯著下調(diào)。猜測lnc-920的表達量可能與脂質(zhì)相關(guān)基因的表達量呈負相關(guān)，具體調(diào)控機制還需要進一步研究。因此，相較于大恒肉雞，藏雞的肉品質(zhì)較好，其脂質(zhì)代謝相關(guān)基因的表達量高可能是原因之一。同時值得注意的是，已有研究表明上述幾個基因可能與雞的脂肪沉積有關(guān)[47]，加之人們對畜禽產(chǎn)品需求的改變，提示藏雞作為高原養(yǎng)殖業(yè)的優(yōu)良雞種還有較大的遺傳改良潛質(zhì)。脂肪酸的含量及類型會影響肉的風(fēng)味，富集在脂肪酸生物合成通路上的基因?qū)?yīng)的lncRNA中MSTRG.21490.5(lnc-90.5)的表達量在藏雞和大恒肉雞腿肌中分別是其它DELs表達量的約7和20倍，推測該lncRNAs可能對肉質(zhì)有一定的影響，且該lnc-90.5是新發(fā)現(xiàn)的lncRNA，具有極大潛在的研究價值。

本研究基于DEMs和DELs靶基因富集分析得到的富集于骨骼肌生長分化以及脂質(zhì)代謝等條目上的基因，通過對相關(guān)lncRNAs以及miRNAs進行預(yù)測，構(gòu)建lncRNA-miRNA-mRNA競爭性調(diào)控網(wǎng)絡(luò)。近年來，許多研究發(fā)現(xiàn)lncRNA可以通過多種作用方式參與生物體的各種生物學(xué)過程。其中，lncRNA作為ceRNA與miRNA相結(jié)合進而發(fā)揮作用的報道已經(jīng)屢見不鮮。例如：Sun等[15]通過轉(zhuǎn)錄組測序研究秦川牛骨骼肌時發(fā)現(xiàn)lncMD1作為ceRNA與骨骼肌分化相關(guān)基因IGF2相競爭結(jié)合于miR-125b，最后提高了IGF2的表達，進而促進了骨骼肌的分化；巨曉軍[51]以清遠麻雞的胸大肌和縫匠肌為樣本，構(gòu)建lncRNA-miRNA-mRNA互作網(wǎng)絡(luò)，發(fā)現(xiàn)lncXR_003074785.1-miR-193-3p-PPARGC1A互作通路，其中PPARGC1A被證明調(diào)節(jié)雞骨骼肌纖維的形成以及肌細胞的分化，且miR-193-3p與PPARGC1A的表達呈負相關(guān)，暗示lncXR_003074785.1(lnc-85.1)可能作為ceRNA與miR-193-3p相結(jié)合從而影響肌肉的生長發(fā)育。在本研究中，一些lncRNAs直接作用于基因，如ENSGALT00000098971(lnc-971)-EGR1、XR_003073121.1 (lnc21.1)-ENSGALG00000005043以及ENSGALT00000098920(lnc-920)既作用于與骨骼肌細胞發(fā)育相關(guān)的ANKRD1又作用于與脂質(zhì)代謝相關(guān)的LSS、DHCR24和SQL，猜測這些lncRNAs可能通過干擾基因表達的作用方式來參與肌肉的生長發(fā)育，具體的調(diào)控方式還需要進一步研究。此外，還有一些lncRNAs不直接作用于基因，而是作用于miRNAs，而這些miRNAs又靶向于mRNAs，如ENSGALT00000095670(lnc-670)和MSTRG.21490.5(lnc-90.5)作用與gga-miR-1635，而gga-miR-1635靶向于與骨骼肌細胞發(fā)育的ASB2、FOS以及ANKRD1，暗示lnc-670和lnc-90.5可能作為ceRNAs與相關(guān)基因競爭結(jié)合于gga-miR-1635，從而調(diào)控雞骨骼肌的生長發(fā)育。從調(diào)控網(wǎng)絡(luò)來看，與骨骼肌細胞生長發(fā)育以及脂質(zhì)代謝相關(guān)基因都受到多個miRNAs的調(diào)控，已有大量研究揭示miRNA與基因的相互作用，如本研究中靶向于FOS的gga-miR-1625-5p，可能與胚胎期肌肉的生長發(fā)育有關(guān)，這與Zhou等[52]在研究不同胚齡扁雞骨骼肌miRNA時發(fā)現(xiàn)gga-miR-1625-5p在扁雞肌肉的生長發(fā)育中發(fā)揮重要的作用相一致。本研究表明，ENSGALT00000098971(lnc-971)、XR_003073121.1(lnc-21.1)、ENSGALT00000098920(lnc-920)、ENSGALT00000095670(lnc-670)以及MSTRG.21490.5(lnc-90.5)可能在調(diào)節(jié)雞肌肉生長發(fā)育以及脂質(zhì)代謝方面發(fā)揮重要的作用，可作為候選lncRNAs進行深入研究。

4 結(jié) 論

本研究通過轉(zhuǎn)錄組測序，獲得藏雞和大恒肉雞胚胎期腿肌中差異表達的106個mRNAs和28個lncRNAs，通過GO和KEGG富集分析篩選得到了與骨骼肌發(fā)育相關(guān)的mRNAs和lncRNAs，最后構(gòu)建了胚胎期骨骼肌發(fā)育相關(guān)的lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)。該結(jié)果為進一步解析不同雞品種骨骼肌生長發(fā)育的差異性提供了理論依據(jù)。