不同鸡种胚胎腿肌差异表达mRNA和lncRNA的筛选及其竞争性调控网络的构建

李 洁，赵瑞鹏，陈楚雯，杨朝武，吴锦波，李志雄*

(1.西南民族大学 青藏高原动物遗传资源保护与利用教育部重点实验室，成都 610041；2.西南民族大学畜牧兽医学院，成都 610041； 3.西南民族大学青藏高原研究院，成都 610041；4.四川省畜牧科学研究院，成都 610066； 5.阿坝藏族羌族自治州畜牧科学技术研究所，阿坝 624402)

藏鸡是我国特有的鸡种，于2020年5月29日被列入《国家畜禽遗传资源品种名录》，主要分布在西藏、青海等海拔较高的地区，是当地非常重要的畜禽品种之一。该品种具有抗病力强且耐粗饲的优点，其肉质鲜美，肌肉中的必需氨基酸和鲜味氨基酸含量丰富[1]。经过长期的自然选择，藏鸡表现出对高海拔环境稳定的遗传适应性[2]。因此，藏鸡是发展高原养殖业不可替代的品种资源。虽然藏鸡肉品质好，但生长速度慢，与普通的肉鸡相比体重和体尺较小，产肉率和屠宰性能不佳。大恒肉鸡是由四川省畜牧科学研究院培育的肉鸡新品种，该品种不仅保持了地方鸡种的优质风味和外观性状，且生长速度快，生产性能和屠宰率高[3-4]。根据2018年的统计显示，商品代大恒肉鸡70日龄公鸡和母鸡的体重分别为2.32和1.79 kg，而成年藏鸡公鸡和母鸡的体重分别仅为1.44和1.15 kg[5]。禽体的肌肉约占体重的30%～40%，肌纤维是构成肌肉的基本单位，肌纤维的数量在雏鸡孵化前就已基本固定。在鸡胚胎期肌肉发育的两个主要阶段中，对肌肉生长发育有决定作用的肌细胞在胚胎期间不可逆的退出细胞周期并不再增殖[6]。除此之外，在胚胎末期，成熟的肌纤维周围会出现参与肌纤维肥大和成熟的肌卫星细胞[7]，因此，比较分析藏鸡和大恒肉鸡这两个品种在胚胎期骨骼肌生长发育上的差异对解决藏鸡生长速度慢，产肉率低等问题具有重要的意义。

影响动物肌肉生长发育的因素有很多，例如：基因、营养、饲养管理、环境、品种以及年龄等。其中基因的影响是最根本、最具有决定性的。近年来，随着分子生物学的快速发展以及生物信息技术的不断进步，基因在调控畜禽经济性状中的功能以及机制不断被挖掘和破解。例如：肌分化因子MyoD1[8]通过激活肌肉特异性基因来调控肌肉的生长发育、生肌调节因子(myogenic regulatory factors，MRFs)家族，如肌细胞生成素MyoG参与调控骨骼肌细胞的增殖、分化以及融合[9]，MyoG、MyoD和MRF4是成肌细胞分化所必需的[10]、肌肉骨骼胚胎核蛋白1(MUSTN1)基因[11]与肌肉的分化和再生有关等，这些都是在肌肉生长发育过程中起重要作用的基因。

