鸭胚胎发育中后期胸肌发育阻滞的RNA-seq分析

刘宏祥，徐文娟，朱春红，陶志云，宋卫涛，章双杰，李慧芳



鸭胚胎发育中后期胸肌发育阻滞的RNA-seq分析

刘宏祥，徐文娟，朱春红，陶志云，宋卫涛，章双杰，李慧芳

（江苏省家禽科学研究所，江苏扬州 225125）

【目的】选择两个中国地方品种高邮鸭和金定鸭，对鸭胚胎发育中后期胸肌进行转录组分析研究，旨在探明胸肌发育阻滞的分子变化机制，为鸭骨骼肌调控机理研究打下基础。【方法】在21胚龄和27胚龄两个时间点，分别解剖高邮鸭、金定鸭各3只，采集胸大肌，提取总RNA构建文库，利用Illumina的HiseqTM2000进行高通量测序，并利用生物信息学方法进行差异表达基因挖掘、基因功能注释等分析，探讨21胚龄和27胚龄两个时间点之间胸肌发育阻滞的分子机制。【结果】高邮鸭、金定鸭21胚龄和27胚龄胸大肌组织RNA-seq质量Q20均在94%以上，Q30均在89%以上，测序得到的结果可靠，可用于后续分析。RNA水平相关性检查和基因mRNA表达量聚类图结果都表明，21胚龄（27胚龄）高邮鸭和金定鸭之间表达模式的相关性高于高邮鸭（金定鸭）21胚龄和27胚龄之间的表达模式。不同品种内时间点之间的差异基因数量（高邮鸭6 128个，金定鸭6 452个）远多于同一时间点不同品种间的差异基因数量（21胚龄522个，27胚龄299个）。qRT-PCR验证试验结果与RNA-seq分析结果相关性较强。通过GO和KEGG富集分析发现，高邮鸭、金定鸭胸肌在21胚龄到27胚龄阶段，能量代谢相关基因（主要为辅酶Q相关基因、ATP酶合成相关基因和细胞色素C相关基因）均显著上调，DNA复制和细胞周期相关基因（主要为微型染色体维持蛋白（MCM）相关基因、复制因子C（RFC）相关基因）均显著下调。相关基因表达的变化可能与此阶段成肌细胞增殖速度减慢，逐渐退出细胞周期开始准备下一阶段融合成多核肌管并形成肌纤维有关。对肌肉生长发育相关的关键基因分析发现，促进肌肉生长的1和诱导成肌细胞末端分化的显著下调，促进肌纤维分化融合的1基因、诱导肌祖细胞向成肌细胞转化的1基因显著上调。【结论】鸭胚胎中后期胸肌发育过程中大量基因差异表达。其中能量代谢相关基因的上调和DNA复制、细胞周期相关基因的下调，以及肌肉发育相关基因1显著上调，1、等显著下调，可能与鸭胚胎中后期胸肌发育阻滞现象密切相关。

鸭；胸肌发育；转录组测序；差异表达基因

0 引言

【研究意义】动物肌肉的生长发育大体可以分为两个阶段：出雏（生）前的胚胎期阶段和出雏（生）后阶段。动物肌纤维数量在出雏（生）前就已经固定[1-3]，后期肌肉量的增加主要由于肌纤维的肥大，因此禽类上市日龄时的产肉量取决于胚胎中后期肌纤维的数量。胸肌是禽类肌肉的重要组成部分，因此研究禽类胚胎期胸肌发育的调控机制，挖掘影响肌纤维发育的关键基因，对育种实践中禽类产肉量的选育具有重要指导意义。【前人研究进展】动物在胚胎发育期间，成肌细胞增殖、分化成为多核的肌管，最终形成成熟的肌纤维，此过程涉及肌纤维的增生和肥大[2]。哺乳动物胚胎期胸肌重量持续增长[4]，而禽类胚胎期胸肌发育过程较为复杂，表现为胚胎前期发育迅速而中后期发育迟缓甚至萎缩。MOORE等[5]报道火鸡胸肌肌纤维横截面积在胚胎后期减小。CHEN等[6-7]对鸭的研究发现，鸭胸肌重量和肌纤维横截面积从22胚龄到出雏期呈减小趋势。本课题组前期研究发现高邮鸭和金定鸭胚胎期21胚龄后胸肌重量没有增加[8]。这可能与此时胸肌肌纤维主要表现为融合，而非肥大有关[7]，但具体分子机制还未完全了解。近年来，随着鸭基因组测序的完成[9]（基因组草图见http://www.ensembl.org/Anas_platyrhynchos/Info/ Index），以及高通量测序、分析技术的发展，RNA-seq分析鸭不同发育时期基因mRNA表达水平变化情况及其调控机制成为可能。【本研究切入点】鸭胚胎中后期胸肌发育阻滞现象已有多方报道，但具体分子机制还未完全明晰。【拟解决的关键问题】高邮鸭和金定鸭分别为中国地方体型差异较大的品种，其肌肉生长具有明显的表型差异，这为研究肌肉生长发育的分子机制提供了良好的素材。本研究拟通过RNA-seq技术方法，探讨胚胎发育中后期高邮鸭、金定鸭的胸肌转录组表达差异情况，为深入了解此阶段胸肌生长发育的分子机制奠定基础。

