基于全转录组测序的绵羊胚胎不同发育阶段骨骼肌circRNA的...

奥鹏网院作业

基于全转录组测序的绵羊胚胎不同发育阶段骨骼肌circRNA的分析与鉴定
石田培，王欣悦，侯浩宾，赵志达，尚明玉，张莉

（中国农业科学院北京畜牧兽医研究所，北京 100193）

摘要：【目的】肉用性状是家养动物的重要性状，尤其是产肉性能与骨骼肌的发育密切相关，对多数动物而言，骨骼肌的生长取决于肌纤维的早期发育。本研究将揭示绵羊（Ovis aries）胚胎期骨骼肌组织发育重要节点、肌纤维的形成及转换机制，探究骨骼肌在动物胚胎期的生长对后续发育潜力的影响。【方法】在前期绵羊胚胎期组织结构试验的基础上，选择肌纤维发育及类型转换的3个重要时间节点：胎龄85 d（D85）、105 d（D105）、135 d（D135），并对处于这些节点期的中国美利奴绵羊胚胎背最长肌进行全转录组测序。通过生物信息学分析和功能预测筛选出差异表达显著的circRNA，采用实时荧光定量 PCR（qRT-PCR）对其进行验证。【结果】通过条件筛选（|log2| ≥1且P≤0.05）获得差异表达circRNA 1 126个。将3组进行比较，在各时间节点都发现较多的特异性表达circRNA。在D85 vs D135比较组中，特异性差异表达数量最多：共发现差异表达circRNA 374个，上调表达201个，下调表达173个，其中具有4倍以上差异的基因167个，占差异基因的44.7%。对这些差异性表达的circRNA进行GO和KEGG功能分析和靶向预测，富集到与肌肉分化和肌纤维发育相关的能量代谢和信号转导通路，包括MAPK、PI3K-Akt、Ras、Regulation of actin cytoskeleton等。结果表明：在D85至D105期间富集到的差异circRNA多与肌细胞发育调控、细胞增殖和存活、细胞周期等相关，而在D105至D135期间富集到的通路则主要与能量转换、物质运输、RNA转运和DNA修复有关。靶向预测结果通过Cytoscape软件绘制出可视化共表达网络，筛选到circRNA8239、circRNA19073，circRNA2765、circRNA1616等核心调控转录物，并在D105找到调节快慢肌类型转换的关键因子circRNA7527，其靶向作用于bta-miR-135a、bta-miR-615、chi-miR-133a-5p进而调控MEF2C基因。通过3个时间节点的差异表达情况和功能预测描述选取了与肌肉发育相关的4个核心circRNA和其靶向的4个miRNA进行qRT-PCR分析，其基因表达趋势与测序数据一致。【结论】本试验在转录水平上验证了绵羊胚胎发育D85至D105是肌纤维数量的稳定发生期，而D105至D135则是肌纤维的肥大期，由此推断D105可能是绵羊胚胎发育的关键时间窗口。本研究首次基于全转录组测序构建了绵羊胚胎期骨骼肌发育的circRNA图谱，揭示了不同发育阶段的转录组差异，并对绵羊胚胎骨骼肌发育期间重要的调控因子进行挖掘和验证，发现了以MEF2C为靶基因的多个circRNA和miRNA参与MAPK信号通路，为绵羊肌纤维发育分子机制以及其它家畜非编码RNA的研究提供了参考。

关键词：绵羊；胚胎；全转录组；骨骼肌；生长发育；circRNA

0 引言
【研究意义】肉用性状是家畜的重要经济性状，研究骨骼肌发生机制有助于提升家畜产肉性能及肌肉品质。骨骼肌的形成根据发育特点可以分为两个阶段：成肌细胞的增殖、迁移和分化，以及肌纤维的延伸、增粗和类型转换。研究表明成肌细胞在动物妊娠中后期融合形成梭形排列的肌管，最终形成终末分化的肌原纤维[1]。家畜肌纤维数量在妊娠后期达到一定数目后就不再发生变化，而成年后肌肉量的增加主要依靠肌纤维长度和周径的增加。由此可见，胚胎期骨骼肌的生长和发育对家畜的肉用性状至关重要。【前人研究进展】目前，已有大量关于骨骼肌转录组和蛋白质组发生机制的研究，一些蛋白质和转录物作为调节因子在胚胎肌细胞发育过程中发挥着重要作用，如肌源性调节因子（myogenic regulatory factors, MRFs）通过肌源性因子5（myogenic factor 5, Myf5）、肌源性分化1（myogenic differentiation 1, MyoD）、肌源性调节因子4（myogenic regulatory factor 4, Mrf4）和肌细胞生成素（Myogenin）等协同调控胚胎骨骼肌的发育和分化[2-3]，其中每个MRF都可单独作为肌生成的主要调节因子，精细调节胚胎骨骼肌的发育过程[4]。Six家族蛋白质作为转录因子存在两个特征性结构域：一个是可与DNA结合的Six型同源结构域，另一个是与转录激活或抑制因子相互作用的氨基末端结构域。Six1（The sine oculis–related homeobox 1）和Six4通常被认为是遗传调节级联的制高点[5-6]，能够诱导肌源性祖细胞发育成为肌源性细胞系。除经典调控因子外，Pax7能与Wdr5-Ash2L-MLL2组蛋白甲基转移酶（HMT）复合物结合，在Myf5的作用下，致使染色质周围的H3K4三甲基化[7]。因此，通过Pax因子富集HMT复合物的途径被普遍认为是可能重塑染色质结构以控制谱系特异性基因表达的保守机制[8]。近年来，有关骨骼肌发育非编码RNA的研究较多[9-11]，但受样品、设备及分析方法等因素的限制，在全转录组水平上同时对多种非编码RNA的关联分析还未见报道。【本研究切入点】肌肉形成与发生机制是一个复杂的动态过程，目前对该生理过程还缺乏系统的研究。本研究将使用二代测序技术，从基因转录组层面对肌肉发生机制进行系统的、全面的特征化和量化分析[12]。【拟解决的关键问题】应用RNA sequencing技术对绵羊不同发育时期的胚胎骨骼肌进行全转录组测序和生物信息学分析，对骨骼肌早期发育过程中的关键调控机制进行探究，发掘重要的功能性调控非编码RNA。

