不同雞種胚胎腿肌差異表達mRNA和lncRNA的篩選及其競爭性調控網絡的構建

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

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

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

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

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

1 材料與方法

1.1 試驗動物及樣品采集

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

1.2 文庫構建及轉錄組測序

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

1.3 測序數據質控及序列比對

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

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

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

以大恒肉雞為對照組，藏雞為處理組，使用RSEM軟件(http://deweylab.github.io/RSEM/)[21]通過單端或雙端測序數據來計算基因或者轉錄本的表達量，并以TPM(Transcripts Per Million reads)作為表達衡量指標。使用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為標準，篩選顯著富集的GO條目和KEGG通路。對于篩選得到的差異表達lncRNAs(differentially expressed lncRNAs，DELs)，首先根據lncRNA不同的作用模式進行順式(cis)作用靶基因預測和反式(trans)作用靶基因預測，然后將預測得到的DELs的靶基因進行GO富集分析和KEGG通路分析，最后篩選顯著富集的GO條目和KEGG通路。挑選得到與肌肉生長發(fā)育和脂質代謝相關的GO條目和KEGG通路。

1.5 lncRNA-miRNA-mRNA競爭性調控網絡的構建

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

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

為保證測序結果的準確可靠性，隨機選取基因和lncRNAs各5個，以GAPDH為內參基因，使用SYBR Green法進行qRT-PCR。qRT-PCR驗證過程中所用到的lncRNAs、基因以及引物信息見表1。反應的體系為：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。每個樣本設置3個技術重復，采用2-ΔΔCt法計算基因和lncRNAs的相對表達量。

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

2 結 果

2.1 雞胚胎腿肌轉錄組測序數據統計

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

表2 測序數據質量分析

2.2 雞胚胎腿肌轉錄組測序結果比對分析

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

表3 比對統計結果

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

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

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

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

圖1 藏雞和大恒肉雞腿肌lncRNA的長度分布 (A) 及分類統計 (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為處理組，差異表達分析結果顯示，差異表達的mRNAs有106個，其中上調的有48個，下調的有58個(圖3A)；得到差異表達的lncRNAs有28個，其中上調的有10個，下調的有18個(圖3B)。之后對DELs的靶基因進行預測，得到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富集結果顯示，共有8個GO條目上出現顯著富集(圖4)，其中與骨骼肌生長發(fā)育有關的條目是骨骼肌細胞分化(GO： 0035914)，有較多DEMs富集在該條目。KEGG富集結果顯示，共有9個KEGG通路顯著富集(圖5)，其中與脂質代謝相關的通路是類固醇生物合成信號通路(map00100)。

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

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

2.6 lncRNA-miRNA-mRNA競爭性調控網絡分析

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

2.7 實時熒光定量PCR驗證

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

3 討 論

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

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

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

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

圖中矩形節(jié)點代表lncRNAs；圓形節(jié)點代表miRNAs；三角形節(jié)點代表mRNAs，上方4個基因是與骨骼肌細胞分化相關基因，下方5個基因是與脂質代謝相關基因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調控網絡分析Fig.7 Regulatory network analysis of lncRNA-miRNA-mRNA

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

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

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

4 結 論

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