1 材料与方法

1.1 试验动物及样品采集

本试验研究对象为中国地方品种鸭——蛋肉兼用型品种高邮鸭和蛋用型品种金定鸭，于2014年饲养于江苏省高邮鸭集团。选择正常的高邮鸭和金定鸭种蛋（重量接近）各10枚于微电脑全自动孵化器（山东德州，科裕孵化机）中孵化，21胚龄和27胚龄时随机挑选高邮鸭、金定鸭胚蛋各3枚。解剖胚胎，采集右侧胸大肌样品，迅速置于液氮中冻存。

1.2 试验方法

1.2.1 总RNA提取与质量鉴定 使用Trizol试剂（Invitrogen）从约30 mg胸大肌样品中提取总RNA，使用1%甲醛凝胶对RNA样品进行电泳，鉴定RNA完整性和降解情况；使用Qubit RNA试剂盒（Invitrogen）测定RNA浓度；使用Nano Photometer分光光度计（德国Implen）检测RNA纯度。

1.2.2 测序文库的构建 RNA检测合格后，将mRNA打断成短片段并反转录成cDNA，选择合适大小的片段进行PCR富集得到最终的cDNA文库。使用qRT-PCR方法对文库的有效浓度进行准确定量（文库有效浓度＞2 nM），以保证文库质量。

1.2.3 Illumina测序 库检合格后，交由北京诺禾致源公司使用HiseqTM2000测序仪（Illumina）进行转录组测序。

1.2.4 qRT-PCR验证 根据转录组测序结果，挑选16个差异表达基因进行qRT-PCR验证。使用Primer 3在线工具设计引物（http://fokker.wi.mit.edu/primer3/ input.htm）（所选基因和设计的引物见表1）。

表1 用于qRT-PCR验证的基因及其引物

*看家基因Housekeeping genes

1.3 测序数据的处理分析

1.3.1 测序数据的处理 为保证信息分析质量，对原始序列（Raw reads）进行过滤去杂（去除测序接头序列、重复冗余序列、低质量序列）得到高质量的序列数据（clean reads）。采用TopHat2方法[10]将clean reads与鸭参考基因组进行比对。

1.3.2 差异基因的筛选 本试验每个胚龄点都有3个生物学重复。通过DESeq软件包，对前步得到的clean reads进行标准化，并评估read count的分布，进行BH检验（Benjamini和Hochberg于1995年提出的控制假阳性结果错误率的多重比较方法[11]），最终得到不同组别之间的差异基因比较结果。

1.3.3 mRNA差异表达基因GO分类和KEGG富集分析 使用R语言中的goseq包[12]将mRNA差异表达基因序列与GO（gene ontology）数据库进行比对分析，获得GO功能注释；将mRNA差异表达基因与KEGG（kyoto encyclopedia of genes and genomes）数据库进行BLASTX比对获得mRNA差异表达基因相对应的Pathway注释信息。

2 结果

2.1 RNA-seq数据总体评价

本研究对高邮鸭、金定鸭21胚龄和27胚龄胸大肌组织转录组测序，高邮鸭和金定鸭两个时间点均获得了43+ M条的reads，总碱基数5+ Gb，Q20接近95%，Q30接近90%（表2），说明RNA-seq测序结果可靠，可用于后续分析。高邮鸭、金定鸭21胚龄和27胚龄样品比对到参考基因组上的reads均达到总reads数的60%左右（表3），比对率均较高。

表2 测序数据评估统计

Q20、Q30分别指质量值大于等于20或30的碱基所占百分比。GC为测序结果中G和C碱基所占总碱基含量的百分比

Q20 and Q30 means the percentage of quality value grater or equal to 20 or 30 base, respectively; GC means the percentage of G and C bases in total bases sequenced

表3 Reads与参考基因组比对情况

2.2 RNA-seq水平相关性检查

使用皮尔逊相关系数对样品进行相关性分析。同一时间点、品种内的3个生物学重复的相关系数平方（2）均大于0.95，同一时间点不同品种之间的相似性高于同一品种不同时间点之间的相似性（图1）。