1 材料与方法
1.1 动物道德伦理声明
本研究方法严格按照中华人民共和国农业农村部制定的相关指导方针和规定进行，所有试验方案均经中国农业科学院北京畜牧兽医研究所批准。

1.2 实验动物和样品制备
试验选择年龄相近、体况良好、体重在 55—60 kg的成年中国美利奴羊（购自新疆巴州种畜场），经埋栓同期发情处理后进行人工授精配种，通过超数排卵获得胚胎并移植到体况相近的受体母羊，妊娠母羊均按照国家农业行业肉羊饲养标准（NY/T816- 2004）进行饲喂。胚胎分别在妊娠第85天、105天和135天时通过外科手术取出（每组采样个体为 3 只），采集背最长肌，并分成小块于-80℃冰箱保存备用。

1.3 方法
1.3.1 RNA 的提取及质控 准备足量样品，按照Trizol试剂盒提取样品总RNA，提取后使用微量核酸蛋白测定仪（scandrop100）检测OD 260 nm /OD 280 nm 值以鉴定RNA样品浓度，排除RNA污染可能性。Bioanalyzer 2100（Agilent, Santa Clara, CA）用于评估总 RNA 质量，以 RIN≥7 且260/280≥1.8为阈值，RNase-free DNase I (Ambion Inc., Texas, USA)用于消除潜在的基因组 DNA 污染。分离的RNA 样本保存在-80℃环境下。

1.3.2 cDNA 文库构建和 RNA 测序 RNA质检合格后，分别从3个阶段肌肉样品中各取出约10 μg RNA，分别构建测序文库，主要包括酶促RNA 片段化、cDNA 合成、测序接头衔接及 PCR 扩增等过程。首先，根据EpicentreRibo-Zero Gold试剂盒（Illumina，San Diego，USA）的操作说明，消耗样品中的核糖体RNA、线性RNA并纯化后，在高温下使其片段化，随后转录成cDNA文库，最后分别利用Illumina Hiseq 4000进行测序。

1.3.3 数据处理及转录水平分析 对测序产生的原始数据进行Cutadapt[13]过滤，除去低质量、含N的序列后通过Bowtie2[14]和Tophat2[15]将读数映射到绵羊的Oar 3.1基因组，其中获得的注释文件包含转录本的分布信息。使用StringTie[16]识别新的转录本，并在此基础上利用Ballgown [17]估计转录物的表达水平。然后通过FPKM算法消除偏好影响，进行基因表达定量。在本研究中，采用edgeR包确定差异表达基因[18]，同时使用上四分位数算法校正数据，即设定| log2（两个样本的FPKM比率）| ≥1且P≤0.05为阈值。TopHat- fusion[19]和CIRCExplorer2[20]可以在未映射序列中识别反向剪接读数来估计circRNA的表达水平[21]，circRNA及其相应亲本基因的共表达关系通过Pearson相关系数R2值估算确定。miRNA序列的处理需使用专有脚本，同时使用多个过滤器除去不可映射的测序读数。通过比对miRBase 21.0中pre-miRNA（mir）和成熟miRNA（miR）的独特序列或已公开数据库进行“映射”，剩余序列进行BLAST。最后通过RNA折叠软件从侧翼80nt序列处预测含有发夹RNA结构的序列。本研究通过TMeV软件对差异表达的circRNA进行热图绘制，直观展示各样本的表达一般规律、差别及变化，从而通过聚类预测新型circRNA的功能及作用。CircRNA-miRNA互作分析主要采用Targetscan与miRanda软件进行靶标预测后取交集，将得到的数据上传至Cytoscape绘制可视化互作网络图。

1.3.4 差异表达基因 GO 及 KEGG 通路富集分析 经定量确定的差异表达RNA仍需深入了解其生物学功能以发掘重要的调控因子，在生物信息学分析中，常采用基因本体论（gene ontology，GO）[22]进行分析。在GO体系中，可以依次从细胞成分、生物学过程和分子功能3个方面定位基因功能，使用其注释富集功能将具有相似注释的转录物整合归类。在预测各类未知RNA时，除了通过RNA的宿主基因进行功能预测外，后期还可间接通过它们在生物信号通路中的参与情况进行功能分析。生物通路富集分析主要依赖于京都基因和基因组百科全书（kyoto encyclopedia of genes and genomes , KEGG）[23]数据库，其中途径富集分析可以精确定位差异表达基因的主要代谢转导和信号传导途径。本研究使用Scripts in house软件进行GO（http://geneontology.org）和KEGG（http://www.kegg. jp/kegg）富集分析，设置统计P≤0.05为结果显著阈值。

1.3.5 实时荧光定量 PCR 验证 为确保测序结果的准确性和有效性，需要通过qRT-PCR 验证基因表达水平。选择验证的RNA至少在一组对照组中表达差异显著，且与骨骼肌的生长发育相关。经过组间对比和筛选，最终分别选择circRNA和miRNA各4个进行验证。根据转录本数据库设计转录物引物（表1），使用CFX-Connect TM实时系统（BioRad）上的SYBR green I master mix（Roche，GmBH，Basel，Switzerland）进行实时PCR。试验结果用目的基因值减去内参基因值得到ΔCt，再以相对参照样品的ΔCt均值为对照组，用各组ΔCt值减去对照组得到ΔΔCt，以2-ΔΔCt 平均值表达出各个基因的相对表达量，利用SPSS 22.0软件对所得表达量进行单因素方差分析(One-way ANOVA)，结果分为差异显著（P≤0.05）和差异极显著（P≤0.01）。

表1 荧光定量引物基因列表

Table 1 The list of gene primers of qRT-PCR


2 结果
2.1 RNA-seq 数据分析
采用RNA-Seq和Small RNA-Seq深度测序，每个样本文库内获得约90 000 000的原始数据，质量控制后平均得到约 88 000 000的有效读段。Q-score＞20的reads 在 99%以上，超过97%的reads Q-score＞30。在有效读段中，约85% reads比对到绵羊参考基因组上，其中超过82%是唯一比对，除此之外还进行了非剪接/剪接比对计数和正/反义链计数。综上，充分证明样本数据处理得当，质量可靠，符合下游数据分析要求。在深度分析样本注释结果，评估骨骼肌中转录本丰度及表达差异情况前，计算评估了各生物学重复间皮尔森积矩相关系数平方（R2）均大于 0.98，均满足统计学要求。

