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

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

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

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

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

1 材料與方法

1.1 試驗(yàn)動(dòng)物及樣品采集

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

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

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

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

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

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

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

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

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

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

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

1.6 實(shí)時(shí)熒光定量PCR(qRT-PCR)驗(yàn)證

為保證測(cè)序結(jié)果的準(zhǔn)確可靠性，隨機(jī)選取基因和lncRNAs各5個(gè)，以GAPDH為內(nèi)參基因，使用SYBR Green法進(jìn)行qRT-PCR。qRT-PCR驗(yàn)證過(guò)程中所用到的lncRNAs、基因以及引物信息見(jiàn)表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。擴(kuò)增程序?yàn)椋?5 ℃ 30 s；95 ℃ 10 s，60 ℃ 30 s，30個(gè)循環(huán)；熔解曲線程序：95 ℃ 15 s，60 ℃ 1 min，95 ℃ 15 s。每個(gè)樣本設(shè)置3個(gè)技術(shù)重復(fù)，采用2-ΔΔCt法計(jì)算基因和lncRNAs的相對(duì)表達(dá)量。

表1 實(shí)時(shí)熒光定量PCR引物信息

2 結(jié) 果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.7 實(shí)時(shí)熒光定量PCR驗(yàn)證

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

3 討 論

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

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

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

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

圖中矩形節(jié)點(diǎn)代表lncRNAs；圓形節(jié)點(diǎn)代表miRNAs；三角形節(jié)點(diǎn)代表mRNAs，上方4個(gè)基因是與骨骼肌細(xì)胞分化相關(guān)基因，下方5個(gè)基因是與脂質(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)錄組測(cè)序中DEMs (A) 和DELs (B) 的qRT-PCR驗(yàn)證Fig.8 Validation of DEMs (A) and DELs (B) in RNA-seq by qRT-PCR

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

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

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

4 結(jié) 論

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