2.3 基因mRNA表达聚类图

根据高邮鸭、金定鸭21胚龄、27胚龄各基因RPKM进行层次聚类分析。同一时间点不同品种之间的表达模式相近，而同一品种不同时间点的表达模式差异较大（图2）。

2.4 mRNA差异表达基因的筛选

表4显示了21胚龄和27胚龄样品的RPKM值。21胚龄高邮鸭、金定鸭胸大肌组织高表达基因较多，RPKM值大于3的基因均占到60%左右。

表4 不同胚龄胸大肌RPKM值

图1 样品间RNA-seq水平相关性检查

图2 基因mRNA差异表达聚类图

高邮鸭和金定鸭不同时间点之间的mRNA差异表达基因均较多。高邮鸭21胚龄和27胚龄之间、金定鸭21胚龄和27胚龄之间的mRNA显著差异表达基因数量均超过6 000个。21胚龄金定鸭和高邮鸭mRNA显著差异表达基因、27胚龄金定鸭和高邮鸭mRNA显著差异表达基因分别为522个和299个（图3），明显少于品种内不同时间点之间的mRNA差异表达基因数量。

纵坐标1.3以上表示padj<0.05 Numbers higher than 1.3 in Y-axis indicate padj＜0.05

2.5 qRT-PCR验证

从差异表达分析结果中挑选出16个基因进行qRT-PCR验证。结果显示，所选择的16个基因qRT-PCR结果（∆∆Ct值）与RNA-seq结果（RPKM比值）具有较强的相关性（> 0.70）（表5），表明RNA-seq结果可靠。

2.6 骨骼肌相关基因差异表达情况

挑选出与骨骼肌生长发育相关的关键基因1、1、1、和，发现在21胚龄和27胚龄两个时间点，高邮鸭和金定鸭之间1、1、、1 mRNA表达水平均没有显著变化。两个品种21胚龄到27胚龄1、均显著下调，1、1均显著上调，均没有显著变化（表6）。

2.7 mRNA差异表达基因GO富集分析

分别对高邮鸭、金定鸭21胚龄和27胚龄胸大肌组织的上调差异基因、下调差异基因进行GO（Gene Ontology）分析，高邮鸭和金定鸭分别富集到24个（上调和下调分别2个和22个）和26个显著GO条目（上调和下调分别1个和25个（图4）。

B1：细胞成分运动（cellular component movement），B2：基于微管的过程（microtubule-based process），B3：基于微管的运动（microtubule-based movement），B4：细胞周期（cell cycle），B5：分解代谢过程的正向调控（positive regulation of catabolic process），B6：自噬调控（regulation of autophagy），B7：自噬正调控（positive regulation of autophagy），B8：细胞膜组织（cellular membrane organization），B9：（evasion or tolerance of host defenses by virus），B10：细胞分解过程的正向调控（positive regulation of cellular catabolic process），B11：宿主防御逃逸（avoidance of host defenses），B12：宿主防御回避或耐受（evasion or tolerance of host defenses），B13：其它共生互作生物的防御逃逸avoidance of defenses of other organism involved in symbiotic interaction），B14：其他共生互作生物的防御回避或耐受（evasion or tolerance of defenses of other organism involved in symbiotic interaction），B15：其他共生互作生物的防御响应（response to defenses of other organism involved in symbiotic interaction），B16：对宿主防御的响应（response to host defenses），B17：对宿主的响应（response to host）

C1：线粒体（mitochondrion），C2：核糖体（ribosome），C3：核糖核蛋白复合物（ribonucleoprotein complex），C4：细胞外基质（extracellular matrix），C5：大分子复合物（macromolecular complex），C6：非膜结合细胞器（non-membrane-bounded organelle），C7：胞内非膜结合细胞器（intracellular non-membrane-bounded organelle）

M1：核糖体结构成分（structural constituent of ribosome），M2：肌动活性（motor activity），M3：微管肌动活性（microtubule motor activity），M4：结构分子活性（structural molecule activity），M5：微管结合（microtubule binding），M6：细胞支架蛋白结合（cytoskeletal protein binding），M7：DNA依赖的ATP酶活性（DNA-dependent ATPase activity），M8：微管蛋白结合（tubulin binding），M9：焦磷酸酶活性（pyrophosphatase activity），M10：氧化还原酶活性（oxidoreductase activity），M11：作用于酸酐的水解酶活性（hydrolase activity, acting on acid anhydrides），M12：作用于含磷酸酐的水解酶活性（hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides），M13：核苷三磷酸酶活性（nucleoside-triphosphatase activity）