2.2 差异表达转录本鉴定
在3个对照组中，共鉴定37 474条circRNA，利用 FPKM 算法标准化基因转录水平差异，筛选出显著差异circRNA 1 126条。在转录本识别过程中，circRNA表达量FPKM值结果显示转录本表达水平存在明显差异，但仍有相当一部分没有表达或者表达量极低（图1）。对3组数据进行比对，circRNA的表达情况表明各阶段都存在特异性表达的非编码RNA，但在D85 vs D135的发育阶段中，特异性表达数量最多。两组比对出374个差异表达基因，其中上调表达201个，下调表达173个，具有4倍以上差异的基因167个，约占差异基因的44.7%（图2）。miRNA识别元件（miRNA response element, MRE）的存在可介导miRNA与靶标转录物的部分互补序列结合，导致机体相应靶标基因的表达抑制。因此，通过TargetScan和miRanda预测miRNA和mRNA之间的靶向作用，共找到76 759个互作关系，根据GO term和KEGG pathway功能描述，筛选出与骨骼肌发育相关的互作关系497对。在D85 vs D105 vs D135的比较组中，将基因按照GO匹配显著度排名，其中GJA1、NACA和PAX7极为显著。

2.3 差异表达circRNA的功能富集分析
基于GO数据库的注释结果表明，Unigenes中共有20 840条转录本得到归类注释，显著富集于353个生物学过程条目，117个分子功能条目，144个细胞组分条目，主要富集于纤维组装、细胞增殖与迁移、细胞分裂、骨骼肌收缩等生物学过程。3个阶段胚胎背最长肌circRNA显著富集于124个生物学过程条目，36个分子功能条目和40个细胞组分条目。3个发育阶段都富集到与骨骼肌发育相关的生物学过程，但显著高表达基因注释结果略有差异（图3）。在D85 vs D105中，信号转导调节、转录活性调节和骨骼肌细胞分化相关的基因富集最为明显；在D105 vs D135中，肌纤维发育、肌细胞增殖和胚胎发育调节相关的基因富集最为明显；在D85 vs D135中，显著富集并与其它阶段重合的通路有细胞增殖与转运，肌细胞分化和信号转导等。通过KEGG 通路富集分析，发现差异表达基因在多个信号转导通路中被显著富集，其中富集于心肌病、Wnt信号通路和Ras信号通路的基因比例较大。结合RNA的GO分析和 KEGG分析结果发现，绵羊胚胎骨骼肌的发育同时受到多种机制的调控，在不同阶段涉及的信号通路略有差别。差异表达基因被显著富集于69个通路，其中涵盖了多个基因。如图4所示，MAPK、PI3K-Akt、Ras、Regulation of actin cytoskeleton等与肌肉发育相关能量代谢及转运通路显著富集。在所有被富集途径中，PI3K-Akt信号传导途径是参与细胞凋亡、蛋白质合成和细胞周期等细胞基本功能的较活跃途径。在信号系统中常通过信使作用产生级联反应，首先生长因子与受体酪氨酸激酶（RTK）或G蛋白偶联受体（GPCR）结合分别刺激Ia和Ib类PI3K活化，继而催化细胞膜上磷脂酰肌醇-3,4,5-三磷酸（PIP3）的产生。然后PIP3作为第二个信使激活Akt，Akt便可以通过磷酸化调控关键的细胞过程。活化产生的PIP激活了PDK1、Akt等第二信使导致PKCs、PKN[24]、SGK等在细胞增殖过程中上调表达，Raf-1[25]、ERK等下调表达，稳定维持肌细胞的增殖合成，血管生成以及DNA损伤修复。在绘制的信号转导通路中清晰可见，大部分与细胞发育增殖、蛋白合成和DNA损伤修复相关的基因表达显著，呈现橙色（图4）。此外，通过对D85 vs D105、D105 vs D135和D85 vs D135等3组差异表达circRNA进行聚类筛选，最终挑选出37个具有骨骼肌细胞分化、骨骼肌纤维发育和骨骼肌细丝组装等功能描述的circRNA，根据转录水平差异对其绘制聚类热图（图5），由表达模式聚类可预测未能得到注释的转录本。在热图中，这些circRNA都至少在一个时间节点显著上调或下调。

width=226.8,height=249.5
图1 不同样本circRNA的FPKM值分布

Fig. 1 Distribution of circRNAs based on FPRM in various samples

width=206.5,height=131.35
图2 不同样本circRNA表达统计图

Fig. 2 Distribution of circRNAs based on FPRM in various samples

width=407.15,height=373.35
width=386.85,height=295.15

width=332.45,height=309.45

A：D85 vs D105；B：D105 vs D135；C：D85 vs D135；Y轴：pathway名称；X轴：富集因子；气泡大小：富集到某个通路中差异表达的数量/背景中所有数量；颜色反映了富集在某通路差异表达数及富集显著性

A: D85 vs D105; B: D105 vs D135; C: D85 vs D135; Y-axis: pathway name; X-axis: enrichment factor; Bubble size reflect the number of differential expressions enriched to a certain pathway/all quantities in the background; Color reflect the significance of differential expression in a certain process

图3 差异表达circRNA的KEGG通路分析

Fig. 3 KEGG pathway analysis of differentially expressed circRNAs

2.4 差异表达miRNA-circRNA互作网络
在转录调控中，网络调控分析可以鉴定出潜在调控重要性状的关键候选基因[26]。有研究表明circRNA可作为参与miRNA表达调控的竞争性内源RNA（competing endogenous RNA, ceRNA），本研究采用Targetscan和miRanda软件来分析circRNA-miRNA相互作用关系。Targetscan基于种子区域进行miRNA目标预测，而miRanda基于circRNA和miRNA的自由能的组合，根据软件评分结果，确定了877个显著差异circRNA，它们均与多个miRNA之间存在互作，共鉴定到6 340个miRNA-circRNA相互作用对，其中194个相互作用对与骨骼肌发育相关，我们将这些与骨骼肌发育相关的miRNA和circRNA构建了相互作用网络。在网络图中，某一节点与另一节点的连线称之为度，度的大小可直观反映基因在网络中的重要程度，拥有相对较大度值的基因可能具有更重要的生物学价值。计算共表达网络中所有基因的度，根据共表达网络中所有基因度值排名，提取网络中连接度较大的circRNA、miRNA，对其表达情况进行验证。circRNA5502、circRNA5505、circRNA5561、circRNA30828和circRNA23584在互作网络中度值较大，表明其宿主基因可能在骨骼肌细胞生长代谢及发育过程中起重要作用（图6）。