lncRNA是长度大于200 nt的非编码RNA[12]。迄今为止，已有大量的研究证明lncRNA可以通过干扰基因的表达、与蛋白质相结合、作为竞争性内源RNA(competitive endogenous RNA，ceRNA)以及编码肽类等多种作用方式直接或间接的参与家畜生命活动的几乎每一个过程[13-16]。在畜禽肌肉的生长发育的过程中，lncRNA常作为ceRNA与肌肉生长发育相关基因竞争结合同一个miRNA从而间接调控肌肉生长发育相关基因的表达，最终影响畜禽肌肉的生长发育，参与调控。也有越来越多特定的lncRNA被发现参与调控畜禽肌肉的生长发育，例如有研究表明lncRNA MEG3参与猪骨骼肌分化的调控，并且作为提高猪肉产量的候选基因[17]，lncRNA-Six1能够通过编码微肽间接的促进鸡骨骼肌细胞的增殖[18]，lncMD1可作为ceRNA与miR-133和miR-135相结合，进而调节与晚期肌细胞增殖分化相关的转录因子表达，从而调控牛肌肉的生长发育[15]，以及lnc-002783能够抑制山羊骨骼肌卫星细胞的增殖和分化[19]。大量的研究证实，lncRNA在很多畜禽经济性状的调控中普遍存在，并发挥着重要的作用。但lncRNA在藏鸡肌肉生长发育中的相关研究却很少。因此，本研究通过转录组测序筛选了胚胎期藏鸡和大恒肉鸡腿肌中差异表达的mRNAs和lncRNAs，挑选与肌肉生长发育相关的mRNAs和lncRNAs构建竞争性调控网络，为后续深入分析两鸡种在肌肉生长发育方面呈现差异的分子机制奠定基础。

1 材料与方法

1.1 试验动物及样品采集

本试验所用藏鸡和大恒肉鸡的种蛋分别由阿坝州九顶原生态畜禽养殖有限责任公司和四川大恒家禽育种有限公司(国家肉鸡核心育种场)提供。将所有种蛋按照常规程序在相同环境条件下进行孵化，胚胎期第18天，在两个品种的种蛋中分别挑选3个发育正常的雄性鸡胚，快速采集腿肌组织并装入无RNA酶的冻存管中，藏鸡的腿肌组织标记为TC_E18，大恒肉鸡标记为DH_E18，之后迅速置于液氮中保存备用。

1.2 文库构建及转录组测序

按照TRIzol试剂盒(Invitrogen)说明书提取藏鸡和大恒肉鸡腿肌组织的总RNA。之后通过1%琼脂糖凝胶电泳检测RNA的质量，利用NanoDrop ND-2000(Agilent)检测RNA的浓度和纯度。将检测合格后的总RNA使用PrimeScriptTMRT reagent Kit with gDNA Eraser(Perfect Real Time)反转录合成cDNA链，置于-80 ℃保存备用。两个品种采用去除核糖体RNA的方法共构建6个特异性文库并进行转录组测序(转录组文库的构建和测序由上海美吉生物医药科技有限公司完成)。

1.3 测序数据质控及序列比对

通过Illumina二代测序技术进行测序后产生大量以fastq格式储存起来的原始数据(raw data)，对每个样本的原始测序数据进行相关质量评估。通过使用fastp软件对碱基质量分布、碱基错误率分布以及碱基含量分布进行统计以保证之后生物信息分析的准确可靠性。由于原始数据中会包含测序接头序列、低质量读段、不确定碱基信息(N率)较高的序列以及长度过短序列等，这些都会严重影响后续分析的准确性。因此统计结果符合要求之后，使用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/)软件进行转录本组装获得每个样本的转录本信息以及发现新的转录本或新基因。

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通路。挑选得到与肌肉生长发育和脂质代谢相关的GO条目和KEGG通路。

1.5 lncRNA-miRNA-mRNA竞争性调控网络的构建

收集富集到与肌肉生长发育和脂质代谢相关的GO条目和KEGG通路上的基因。通过lncRNA靶基因预测，发现有些lncRNA可以同时作用于两个基因；之后通过miRNA靶基因预测发现，有miRNA靶向于与肌肉生长发育和脂质代谢相关的lncRNA和基因上；还有些lncRNA虽然没有直接作用于与肌肉生长发育和脂质代谢相关的基因上，但通过作用于miRNA可间接作用于相关基因。因此，通过lncRNA和miRNA靶基因预测，收集与肌肉生长发育和脂质代谢相关的lncRNAs、miRNAs和基因整理其联系。利用Cytoscape V-3.6.0软件进行lncRNA-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个循环；熔解曲线程序：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。表明测序的质量高，数据可靠准确，可用于后续分析。