图4 高邮鸭及金定鸭21胚龄和27胚龄mRNA差异表达基因GO功能分类图

Fig. 4 GO function classification of mRNA differentially expressed genes of GY and JD between 21 ed and 27 ed

两个品种富集到13个共有GO条目（表7），其中，C1（GO:0005739，mitochondrion）为上调基因所富集，与线粒体能量代谢有关；其余均由下调基因所富集，主要与细胞周期、DNA复制有关。每个共有GO条目下的共有基因数量见表7。

2.8 mRNA差异表达基因KEGG富集分析

本试验分别对高邮鸭、金定鸭21胚龄和27胚龄胸大肌组织的上调差异基因、下调差异基因进行KEGG分析，高邮鸭和金定鸭分别富集到8个和7个显著KEGG通路（图5）。高邮鸭上调基因富集到4个显著KEGG通路，下调基因富集到4个显著KEGG通路；金定鸭上调基因富集到4个显著KEGG通路，下调基因富集到3个显著KEGG通路。

表5 qRT-PCR实验得到的基因表达∆∆Ct值与RPKM比值之间的相关

表6 骨骼肌相关基因mRNA表达水平变化

表7 高邮鸭、金定鸭21胚龄和27胚龄之间显著富集的共有GO分类

Term2：柠檬酸循环（三羧酸循环）（Citrate cycle (TCA cycle)），Term6：氧化磷酸化（Oxidative phosphorylation），Term7：缬氨酸、亮氨酸和异亮氨酸降解（Valine, leucine and isoleucine degradation），Term9：碳新陈代谢（Carbon metabolism），Term11：核糖体（Ribosome），Term12：DNA复制（DNA replication），Term14：PPAR信号途径（PPAR signaling pathway），Term15：细胞周期（Cell cycle），Term17：ECM受体互作（ECM-receptor interaction），Term18：间隙连接（Gap junction）

两个品种富集到5个共有KEGG通路（表8），其中，Term2（apla00020，Citrate cycle (TCA cycle)）、Term6（apla00190，Oxidative phosphorylation）、Term9（apla01200，Carbon metabolism）由上调基因所富集，主要与能量代谢有关；Term12（apla03030，DNA replication）、Term15（apla04110，Cell cycle）由下调基因所富集，主要与DNA复制和细胞周期有关。每个共有KEGG通路下的共有基因数量见表8。

表8 高邮鸭、金定鸭21胚龄和27胚龄之间显著富集的共有KEGG通路

2.9 GO、KEGG联合分析

对共有GO条目的基因和共有KEGG通路的基因取交集，筛选出16个上调基因和15个下调基因（表9）。上调基因大多为辅酶Q相关基因、ATP酶合成相关基因、细胞色素C相关基因，下调基因大多为微型染色体维持蛋白（MCM）相关基因、复制因子C（RFC）相关基因。

3 讨论

肌肉发育是一个持续不断累积的过程，本课题组前期对高邮鸭、金定鸭的研究[8,13]以及其他课题组的多个试验均报道鸡[14]、鸭[6-7]胚胎发育中后期，胸肌生长迟缓，甚至出现萎缩现象。火鸡胚胎后期胸肌横截面积减小、肌肉萎缩的同时，肌肉卫星细胞的有丝分裂活性也有降低[5]。为了解这一现象背后的分子机制，本试验以高邮鸭、金定鸭21胚龄、27胚龄胚胎胸肌为试验素材，通过RNA-seq比较品种之间以及21胚龄、27胚龄两个时间点之间的mRNA差异表达基因。挑选部分显著差异表达基因进行qRT-PCR验证，结果表明RNA-seq结果与荧光定量PCR结果相似，高通量测序具有较高的可信度。

RNA-seq水平相关性检查发现，同一时间点、品种内的3个生物学重复之间的表达模式相似度较高，能够保证后续差异基因的可靠分析。同一时间点不同品种之间的相关性高于同一品种不同时间点之间的相关性。高邮鸭21胚龄、27胚龄的差异基因和金定鸭21胚龄、27胚龄的差异基因数量远多于21胚龄金定鸭、高邮鸭之间的差异基因数量和27胚龄金定鸭、高邮鸭之间的差异基因数量。对高邮鸭、金定鸭两个时间点的基因按照其RPKM值进行层次聚类分析，21胚龄高邮鸭、金定鸭聚为一类，27胚龄高邮鸭、金定鸭聚为一类。以上结果均表明，在21胚龄到27胚龄的发育过程中，高邮鸭与金定鸭的胸肌表达模式具有相似性，在胸肌发育过程中时间因素的影响要大于品种之间遗传因素的影响。