width=484.25,height=401.95
白色小框是根据已有的数据绘制的具有一般参考意义的，概括详尽的代谢图；绿色小框为该物种特有的基因或酶，具有更详细的信息；橙色方框表示该基因在此时差异性表达

The white box is a generalized and detailed metabolic map based on the existing data; the greenbox is the unique gene or enzyme of the species with more detailed information; the orange box indicates that the gene differential expressed in here

图4 PI3K-Akt信号通路

Fig. 4 PI3K-Akt Signal Pathway

2.5 构建circRNA-miRNA-mRNA共表达调控网络
将筛选得到的与肌肉发育及肌细胞增殖分化相关的差异表达circRNA及其靶向的miRNA和mRNA进行转录调控分析和共表达网络分析。在互作网络中，4个特异性差异表达circRNA和其互作的miRNA与mRNA共形成208个互作关系，根据度值发现MAPK9、IGF2、MYO5A等在网络中具有重要调控地位（图7）。TTN表达的circRNA1616受chi-miR- 30f-3p反式调控作用于MYO5A和MYO5C，MEF2C表达的circRNA19073，circRNA2765作用于MYO5C、DMD、IGF2、MYO5A等。结合三者互作的GO分析结果（图8），发现这些表达性基因主要被富集在胞质、核质和ATP结合等能量代谢相关的生物学过程，KEGG通路中显著富集到的通路都是参与代谢和细胞生长的途径，如Ubiquitin mediated proteolysis、Wnt signaling pathway、Insulin signaling pathway、Focal adhesion、Cell cycle。根据以上结果可推测，在骨骼肌发育的D85至D135过程中，骨骼肌发生重点进行了转移，circRNA后期主要参与能量代谢和物质运输的调控。

width=420.4,height=529.95
横轴：胚胎日龄，纵轴：circRNA名称；蓝色：下调，红色：上调

Horizontal axis: embryo age; Vertical axis: circRNA name; Blue: down; Red: up

图5 差异表达circRNA热图

Fig. 5 The heat map of differentially expressed

2.6 实时荧光定量检测
实时荧光定量结果表明，获得的所有基因的溶解曲线均为单峰，说明定量试验设计的引物较优，并且在定量检测过程中不存在非特异性扩增与引物二聚体（图9）。通过对8条RNA的qPCR检测结果进行校正后，检测差异表达水平与深度测序结果一致，说明高通量测序数据可靠。继续对功能聚类后筛选出的关键RNA进行分析，以D85 vs D105、D105 vs D135和D85 vs D135进行组间显著性检验。设定P≤0.05为显著，P≤0.01代表极显著，差异具有统计学意义。结果显示：circRNA17939、chi-miR-133a-5p和oar-miR- 150在组间差异极显著，circRNA17939在D85 vs D105阶段差异极显著，bta-miR-365-3p在D85 vs D105阶段差异显著，在D105 vs D135差异极显著，符合数据鉴定和功能分析的结果。

3 讨论
本项目前期动物试验表明：在美利奴绵羊胚胎期，从D75至D135，各体尺指标（体重、体长和体高）均随日龄增加而增长。在D75骨骼肌可见中空管状的初级肌纤维典型结构，且在初级肌纤维四周围绕着大量的次级肌纤维；D105 和D135切片中只有次级肌纤维结构而无初级肌纤维结构特征。可见在D75后，初级肌纤维与次级肌纤维发生融合，肌纤维的形成基本完成；在D105时肌纤维数量基本恒定；到D135时肌纤维的发育均表现为肌纤维面积不同程度的增大，肌纤维直径显著增加[9]。李雪娇等在绵羊胎儿骨骼肌组织学结构发育特征研究中表明胎儿D105是肌纤维从肌增生到肌肥大的关键时间[27]。基于此，对绵羊胚胎D75至D105之间进行了进一步细化分析，结果发现D75至D85肌纤维开始相互融合，D85时已无显著特征的初级肌纤维，形成成熟的肌纤维结构，肌纤维数量达到整个胎儿发育阶段的峰值。因此选择D85、D105 和D135时间节点构建绵羊骨骼肌发育阶段的全转录组文库，通过 RNA-seq 测序和生信分析筛选出在转录水平上可能参与调控的RNA，结果共得到显著差异表达circRNA1 126 个。对筛选结果进行GO和KEGG富集分析，发现3个阶段都显著富集到大量与肌肉发育相关的基因和通路，如RNA信号转导、物质代谢和转录活性通路等。GO分析结果显示，circRNA及其互作的miRNA和mRNA在骨骼肌的生长发育中都具有重要功能，主要包括肌细胞发育（GO:0055001）、骨骼肌卫星细胞分化（GO:0014816）、肌纤维发育（GO:0048747）、骨骼肌收缩（GO:0003009）、骨骼肌纤维发育（GO:0006979）和骨骼肌组织生长的积极调节（GO:0048633）等生物学过程。在D85 vs D105比较组中，鉴定到circRNA23653、circRNA19073和circRNA6142等与肌细胞增殖和调控相关的circRNA，还发现circRNA7527与慢速和快速纤维之间的转换密切相关。在D85 vs D135比较组中，发现与肌纤维发育相关的circRNA15个：circRNA1615、circRNA5174、circRNA3179、circRNA15129、circRNA7311、circRNA1616、circRNA6291、circRNA2324、circRNA18785、circRNA19839、circRNA1602、circRNA19702、circRNA31369、circRNA7569、circRNA489。比较分析发现circRNA7527未在D85 vs D135差异表达，说明circRNA7527只在D105时期特异性表达，推测其功能可能是调控肌纤维的转换。根据转录水平的对比，筛选到MSTN、MYOD、MYOG、MYF5、PAX7、NF1、MEGF10、GDF11、PLAGL1、S100B、MEF2C、MEF2D等15个重要基因，涉及到肌纤维发育的信号通路包括TGFβ、MAPK、PI3K- Akt、Ras、Akt/mTOR等。circRNA18785、circRNA19839、circRNA1602、circRNA19702、circRNA31369、circRNA7569、circRNA489。比较分析发现circRNA7527未在D85 vs D135差异表达，说明circRNA7527只在D105时期特异性表达，推测其功能可能是调控肌纤维的转换。根据转录水平的对比，筛选到MSTN、MYOD、MYOG、MYF5、PAX7、NF1、MEGF10、GDF11、PLAGL1、S100B、MEF2C、MEF2D等15个重要基因，涉及到肌纤维发育的信号通路包括TGFβ、MAPK、PI3K- Akt、Ras、Akt/mTOR等。