表2 测序数据质量分析

2.2 鸡胚胎腿肌转录组测序结果比对分析

将质控后得到的高质量读段(clean reads)与鸡参考基因组比对，结果显示DH_E18和TC_E18组比对到鸡参考基因组上的clean reads比对率为90.05%～94.35%，在参考序列上有多个比对位置的clean reads比率为12.69%～17.72%，比对效果较好，满足后续分析要求(表3)。

表3 比对统计结果

2.3 鸡胚胎腿肌中lncRNA的长度分布及分类统计

统计文库中所有lncRNA长度分布，大多数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的表达量进行分析，发现共有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)，其中与骨骼肌生长发育有关的条目是骨骼肌细胞分化(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个节点，产生了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 讨 论

藏鸡作为我国特有的地方品种，是西藏青海等高海拔地区发展高原养殖业的优良鸡种。关于藏鸡的研究大多集中在产蛋性能以及对高海拔的适应性上，有关肌肉生长发育以及肉品质方面的研究还较少。本研究对胚胎第18天的藏鸡和大恒肉鸡腿肌进行全转录组测序，筛选DEMs和DELs，然后进一步对其进行功能分析，并建立竞争性调控网络进行联合分析，为进一步探索鸡肌肉生长发育的分子机制提供一定的理论基础。

图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条通路与免疫系统相关。补体活化系统是免疫应答中固有免疫的重要组成部分，研究表明，补体活化途径与骨骼的生长发育以及骨重塑之间有着密切的联系[23]。补体C3是补体活化单位之一，Matsuoka等[24]研究表明，在核因子κB受体活化因子配体(RANKL)注射的高骨转化状态下，补体C3的表达量升高，在每天服用抑制补体C3激活的药物状态下，导致骨形成减弱，骨丢失加剧；动物的产肉量和体型大小与骨骼的生长发育密切相关，成骨细胞和破骨细胞对于骨骼的形成、发育、生长以及修复具有重要的作用，而补体活化途径对成骨细胞和破骨细胞的生长和分化具有重要的调节作用[22]，进而会影响骨骼的形成、发育和生长。提示可以从补体活化的角度进行深入研究，可能会在一定程度上提高藏鸡的生长速度，增大其体格。近年来，越来越多的研究结果都表明，免疫应答对肌肉的生长发育有一定的影响，Li等[25]对泸宁鸡和白羽肉鸡腓肠肌的差异表达基因进行GO富集分析发现主要富集的条目为免疫应答，表明免疫应答对肌肉的生长发育具有一定的影响；Dou等[26]以大体型的武定鸡和小体型的大围山鸡为样本进行比较研究，发现在胸肌组织中，与大围山鸡相比，武定鸡的免疫功能相关基因表达增加；袁茂等[27]对不同发育阶段的藏鸡腿肌的差异表达基因进行KEGG富集分析发现，在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)。