在GO和KEGG富集分析中，上调基因富集到的GO条目、KEGG通路主要与线粒体能量代谢有关；下调基因富集到的GO条目、KEGG通路主要与细胞周期、DNA复制有关。筛选出GO条目和KEGG通路中的共有基因，发现与能量代谢有关的上调基因主要为辅酶Q相关基因、ATP酶合成相关基因以及细胞色素c相关基因，与DNA复制和细胞周期有关的下调基因主要为微型染色体维持蛋白（MCM）相关基因、复制因子C（RFC）相关基因。辅酶Q对体内呼吸链中的质子移位及电子传递起重要作用，是细胞代谢和细胞呼吸的激活剂[15-16]。细胞色素c也是呼吸链中可流动的递氢体或递电子体。上调基因中，辅酶Q相关基因（NDUF*等）分别参与了NADH还原酶1α和1β等复合物的编码合成；ATP合成酶相关基因（ATP5F1、ATP5H和ATP5J）参与了ATP合成酶复合物的编码合成。21胚龄到27胚龄阶段，鸭胸肌组织中这些基因的上调，增强了能量的转化水平和能量代谢频率。真核细胞DNA的复制和有丝分裂需要微型染色体维持蛋白（MCM，Minichromosome maintenance protein complex）[17-19]和复制因子C[20-22]的参与。21胚龄到27胚龄鸭胸肌DNA复制和细胞周期相关基因下调，说明在此阶段胸肌细胞的增殖速度开始减缓。细胞周期阻滞是后续细胞融合和分化的必要条件，细胞的融合、分化也将使组织器官获取相应的功能和合适的尺寸[23]。

表9 GO与KEGG联合分析后的共有基因

动物机体肌肉的形成可分为多个步骤。在胚胎期，成肌细胞增殖、分化形成多核的肌管，肌管再分化形成肌纤维，最后组装成肌肉组织[2,24-25]。肌肉发育过程需要有许多调控因子的参与，其中1、1、1、和等起着重要作用，因此本试验筛选这5个基因进行详细分析。1是由肝脏主要分泌的调控生长的一种生长因子，近来的研究表明，1也可在肌肉组织局部表达，以旁分泌/自分泌的形式，通过IGF1-Akt/PKB途径调节骨骼肌生长[26]，还可通过IGF1-CaN-NFATc3途径促进成肌细胞分化和肌纤维类型的转换[27]。MUSTN1基因在肌纤维分化、融合过程中具有关键作用，并调节下游靶基因1和[28]。北京鸭上的研究发现，1 mRNA表达水平与鸭胸腿肌的相对生长率相关[29]。1和为MRFs（生肌调控因子）的两个重要成员。在动物胚胎期1可以诱导肌祖细胞向生肌细胞系转变[30,31]，体外试验也表明，1可将其他类型细胞转变为成肌细胞[32-34]。在成肌细胞脱离细胞周期并从增殖过程转变到分化过程中起着重要作用，但不影响成肌细胞的增殖[35-37]。是一种重要的负调控骨骼肌生长的关键基因。在功能缺失的牛[38-40]、鼠[41-43]，其肌肉量显著增加。本研究发现，在21胚龄到27胚龄胸肌发育阻滞的同时，1显著上调（高邮鸭和金定鸭基因表达水平分别提高6.5倍和4.7倍），表明1在鸭胸肌发育中起着重要作用。1与1的mRNA表达模式相似，21胚龄到27胚龄也表现显著上调，而1和均显著下调，均没有显著变化。本试验前期研究中对21胚龄、27胚龄高邮鸭、金定鸭胸肌肌肉生长相关基因进行了qRT-PCR定量分析，发现两品种21胚龄到27胚龄期间1均显著下调，均略微上调，表达在21胚龄和27胚龄两个时间点均没有品种差异[13]。在对北京鸭的研究中，Gu等发现21胚龄到27胚龄阶段，表达水平持续下降，而表达水平持续上升[7]。除了基因，1、和1结果都与本试验RNA-seq结果相似。1的上调将促使更多的肌祖细胞向成肌细胞转化，这为出雏后肌纤维的快速生长提供了更多的材料来源。细胞周期与DNA复制相关基因mRNA表达水平显著下降，这可能导致肌肉卫星细胞的有丝分裂活性降低；另外1、的下调，最终使得鸭胚在接近出雏时胸肌发育阻滞。

4 结论