width=484,height=415.1
橙色椭圆表示差异表达的miRNA；绿色箭头表示差异表达的circRNA；边表示两者之间的互相作用关系。颜色由浅到深代表表达差异程度，颜色越深差异越显著，反之则显著程度降低

Orange ellipses indicate differentially expressed miRNAs; Green arrows indicate differentially expressed circRNAs; and Edges indicate interactions between the two. The color from light to dark represent the degree of difference in expression, and the deeper the color, the more significant the difference, and vice versa

图6 miRNA-circRNA互作网络图

Fig. 6 miRNA-circRNA interaction network diagram

width=480.3,height=318.6
橙色椭圆表示差异表达的mRNA；蓝色箭头表示差异表达的circRNA；红色菱形表示差异表达的miRNA；边表示两者之间的互相作用关系

Orange ellipses indicate differentially expressed mRNA; Blue arrows indicate differentially expressed circRNA; Red diamonds indicate differentially expressed miRNAs; and Edges indicate interactions between the two

图7 circRNA- miRNA-mRNA互作网络图

Fig. 7 circRNA-miRNA-mRNA interaction network map

width=332.2,height=322.4
Y轴：生物学过程；X轴：富集因子；气泡大小反映了富集到某个GO_Term的差异表达的数量/背景中所有数量；颜色反映了富集在某过程差异表达数及富集显著性

Y-axis: biological process; X-axis: enrichment factor; Bubble size reflect the number of differential expressions enriched to a certain GO_Term/all quantities in the background; Color reflect the significance of differential expression in a certain process

图8 circRNA-miRNA-mRNA互作网络的GO富集分析

Fig. 8 GO analysis of circRNA-miRNA-mRNA

成熟的miRNA通过特定序列互补配对识别靶mRNA，并根据互补程度沉默复合体降解mRNA或阻抑mRNA的翻译。与miRNA单一的作用机制不同，circRNA可从转录、转录后及表观遗传学等多个水平参与基因的调控。circRNA可以包含核糖体进入位点，可以翻译表达有效蛋白，转录时与mRNA前体发生剪切竞争，还可与基因上游启动子序列结合调控下游基因表达，作为相关分子元件与特定蛋白结合调节相应蛋白活性，因此circRNA的功能研究更为复杂。本研究中发现多个circRNA通过靶向调控或吸附miRNA而顺式调控其宿主基因的表达。在D105时期找到关键性的circRNA7527，并在bta-miR-135a、bta-miR-615、chi-miR-133a-5p和mmu-mir-6240-p5_5序列中发现可与其结合的位点，此外，根据功能性靶向预测，发现这些miRNA都可作用于MEF2C。在D105 至D135时期，bta-miR-135a、bta-miR-615和mmu-mir-6240- p5_5随着circRNA的下调而上调，但chi-miR-133a-5p却随之下调，说明bta-miR-135a、bta-miR-615和mmu- mir-6240-p5_5受到的调控机制可能与chi-miR-133a- 5p不完全相同。除此之外，预测出的miRNA（bta-miR- 135a、bta-miR-615和mmu-mir-6240-p5_5）的靶基因MEF2C上下调水平呈现相反趋势，该表达规律符合哺乳动物体内miRNA与mRNA的作用机制。但chi-miR-133a-5p与MEF2C的下调情况相同，说明其之间可能不存在靶向关系。有研究证明在胚胎发生过程中心肌和骨骼肌组织分化开始时，MEF2C在MEF2家族其他基因开始表达后随即表达，但单个MEF2C的功能缺失会造成肌细胞的分化[28]。由此可见circRNA7527、bta-miR-135a、bta-miR-615、chi-miR- 133a-5p和mmu-mir-6240-p5_5中的一个或几个对肌纤维的分化具有重要调控意义。在比较组中还发现circRNA2324可与bta-miR-219作用靶向调节NF1，但该circRNA只在D85和D105两个阶段上调表达，在D135时显著下调，说明在胚胎骨骼肌发育后期，circRNA2324的作用受到抑制。NF1[29-30]编码神经纤维素蛋白，具有Ras特异性GTP酶激活结构域，可作为Ras /丝裂原活化蛋白激酶（MAPK）信号传导的负调节物。由于胞内表达NF1可诱导细胞退出终止细胞周期，导致肌细胞终末分化[31]。所以早期间充质中NF1的失活可导致肌纤维数量骤减，机体全部肌肉纤维化[32]。此外，NF1的缺失会造成MEK1的减少，而MEK1是早期负向调节肌肉分化的重要因子，MEK1与MyoD核内复合物相互作用并抑制其功能[33]。MEK1上游的Ras和/或Raf激酶的组成型活化可以抑制肌细胞分化[34-35]，Ras的下游MAPK信号传导参与肌原性分化的调节[36]，这与之前进行的功能性pathway分析结果相符。

width=438,height=485.7
图中为RNA-seq结果和qRT-PCR的定量结果；折线图表示测序的FPKM值，柱状图表示相对定量结果；显著性检验：D105柱状图上标符号为D85 vs D105比较结果；D135柱状图上2个标符号，左边为D85 vs D135比较结果，右边为D105 vs D135比较结果；* P≤0.05，** P≤0.01

The figure shows the quantitative results of RNA-seq and qRT-PCR; the line graph shows the FPKM value of the sequencing and the histogram shows the relative quantitative results; the significance test: the D105 histogram is marked with the D85 vs D105 comparison result; 2 test result in D135 histogram: left is D85 vs D135 and right is D105 vs D135; *P≤0.05, ** P≤0.01

图9 实时荧光定量检测

Fig. 9 The results of qRT-PCR