图中矩形节点代表lncRNAs；圆形节点代表miRNAs；三角形节点代表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]，对于肌肉的生长发育具有重要的调控作用。研究表明，MYF6基因参与肌纤维的成熟，并且维持肌纤维的分化状态[29]，与MyoD家族中的其它成员协同调控，在肌肉生长发育的各阶段发挥关键作用[30-32]。汤展毅等[33]通过克隆MYF6基因并构建载体，运用脂质体转染法转染鲁西黄牛成纤维细胞，检测MYF6对成肌纤维细胞的影响，发现与肌肉生长发育相关基因的表达量显著升高，且转染后的成纤维细胞发育为肌管，表明了MYF6具有促进肌细胞分化的功能；Li等[34]通过研究MYF6基因在家鸭胚胎和孵化后胸肌和腿肌中的表达模式证明，MYF6与肌肉生长发育有关，是骨骼肌生长发育的一个组成部分。FOS和EGR1基因属于立早基因(IEGs)，在细胞的增殖分化中具有重要的作用[35-36]。FOS参与肌肉的生成[37]，李志雄等[38]研究发现，在注射了使肉鸡腓肠肌纤维直径显著减小的miR-499-5p的处理组中FOS的表达量显著高于对照组，提示FOS基因会抑制肌肉纤维的生长发育；但也有研究表明，FOS下调，通过调节TNF-α的表达进而激活NF-κB以及IGF-1信号通路的Akt酶从而抑制肌管的形成以及肌源性分化[39]，通过抑制酪氨酸磷酸化从而影响肌生成素的表达[40]，提示FOS可能促进肌肉的生长发育。表明肌肉的生长发育机制是一个极其复杂的过程。EGR1主要在肌肉中表达，已有研究表明，EGR1调控肌肉细胞的生长和分化[41]。Yang等[42]通过抑制靶向于EGR1的lncRNA-COX2，使成纤维细胞中的EGR1显著增加，纤维生成显著促进，并且成纤维细胞显示出比对照组更多的迁移细胞，而在EGR1被抑制后，这些迁移被阻断，说明EGR1基因具有促进成纤维细胞活化以及肌纤维生成的作用；崔亚凤[43]利用CRISPRi技术激活和抑制EGR1，发现EGR1的表达量与小鼠成肌细胞C2C12的分化呈显著的正相关，证明了EGR1能够促进肌肉的分化。ASB2基因是ASB家族成员之一，有研究发现该基因可能参与骨骼肌的分化、参与肌肉再生以及肌肉损伤修复等过程[44]。康亭亭[45]在研究猪骨骼肌卫星细胞时，发现ANKRD1在其分化过程中的表达量较高且在约克夏猪和监利猪中的表达量差异较大，推测ANKRD1可能参与调控肌卫星细胞的分化过程；除此之外，也有研究表明ANKRD1可能也参与调控成肌纤维细胞的生长[46]。这些基因都能够促进肌肉的生长发育，但其表达量均显著下调，由此提示相较于大恒肉鸡来说，藏鸡的生长速度较慢、产肉率低可能与这些基因的表达水平较低有关。可以对这些基因进行深入探索，揭示其具体的调控机制，为提高藏鸡生长速度，改良其产肉性能的研究提供了思路。

lncRNA靶基因的KEGG显著富集的结果中最可靠的是脂质代谢相关通路，其中类固醇生物合成是最显著富集的KEGG通路，该通路同样出现在DEMs的KEGG显著富集通路中。富集在该通路的基因分别是DHCR24、SQLE以及LSS，且ENSGALT00000098920(lnc-920)共同靶向于这3个基因。这些基因都与胆固醇的合成有关，分别是催化胆固醇合成的起始、最终形成以及胆固醇前体的重要因子[47-49]。Mirza等[50]敲除小鼠肌肉组织中的DHCR24后，敲除小鼠的体积较同窝中未敲除小鼠的体积小且在几小时后死亡，表示DHCR24可能对肌肉的生长发育也有一定的影响。查询差异表达量，发现这3个基因均显著上调，而lnc-920显著下调。猜测lnc-920的表达量可能与脂质相关基因的表达量呈负相关，具体调控机制还需要进一步研究。因此，相较于大恒肉鸡，藏鸡的肉品质较好，其脂质代谢相关基因的表达量高可能是原因之一。同时值得注意的是，已有研究表明上述几个基因可能与鸡的脂肪沉积有关[47]，加之人们对畜禽产品需求的改变，提示藏鸡作为高原养殖业的优良鸡种还有较大的遗传改良潜质。脂肪酸的含量及类型会影响肉的风味，富集在脂肪酸生物合成通路上的基因对应的lncRNA中MSTRG.21490.5(lnc-90.5)的表达量在藏鸡和大恒肉鸡腿肌中分别是其它DELs表达量的约7和20倍，推测该lncRNAs可能对肉质有一定的影响，且该lnc-90.5是新发现的lncRNA，具有极大潜在的研究价值。