本试验利用高通量测序技术对高邮鸭、金定鸭胚胎中后期胸肌组织的转录组进行了测序分析，结果表明在此阶段胸肌组织中与能量代谢相关的基因显著上调，与DNA复制、细胞周期相关的基因显著下调，一些肌肉发育相关的关键基因，如1、显著下调，最终导致胸肌肌肉细胞增殖减缓，胸肌发育阻滞。该发现将为下一步深入探索鸭胸肌发育阻滞机制奠定基础。

[1] SMITH J H. Relation of body size to muscle cell size and number in the chicken., 1963, 42(2): 283-290.

[2] PICARD B, LEFAUCHEUR L, BERRI C, DUCLOS M J. Muscle fibre ontogenesis in farm animal species., 2002, 42(5): 415-431.

[3] REHFELDT C, STICKLAND N C, FIEDLER I, WEGNER J. Environmental and genetic factors as sources of variation in skeletal muscle fibre number., 1999, 9(5): 235-254.

[4] SWATLAND H J. Muscle growth in the fetal and neonatal pig., 1973, 37(2): 536-545.

[5] MOORE D T, FERKET P R, MOZDZIAK P E. Muscle development in the late embryonic and early post-hatch poult., 2005, 4(3): 138-142.

[6] CHEN W, TANGARA M, XU J, PENG J. Developmental transition of pectoralis muscle from atrophy in late-term duck embryos to hypertrophy in neonates., 2012, 97(7): 861-872.

[7] GU L H, XU T S, HUANG W, XIE M, SHI W B, SUN S D, HOU S S. Developmental characteristics of pectoralis muscle in Pekin duck embryos., 2013, 12(4): 6733-6742.

[8] 胡艳, 刘宏祥, 单艳菊, 姬改革, 束婧婷, 徐文娟, 朱春红, 陶志云, 李慧芳. 鸭发育早期骨骼肌异步发育和IGF-1/MSTN-A表达的相关性. 中国农业科学, 2016, 49(2): 361-370.

HU Y, LIU H X, SHAN Y J, JI G G, SHU J T, XU W J, ZHU C H, TAO Z Y, LI H F. Correlation of the relative levels of insulin-like growth factor-1 and myostatin mRNA expression and asynchronous development of skeletal muscle development in ducks during early development., 2016, 49(2): 361-370. (in Chinese).

[9] HUANG Y, LI Y, BURT D W, CHEN H, ZHANG Y, QIAN W, KIM H, GAN S, ZHAO Y, LI J, YI K. The duck genome and transcriptome provide insight into an avian influenza virus reservoir species., 2013, 45(7): 776-783.

[10] TRAPNELL C, PACHTER L, SALZBERG S L. TopHat: discovering splice junctions with RNA-Seq., 2009, 25(9): 1105-1111.

[11] BENJAMINI Y, HOCHBERG Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing., 1995, 57(1): 289-300.

[12] YOUGN M D, WAKEFIELD M J, SMYTH G K, OSHLACK A. Gene ontology analysis for RNA-seq: Accounting for selection bias., 2010, 11(2): R14.

[13] HU Y, LIU H X, SONG C, XU W J, JI G G, ZHU C H, SHU J T, LI H F. Profiles of mRNA expression of related genes in the duck hypothalamus-pituitary growth axis during embryonic and early post-hatch development., 2015, 559(1): 38-43.

[14] GUERNEC A, BERRI C, CHEVALIER B, WACRENIER N, LE BIHAN-DUVAL E, DUCLOS M. Muscle development, insulin-like growth factor-I and myostatin mRNA levels in chickens selected for increased breast muscle yield., 2003, 13(1): 8-18.

[15] ECHTAY K S, WINKLER E, KLINGENBERG M. Coenzyme Q is an obligatory cofactor for uncoupling protein function., 2000, 408(6812): 609-613.

[16] TURUNEN M, OLSSON J, DALLNER G. Metabolism and function of coenzyme Q., 2004, 1660(1-2): 171-199.

[17] NISHITANI H, LYGEROU Z. Control of DNA replication licensing in a cell cycle., 2002, 7(6): 523-534.

[18] KEARSEY S E, LABIB K. MCM proteins: Evolution, properties, and role in DNA replication., 1998, 1398(2): 113-136.

[19] TYE B K. MCM proteins in DNA replication., 1999, 68(68): 649-686.

[20] ALLEN B L, UHLMANN F, GAUR L K, MULDER B A, POSEY K L, JONES L B, HARDIN S H. DNA recognition properties of the N-terminal DNA binding domain within the large subunit of replication factor C., 1998, 26(17): 3877-3882.

[21] UHLMANN F, GIBBS E, CAI J, O’DONNELL M, HURWITZ J. Identification of regions within the four small subunits of human replication factor C required for complex formation and DNA replication., 1997, 272(15): 10065-10071.