KEGG富集分析发现大多数差异表达circRNA主要富集在与肌细胞发育周期和分化的生理过程中，介导了细胞对外界的反应。p38MAPK、JNK、ERK1/2和ERK5/BMK1是MAPK通路中重要亚族[37-38]，有研究证明它们是调控肌肉细胞和脂肪细胞增殖和分化的重要候选基因[39]。在绘制的信号通路图中，发现由激活的第二信使调控的PKCs、PKN、SGK等基因在细胞增殖、存活和肌动蛋白的组装过程中起着重要作用。根据差异表达情况可知，这3个基因在该通路中的受调控作用显著（图4）。PKN常常通过与骨骼肌和非骨骼肌型α-辅肌动蛋白的第三血影蛋白重复序列结合，并以Ca2 +不敏感的方式与骨骼肌型α-辅肌动蛋白的的EF-手状基序的区域结合而直接关联与细胞骨架网络的连接，在肌动蛋白骨架组装中具有重要作用。Ras蛋白通过在GDP结合状态和GTP结合状态间循环转换来转导胞外生长因子的信号，活化的Ras（Ras-GTP）主要可激活两种蛋白激酶，即Raf-1和MEK激酶（MEKK）。激活丝裂原活化蛋白激酶（MAPK），包括细胞外信号调节激酶（ERK）和Jun激酶（JNK），活化的Raf-1直接促进ERK活化但不促进JNK活化，而MEKK参与JNK活化，但仅在过表达后引起ERK活化，即两种不同的Ras依赖性MAPK级联：一个由Raf-1引发导致的ERK活化，另一个是由MEKK引发导致的JNK活化。在富集到的PI3K-Akt，Rap1和Ras信号通路中，显著体现了Raf-1引发的ERK活化，促使基因表达的级联模式。以上结果均表明MAPK等信号通路是参与调控肌细胞增殖分化的关键通路，直接或间接地参与绵羊胚胎肌肉的发育调控。

值得一提的是，LIU等对初生大白猪（Sus domesticus）的肌肉进行了转录组和蛋白组共分析，发现肌肉中差异表达基因主要被富集在分泌糖蛋白的Wnt、MAPK、TGF-β和 cAMP信号通路中，并验证了TGF-β可通过促进成纤维细胞增殖和神经细胞分化而影响肌肉发育及出生后肌肉肥大[40]，这与本研究中被富集到的显著通路和纤维生长关键基因基本相同。VALENTIN等对不同胚胎期的大白猪骨骼肌发育机制进行了多组学关联分析，结果发现骨骼肌发生后期至肌肉成熟阶段，转录物和蛋白质的功能主要集中在物质能量代谢（氧化磷酸化等）中[41]。此外，HONG等在家禽（Gallus domesticus）中也有相似的胚胎阶段性研究，鸡胚的E11、E16及D1阶段蛋白质互作网络分析表明所有差异表达蛋白质主要参与蛋白质合成、能量代谢、肌肉收缩和氧化磷酸化过程[42]。本研究结果与上述猪和家禽中的试验结论一致，这表明在绵羊胚胎发育中后期，体内细胞充分动员转录因子和蛋白分子参与肌纤维发育，积极调动能量运输维持肌肉生成，后期需深度挖掘氧化磷酸化和泛素化途径以全面系统的解释发生机制。

4 结论
本研究首次获得绵羊胚胎骨骼肌发育的全转录组图谱，分析了circRNAs的差异表达，并找到快肌慢肌转换的关键转录物circRNA7527以及与肌细胞增殖分化相关的核心转录因子circRNA17939、circRNA2324和circRNA19073等。靶基因预测和功能分析表明，在D85至D105时主要调控肌管融合和肌细胞的增殖，在D105至D135时主要涉及肌纤维的分化和类型转换，期间相关差异表达RNA主要参与PI3K-Akt、Rap1、Ras、MAPK等信号通路，与上述RNA的靶基因预测和功能聚类分析一致。

References

[1] BENTZINGER C F, XIN W Y, RUDNICKI M A. Building muscle: Molecular regulation of myogenesis. Cold Spring Harbor Perspectives in Biology, 2012, 4(2): a008342.

[2] WEINTRAUB H. The MyoD family and myogenesis: Redundancy, networks, and thresholds. Cell, 1993, 75(7): 1241-1244.

[3] SABOURIN L A, RUDNICKI M A. The molecular regulation of myogenesis. Clinical Genetics, 2010, 57(1): 16-25.

[4] BUCKINGHAM M, RIGBY P W. Gene regulatory networks and transcriptional mechanisms that control myogenesis. Developmental Cell, 2014, 28(3): 225-238.

[5] KAWAKAMI K, SATO S, H, IKEDA K. Six family genes-structure and function as transcription factors and their roles in development. Bioessays, 2000, 22(7): 616-626.

[6] TESSMAR K, LOOSLI F, WITTBRODT J. A screen for co-factors of Six3. Mechanisms of Development, 2002, 117(1): 103-113.

[7] ZHU C C, DYER M A, UCHIKAWA M, KONDOH H, LAGUTIN O V, OLIVER G. Six3-mediated auto repression and eye development requires its interaction with members of the Groucho-related family of co-repressors. Development, 2002, 129(12): 2835-2849.

[8] MCKINNELL I W, JEFF I, FABIEN L G, PUNCH V G J, ADDICKS G C, GREENBLATT J F, F JEFFREY D, RUDNICKI M A. Pax7 activates myogenic genes by recruitment of a histone methyltransferase complex. Nature Cell Biology, 2008, 10(1): 77-84.

[9] 李雪娇, 刘晨曦, 杨开伦, 刘明军. 德美羊与中美羊胎儿期骨骼肌组织学结构发育特征差异性研究. 草食家畜, 2017(4): 1-6.

LI X J, LIU C X,YNGA K L, LIU M J. Study on differentiation of fetal skeletal muscle development characteristics between German and Chinese Merino Sheep. Grass-feeding Livestock, 2017(4): 1-6.(in Chinese)

[10] 张伟. 绵羊骨骼肌高表达miRNA靶基因的高通量获取及部分候选靶基因生物功能研究[D]. 石河子: 石河子大学, 2015.

ZHANG W. Research on acquiring target genes of miRNA expressed at a high level in sheep skeletal muscle by a high throughput way and function of partial target[D]. Shihezi: Shihezi University, 2015. (in Chinese)

[11] 张世芳, 魏彩虹, 陆健, 张小宁, 周鑫磊, 张淑珍, 王光凯, 曹家雪, 赵福平, 张莉, 杜立新 深度测序鉴定绵羊microRNA转录组. 中国畜牧兽医, 2013, 40(9): 19-22.