本研究基于DEMs和DELs靶基因富集分析得到的富集于骨骼肌生长分化以及脂质代谢等条目上的基因，通过对相关lncRNAs以及miRNAs进行预测，构建lncRNA-miRNA-mRNA竞争性调控网络。近年来，许多研究发现lncRNA可以通过多种作用方式参与生物体的各种生物学过程。其中，lncRNA作为ceRNA与miRNA相结合进而发挥作用的报道已经屡见不鲜。例如：Sun等[15]通过转录组测序研究秦川牛骨骼肌时发现lncMD1作为ceRNA与骨骼肌分化相关基因IGF2相竞争结合于miR-125b，最后提高了IGF2的表达，进而促进了骨骼肌的分化；巨晓军[51]以清远麻鸡的胸大肌和缝匠肌为样本，构建lncRNA-miRNA-mRNA互作网络，发现lncXR_003074785.1-miR-193-3p-PPARGC1A互作通路，其中PPARGC1A被证明调节鸡骨骼肌纤维的形成以及肌细胞的分化，且miR-193-3p与PPARGC1A的表达呈负相关，暗示lncXR_003074785.1(lnc-85.1)可能作为ceRNA与miR-193-3p相结合从而影响肌肉的生长发育。在本研究中，一些lncRNAs直接作用于基因，如ENSGALT00000098971(lnc-971)-EGR1、XR_003073121.1 (lnc21.1)-ENSGALG00000005043以及ENSGALT00000098920(lnc-920)既作用于与骨骼肌细胞发育相关的ANKRD1又作用于与脂质代谢相关的LSS、DHCR24和SQL，猜测这些lncRNAs可能通过干扰基因表达的作用方式来参与肌肉的生长发育，具体的调控方式还需要进一步研究。此外，还有一些lncRNAs不直接作用于基因，而是作用于miRNAs，而这些miRNAs又靶向于mRNAs，如ENSGALT00000095670(lnc-670)和MSTRG.21490.5(lnc-90.5)作用与gga-miR-1635，而gga-miR-1635靶向于与骨骼肌细胞发育的ASB2、FOS以及ANKRD1，暗示lnc-670和lnc-90.5可能作为ceRNAs与相关基因竞争结合于gga-miR-1635，从而调控鸡骨骼肌的生长发育。从调控网络来看，与骨骼肌细胞生长发育以及脂质代谢相关基因都受到多个miRNAs的调控，已有大量研究揭示miRNA与基因的相互作用，如本研究中靶向于FOS的gga-miR-1625-5p，可能与胚胎期肌肉的生长发育有关，这与Zhou等[52]在研究不同胚龄扁鸡骨骼肌miRNA时发现gga-miR-1625-5p在扁鸡肌肉的生长发育中发挥重要的作用相一致。本研究表明，ENSGALT00000098971(lnc-971)、XR_003073121.1(lnc-21.1)、ENSGALT00000098920(lnc-920)、ENSGALT00000095670(lnc-670)以及MSTRG.21490.5(lnc-90.5)可能在调节鸡肌肉生长发育以及脂质代谢方面发挥重要的作用，可作为候选lncRNAs进行深入研究。

4 结 论

本研究通过转录组测序，获得藏鸡和大恒肉鸡胚胎期腿肌中差异表达的106个mRNAs和28个lncRNAs，通过GO和KEGG富集分析筛选得到了与骨骼肌发育相关的mRNAs和lncRNAs，最后构建了胚胎期骨骼肌发育相关的lncRNA-miRNA-mRNA调控网络。该结果为进一步解析不同鸡品种骨骼肌生长发育的差异性提供了理论依据。