[22] ZHANG G, GIBBS E, KELMAN Z, DONNELL M O, HURWITZ J. Studies on the interactions between human replication factor C and human proliferating cell nuclear antigen., 1999, 96(5): 1869-1874.

[23] RAI M, KATTI P, NONGTHOMBA U. Spatio-temporal coordination of cell cycle exit, fusion and differentiation of adult muscle precursors by Drosophila Erect wing (Ewg)., 2016, 141: 109-118.

[24] LIU H H, WANG J W, LI L, HAN C C, HUANG K L, SI J M, HE H, XU F. Molecular evolutionary analysis of the duck MYOD gene family and its differential expression pattern in breast muscle development., 2011, 52(4): 423-431.

[25] BUCKINGHAM M, BAJARD L, CHANG T, DAUBAS P, HADCHOUEL J, MEILHAC S, MONTARRAS D, ROCANCOURT D, RELAIX F. The formation of skeletal muscle: from somite to limb., 2003, 202(1): 59-68.

[26] SCHIAFFINO S, MAMMUCARI C. Regulation of skeletal muscle growth by the IGF1-Akt/PKB pathway: Insights from genetic models., 2011, 1(1): 4.

[27] DELLING U, TURECKOVA J, LIM H W, de WINDT L J, ROTWEIN P, MOLKENTIN J D. A calcineurin-NFATc3-dependent pathway regulates skeletal muscle differentiation and slow myosin heavy-chain expression., 2000, 20(17): 6600-6611.

[28] LIU C, GERSCH R P, HAWKE T J, HADJIARGYROU. Silencing of Mustn1 inhibits myogenic fusion and differentiation., 2010, 298(5): C1100-C1108.

[29] XU T S, GU L H, SUN Y, ZHANG X H, YE B G, LIU X L, HOU S S. Characterization of MUSTN1 gene and its relationship with skeletal muscle development at postnatal stages in Pekin ducks., 2015, 14(2): 4448-4460.

[30] GOLDHAMER D, FAERMAN A, SHANI M, EMERSON C. Regulatory elements that control the lineage-specific expression of myoD., 1992, 256(5056): 538-542.

[31] BERKES C A, TAPSCOTT S J. MyoD and the transcriptional control of myogenesis., 2005, 16(4-5): 585-595.

[32] CHOI J, COSTA M L, MERMELSTEIN C S, CHAGAS C, HOLTZER S, HOLTZER H. MyoD converts primary dermal fibroblasts, chondroblasts, smooth muscle, and retinal pigmented epithelial cells into striated mononucleated myoblasts and multinucleated myotubes., 1990, 87(20): 7988-7992.

[33] DAVIS R L, WEINTRAUB H, LASSAR A B. Expression of a single transfected cDNA converts fibroblasts to myoblasts., 1987, 51(6): 987-1000.

[34] WEINTRAUB H, TAPSCOTT S J, DAVIS R L, THAYER M J, ADAM M A, LASSAR A B, MILLER A D. Activation of muscle- specific genes in pigment, nerve, fat, liver, and fibroblast cell lines by forced expression of MyoD., 1989, 86(14): 5434-5438.

[35] HASTY P, BRADLEY A, MORRIS J H, EDMONDSON D G, VENUTI J M, OLSON E N, KLEIN W H. Muscle deficiency and neonatal death in mice with a targeted mutation in the myogenin gene., 1993, 364(6437): 501-506.

[36] NABESHIMA Y, HANAOKA K, HAYASAKA M, ESUML E, LI S W, NONAKA I, NABESHIMA Y. Myogenin gene disruption results in perinatal lethality because of severe muscle defect., 1993, 364(6437): 532-535.

[37] ARIAS A M, STEWART A. Molecular principles of animal development. Oxford, UK:, 2002.

[38] KAMBADUR R, SHARMA M, SMITH T. Mutations in myostatin (GDF8) in double muscled Belgian Blue and Piedmon tese cattle., 1997, 7(9): 910-916.

[39] MCPHERRON A C, LEE S J. Double muscling in cattle due to mutations in the myostatin gene., 1997, 94(23): 12457-12461.

[40] GROBET L, MARTIN L, PONCELET D, PIROTTIN D, BROUWERS B, RIQUET J, SCHOEBERLEIN A, DUNNER S, MENISSIER F, MASSABANDA J, FRIES R, HANSET R, GEORGES M. A deletion in the bovine myostatin gene causes the double-muscled phenotype in cattle., 1997, 17(1): 71-74.

[41] MCPHERRON A C, LAWLER A M, LEE S J. Regulation of skeletal muscle mass in mice by a new TGF-β superfamily member., 1997, 387(6628): 83-90.