ZHANG S F, WEI C H, LU J, ZHANG X N, ZHOU X L, ZHANG S Z, WANG G K, CAO J X, ZHAO F P, ZHANG L, DU L X. Identification of microRNAome in texel sheep by deep sequencing. China Animal Husbandry ＆Veterinary Medicine, 2013, 40(9): 19-22. (in Chinese)

[12] 黄万龙, 张秀秀, 李嫒, 苗向阳. 利用RNA-seq技术筛选大白猪皮下和肌内脂肪组织差异表达基因. 遗传, 2017, 39(6): 501-511.

HUANG W L, ZHANG X X, LI Y, MIAO X Y. Identification of differentially expressed genes between subcutaneous and intramuscularadipose tissue of Large White pig using RNA-seq. Hereditas(Beijing), 2017, 39(6): 501-511.(in Chinese)

[13] MARTIN M. Cutadapt removes adapter sequences from high- throughput sequencing reads. EMBnet Journal, 2011, doi: 10.14806/ ej.17.1.200.

[14] LANGMEAD B, SALZBERG S L. Fast gapped-read alignment with Bowtie 2. Nature Methods, 2012, 9: 357-359.

[15] KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R, SALZBERG S L. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology, 2013, 14(4): R36.

[16] MIHAELA P, PERTEA G M, ANTONESCU C M, TSUNG-CHENG C, MENDELL J T, SALZBERG S L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature Biotechnology, 2015, 33(3): 290-295.

[17] FRAZEE A C, GEO P, JAFFE A E, BEN L, SALZBERG S L, LEEK J T. Ballgown bridges the gap between transcriptome assembly and expression analysis. Nature Biotechnology, 2015, 33(3): 243.

[18] ANDERS S, HUBER W. Differential expression analysis for sequence count data. Genome Biology, 2010, 11(10): R106.doi:10.1186/gb- 2010-11-10-r106.

[19] KIM D, SALZBERG S L. TopHat-Fusion: an algorithm for discovery of novel fusion transcripts. Genome Biology, 2011, 12(8): R72.

[20] ZHANG X O, WANG H B, ZHANG Y, LU X, CHEN L L, YANG L. Complementary sequence-mediated exon circularization. Cell, 2014, 159(1): 134-147.

[21] Ashwal-Fluss R, Meyer M, Pamudurti N R, Ivanov A, Bartok O, Hanan M, Evantal N, Memczak S, Rajewsky N, Kadener S. circRNA biogenesis competes with Pre-mRNA splicing. Molecular Cell, 2014, 56(1): 55-66.

[22] CAMON E, MAGRANE M, BARRELL D, LEE V, DIMMER E, MASLEN J, BINNS D, HARTE N, LOPEZ R, APWEILER R. The gene ontology annotation (GOA) database: sharing knowledge in uniprot with gene ontology. Nucleic Acids Research, 2004, 32 (Database issue): D262.

[23] HATTORI M, ITOH M, ARAKI M, HIRAKAWA M, KAWASHIMA S, OKUDA S, GOTO S, KATAYAMA T, TOKIMATSU T, YAMANISHI Y, KANEHISA M. KEGG for linking genomes to life and the environment. Nucleic Acids Research, 2007, 36(suppl.1): D480-D484.

[24] MUKAI H, TOSHIMORI M, SHIBATA H, TAKANAGA H, KITAGAWA M, MIYAHARA M, SHIMAKAWA M, ONO Y. Interaction of PKN with alpha-actinin. Journal of Biological Chemistry, 1997, 272(8): 4740.

[25] MINDEN A, LIN A, MCMAHON M, LANGE-CARTER C, DÉRIJARD B, DAVIS R J, JOHNSON G L, KARIN M. Differential activation of ERK and JNK mitogen-activated protein kinases by Raf-1 and MEKK. Science, 1994, 266(5191): 1719-1723.

[26] LEONARDO S, LAURA P, YVONNE T, LEV K, PIER PAOLO P. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell, 2011, 146(3): 353-358.

[27] 李雪娇, 刘晨曦, 孙亚伟, 杨开伦, 刘明军. 德国美利奴羊胎儿期骨骼肌组织学结构发育特征研究. 西北农林科技大学学报(自然科学版), 2018, 46(5): 7-13.

LI X J, LIU C X, SUN Y W, YANG K L, LIU M J. Study on structure development characteristics of German Merino sheep fetal skeletal muscle tissue. Journal of Northwest A＆F University(Natural Science Edition), 2018, 46(5): 7-13.(in Chinese)

[28] WANG D Z, VALDEZ M R, MCANALLY J, RICHARDSON J, OLSON E N. The Mef2c gene is a direct transcriptional target of myogenic bHLH and MEF2 proteins during skeletal muscle development. Development, 2001, 128(22): 4623.

[29] ZHANG Y Y, VIK T A, RYDER J W, SROUR E F, JACKS T, ., SHANNON K, CLAPP D W. Nf1 regulates hematopoietic progenitor cell growth and ras signaling in response to multiple cytokines. Journal of Experimental Medicine, 1998, 187(11): 1893-1902.

[30] HEGEDUS B, DASGUPTA B, SHIN J E, EMNETT R J, HART-MAHON E K, ELGHAZI L, BERNAL-MIZRACHI E, GUTMANN D H. Neurofibromatosis-1 regulates neuronal and glial cell differentiation from neuroglial progenitors in vivo by both cAMP- and Ras-dependent mechanisms. Cell Stem Cell, 2007, 1(4): 443-457.

[31] XIAOHUA W, ESTWICK S A, SHI C, MENGGANG Y, WENYU M, NEBESIO T D, YAN L, JIN Y, REUBEN K, DAVID I. Neurofibromin plays a critical role in modulating osteoblast differentiation of mesenchymal stem/progenitor cells. Human Molecular Genetics, 2006, 15(19): 2837-2845.

[32] NADINE K, SIGMAR S, CHRISTIAN R D, ROBINSON P N, JOHNNY K, CAROLA D, MONIKA O, JIRKO K, STEVENSON D A, THOMAS B. Neurofibromin (Nf1) is required for skeletal muscle development. Human Molecular Genetics, 2011, 20(14): 2697.

[33] PERRY R L, PARKER M H, RUDNICKI M A. Activated MEK1 binds the nuclear MyoD transcriptional complex to repress transactivation. Molecular Cell, 2001, 8(2): 291-301.

[34] RAMOCKI M B, JOHNSON S E, WHITE M A, ASHENDEL C L, KONIECZNY S F, TAPAROWSKY E J. Signaling through mitogen-activated protein kinase and Rac/Rho does not duplicate the effects of activated Ras on skeletal myogenesis. Molecular & Cellular Biology, 1997, 17(7): 3547-3555.

[35] PAGE J L, WANG X L, JOHNSON S E. MEKK1 signaling through p38 leads to transcriptional inactivation of E47 and repression of skeletal myogenesis. Journal of Biological Chemistry, 2004, 279(30): 30966-30972.

[36] GREDINGER E, GERBER A N, TAMIR Y, TAPSCOTT S J, BENGAL E. Mitogen-activated protein kinase pathway is involved in the differentiation of muscle cells. Journal of Biological Chemistry, 1998, 273(17): 10436-10444.

[37] BOST F, AOUADI M, CARON L, BINéTRUY B. The role of MAPKs in adipocyte differentiation and obesity. Biochimie, 2005, 87(1): 51-56.

[38] MYRIAM A, KATHIANE L, MATTHIEU P, YANNICK M B, BERNARD B, FRéDéRIC B. Inhibition of p38MAPK increases adipogenesis from embryonic to adult stages. Diabetes, 2006, 55(2): 281.

[39] FRéDéRIC B, MYRIAM A, LESLIE C, PATRICK E, NATHALIE B, MATTHIEU P, CHRISTIAN D, PAUL H, GILLES P, JACQUES P. The extracellular signal-regulated kinase isoform ERK1 is specifically required for in vitro and in vivo adipogenesis. Diabetes, 2005, 54(2): 402-411.

[40] LIU S, HAN W, JIANG S, ZHAO C, WU C. Integrative transcriptomics and proteomics analysis of longissimus dorsi muscles of Canadian double-muscled Large White pigs. Gene, 2016, 577(1): 14-23.

[41] VOILLET V, SAN CRISTOBAL M, PèRE M-C, BILLON Y, CANARIO L, LIAUBET L, LEFAUCHEUR L. Integrated analysis of proteomic and transcriptomic data highlights late fetal muscle maturation process. Molecular & Cellular Proteomics:MCP, 2018, 17(4): 672-693.

[42] OUYANG H, WANG Z, CHEN X, YU J, LI Z, NIE Q. Proteomic analysis of chicken skeletal muscle during embryonic development. Frontiers in physiology, 2017, 8: 281-317.

Analysis and Identification of circRNAs of Skeletal Muscle at Different Stages of Sheep Embryos Based on Whole Transcriptome Sequencing

SHI TianPei, WANG XinYue, HOU HaoBin, ZHAO ZhiDa, SHANG MingYu, ZHANG Li

(Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193)

Abstract: 【Objective】The meat production of livestock, which is closely related to the development of skeletal muscle, is an important economic trait to measure the quality of livestock. For mammals, the skeletal muscle development depends on the growth and differentiation of embryonic myocyte, which has a significant impact on the subsequent growing potential. In this study, the developmental mode of skeletal muscle, the important transformation nodes, the formation of muscle fibers and the molecular regulation mechanism of transformation were mainly explored. 【Method】Based on the previous research, the important nodes D85, D105 and D135 related to the myotube development were used in the experiment, and the longissimus dorsi muscles were sequenced by whole transcriptome sequencing. The differentially expressed (DE) circRNAs were screened by bioinformatics analysis and verified by quantitative real-time PCR (qRT-PCR). 【Result】1 126 DE circRNAs were obtained by conditional screening (|log2| ≥1 andp≤0.05). The 3 groups were compared and many specific expressions of circRNA were found at each stage, but in the D85 vs D135 group, the amount was the most. 374 DE circRNAs were obtained, which contained 201 up-regulated and 173 down-regulated, and 44.7% of the DE genes were differentially expressed with a difference of more than 4 times. These DE circRNAs were subjected to run GO and KEGG functional analysis and targeted prediction, and they were enriched into some pathways, such as energy metabolism and signal transduction, which involved in muscle differentiation and muscle fiber development, including MAPK, PI3K-Akt, Ras, regulation of actin cytoskeleton and other signal transduction pathways. According to the results, it was confirmed that the DE circRNAs enriched during D85 to D105 were mostly associated with cell proliferation and survival, regulation of myocyte development and cell cycle, while D105 to D135 were mainly related to energy conversion, material transport, RNA transport, and DNA repair. By drawing co-expression visualization network with the targeted prediction results used by Cytoscape, the core regulatory transcripts, such as circRNA8239, circRNA19073, circRNA2765 and circRNA1616, were identified. In the D105 period, a key factor circRNA7527 that regulated the conversion of fast and slow muscle types was found, which targets the bta-miR-135a, bta-miR-615, and chi-miR-133a-5p to regulate the MEF2C gene. According to the differential expression and functional prediction in three comparison groups, 4 circRNAs related to muscle development and 4 target miRNA were selected for qRT-PCR, and the results showed that the gene expression trend was consistent with the sequencing data. 【Conclusion】It was verified that the stabilization of the number of muscle fibers occurred between sheep embryos at D85 and D105, and muscle fiber hypertrophy happened during the D105 to D135 period, which lead to the conclusion that D105 was probably a key time point. In this study, we firstly constructed a circRNA map in sheep embryonic skeletal muscle development based on the whole transcriptome sequencing. The transcriptome differences at key stages were revealed, and multiple circRNAs and miRNAs targeting MEF2Cthat involved in the MAPK signaling pathway were found, which provided reference for livestock myofiber development research and other research on non-coding RNA.

Key words: sheep; embryo; whole transcriptome; skeletal muscle; growth and development; circRNA

收稿日期：2019-03-01；

接受日期：2019-05-30

基金项目：国家自然科学基金联合基金重点支持项目NSFC（U1503285）、 中央级公益性科研院所基本科研业务费专项（Y2017XM02）

联系方式：石田培，E-mail：yangstp@foxmail.com。通信作者张莉，E-mail：zhangli07@caas.cn

开放科学（资源服务）标识码（OSID）：width=42.5,height=42.5

（责任编辑 林鉴非）