[42] CAMPOREZ J P G, PETERSEN M C, ABUDUKADIER A, MOREIRA G V, JURCZAK M J, FRIEDMAN G, HAQQ C M, PETERSEN K F, SHULMAN G I. Anti-myostatin antibody increases muscle mass and strength and improves insulin sensitivity in old mice., 2016, 113(8): 2212-2217.

[43] WHITTEMORE L-A, SONG K N, Li X P, AGHAJANIAN J, DAVIES M, GIRGENRATH S, HILL J J, JALENAK M, KELLEY P, KNIGHT A, MAYLOR R, O'HARA D, PEARSON A, QUAZI A, RYERSON S, TAN X Y, TOMKINSON K N, VELDMAN G M, WIDOM A, WRIGHT J F, WUDYKA S, ZHAO L, WOLFMAN N M. Inhibition of myostatin in adult mice increases skeletal muscle mass and strength., 2003, 300(4): 965-971.

（责任编辑 林鉴非）

RNA-seq Analysis on Development Arrest of Duck Pectoralis Muscle During Semi-Late Embryonic Period

LIU HongXiang, XU WenJuan, ZHU ChunHong, TAO ZhiYun, SONG WeiTao, ZHANG ShuangJie, LI HuiFang

（Jiangsu Institute of Poultry Sciences, Yangzhou 225125, Jiangsu）

【Objective】In order to find the molecular varying mechanism involving the development arrest of pectoralis, Chinese native breeds, including Gaoyou Duck (GYD) and Jinding Duck (JDD), were selected for RNA-seq study using the pectoralis samples during the semi-late embryonic period.【Method】3 ducks of GYD and JDD , respectively, in the 21th embryonic day (ed21) and ed27 were selected to collect pectoralis major muscle and to extract total RNA to build cDNA library for RNA-seq with HiseqTM2000 of Illumina. At last, bioinformatics methods were used to extract differentially expressed genes (DEGs) between different breeds and time points, and to analyze the gene function annotation for studying molecular mechanism of pectoralis development retardation during ed21 and ed27. 【Result】The results showed that the base ratios with quality value higher than 20 in reads (Q20) were more than 94%, and the base ratios of Q30 were more than 89%, which indicated a reliable sequencing result for the following analysis. The RNA level correlation inspection and mRNA expression level cluster graph both manifest that the correlation of mRNA expression patterns of GYD and JDD at ed21 or ed27 were higher than that of GYD (JDD) during ed21 and ed27. The numbers of DEGs between ed21 and ed27 (6128 DEGs for GYD and 6452 DEGs for JDD) were both apparently more than the numbers of DEGs between GYD and JDD in ed21 (522 DEGs) and ed27 (299 DEGs). qRT-PCR results of selected genes showed a strong correlation with RNA-seq results. GO and KEGG enrichment analysis showed the results that the genes involved in energy metabolism (mainly was coenzyme Q related genes, ATP enzymic synthesis related genes, and cytochrome C related genes) up regulated and DNA replication and cell cycle related-genes (mainly was minichromosome maintenance complex related genes and replication factor C related genes) down regulated significantly. This varies of related genes expression might relate to the slow myoblast proliferation and gradually exit the cell cycle to prepare for the next stage of fusing to multi-nuclei myotube and form myofiber. In the key genes involving in muscle growth and development,1 (a gene promoting muscle growth) and(inducing terminal differentiation of myoblast) both down regulated dramatically. However, MUSTN1 gene, accelerating muscle fibre into differentiation and fusing stage, andgene, promoting muscle progenitor cell differentiated to myoblast, were expressed in a higher level in ed27 than in ed21. 【Conclusion】A lot of genes differentially expressed between ed21 and ed27 in pectoralis muscle of duck, among of which the up-regulation of energy metabolism related genes, the down-regulation of DNA replicate and cell cycle related genes, and up-regulated1, down-regulated1 and, might closed relate to arrest phenomenon of pectoralis development during the semi-late stage of duck embryos.

duck; pectoralis muscle development; RNA-seq; differentially expressed genes

2017-08-03；

2018-09-12

国家自然科学基金（31172194）、江苏现代农业（水禽）产业技术体系建设专项资金（JATSI2018J225）、江苏现代农业重点项目（BE2017349）

刘宏祥，E-mail：lhxatyz@foxmail.com。徐文娟，E-mail：xuwj1980@126.com。刘宏祥和徐文娟为同等贡献作者。

李慧芳，E-mail：lhfxf_002@aliyun.com.cn

10.3864/j.issn.0578-1752.2018.22.015