石田培，王欣悅，侯浩賓，趙志達(dá)，尚明玉，張莉

基于全轉(zhuǎn)錄組測(cè)序的綿羊胚胎不同發(fā)育階段骨骼肌circRNA的分析與鑒定

石田培，王欣悅，侯浩賓，趙志達(dá)，尚明玉，張莉

（中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193）

【目的】肉用性狀是家養(yǎng)動(dòng)物的重要性狀，尤其是產(chǎn)肉性能與骨骼肌的發(fā)育密切相關(guān)，對(duì)多數(shù)動(dòng)物而言，骨骼肌的生長取決于肌纖維的早期發(fā)育。本研究將揭示綿羊（）胚胎期骨骼肌組織發(fā)育重要節(jié)點(diǎn)、肌纖維的形成及轉(zhuǎn)換機(jī)制，探究骨骼肌在動(dòng)物胚胎期的生長對(duì)后續(xù)發(fā)育潛力的影響?！痉椒ā吭谇捌诰d羊胚胎期組織結(jié)構(gòu)試驗(yàn)的基礎(chǔ)上，選擇肌纖維發(fā)育及類型轉(zhuǎn)換的3個(gè)重要時(shí)間節(jié)點(diǎn)：胎齡85 d（D85）、105 d（D105）、135 d（D135），并對(duì)處于這些節(jié)點(diǎn)期的中國美利奴綿羊胚胎背最長肌進(jìn)行全轉(zhuǎn)錄組測(cè)序。通過生物信息學(xué)分析和功能預(yù)測(cè)篩選出差異表達(dá)顯著的circRNA，采用實(shí)時(shí)熒光定量 PCR（qRT-PCR）對(duì)其進(jìn)行驗(yàn)證?！窘Y(jié)果】通過條件篩選（|log2| ≥1且≤0.05）獲得差異表達(dá)circRNA 1 126個(gè)。將3組進(jìn)行比較，在各時(shí)間節(jié)點(diǎn)都發(fā)現(xiàn)較多的特異性表達(dá)circRNA。在D85 vs D135比較組中，特異性差異表達(dá)數(shù)量最多：共發(fā)現(xiàn)差異表達(dá)circRNA 374個(gè)，上調(diào)表達(dá)201個(gè)，下調(diào)表達(dá)173個(gè)，其中具有4倍以上差異的基因167個(gè)，占差異基因的44.7%。對(duì)這些差異性表達(dá)的circRNA進(jìn)行GO和KEGG功能分析和靶向預(yù)測(cè)，富集到與肌肉分化和肌纖維發(fā)育相關(guān)的能量代謝和信號(hào)轉(zhuǎn)導(dǎo)通路，包括MAPK、PI3K-Akt、Ras、Regulation of actin cytoskeleton等。結(jié)果表明：在D85至D105期間富集到的差異circRNA多與肌細(xì)胞發(fā)育調(diào)控、細(xì)胞增殖和存活、細(xì)胞周期等相關(guān)，而在D105至D135期間富集到的通路則主要與能量轉(zhuǎn)換、物質(zhì)運(yùn)輸、RNA轉(zhuǎn)運(yùn)和DNA修復(fù)有關(guān)。靶向預(yù)測(cè)結(jié)果通過Cytoscape軟件繪制出可視化共表達(dá)網(wǎng)絡(luò)，篩選到circRNA8239、circRNA19073，circRNA2765、circRNA1616等核心調(diào)控轉(zhuǎn)錄物，并在D105找到調(diào)節(jié)快慢肌類型轉(zhuǎn)換的關(guān)鍵因子circRNA7527，其靶向作用于bta-miR-135a、bta-miR-615、chi-miR-133a-5p進(jìn)而調(diào)控MEF2C基因。通過3個(gè)時(shí)間節(jié)點(diǎn)的差異表達(dá)情況和功能預(yù)測(cè)描述選取了與肌肉發(fā)育相關(guān)的4個(gè)核心circRNA和其靶向的4個(gè)miRNA進(jìn)行qRT-PCR分析，其基因表達(dá)趨勢(shì)與測(cè)序數(shù)據(jù)一致?！窘Y(jié)論】本試驗(yàn)在轉(zhuǎn)錄水平上驗(yàn)證了綿羊胚胎發(fā)育D85至D105是肌纖維數(shù)量的穩(wěn)定發(fā)生期，而D105至D135則是肌纖維的肥大期，由此推斷D105可能是綿羊胚胎發(fā)育的關(guān)鍵時(shí)間窗口。本研究首次基于全轉(zhuǎn)錄組測(cè)序構(gòu)建了綿羊胚胎期骨骼肌發(fā)育的circRNA圖譜，揭示了不同發(fā)育階段的轉(zhuǎn)錄組差異，并對(duì)綿羊胚胎骨骼肌發(fā)育期間重要的調(diào)控因子進(jìn)行挖掘和驗(yàn)證，發(fā)現(xiàn)了以為靶基因的多個(gè)circRNA和miRNA參與MAPK信號(hào)通路，為綿羊肌纖維發(fā)育分子機(jī)制以及其它家畜非編碼RNA的研究提供了參考。

綿羊；胚胎；全轉(zhuǎn)錄組；骨骼??；生長發(fā)育；circRNA

0 引言

【研究意義】肉用性狀是家畜的重要經(jīng)濟(jì)性狀，研究骨骼肌發(fā)生機(jī)制有助于提升家畜產(chǎn)肉性能及肌肉品質(zhì)。骨骼肌的形成根據(jù)發(fā)育特點(diǎn)可以分為兩個(gè)階段：成肌細(xì)胞的增殖、遷移和分化，以及肌纖維的延伸、增粗和類型轉(zhuǎn)換。研究表明成肌細(xì)胞在動(dòng)物妊娠中后期融合形成梭形排列的肌管，最終形成終末分化的肌原纖維[1]。家畜肌纖維數(shù)量在妊娠后期達(dá)到一定數(shù)目后就不再發(fā)生變化，而成年后肌肉量的增加主要依靠肌纖維長度和周徑的增加。由此可見，胚胎期骨骼肌的生長和發(fā)育對(duì)家畜的肉用性狀至關(guān)重要?！厩叭搜芯窟M(jìn)展】目前，已有大量關(guān)于骨骼肌轉(zhuǎn)錄組和蛋白質(zhì)組發(fā)生機(jī)制的研究，一些蛋白質(zhì)和轉(zhuǎn)錄物作為調(diào)節(jié)因子在胚胎肌細(xì)胞發(fā)育過程中發(fā)揮著重要作用，如肌源性調(diào)節(jié)因子（myogenic regulatory factors, MRFs）通過肌源性因子5（myogenic factor 5, Myf5）、肌源性分化1（myogenic differentiation 1, MyoD）、肌源性調(diào)節(jié)因子4（myogenic regulatory factor 4, Mrf4）和肌細(xì)胞生成素（Myogenin）等協(xié)同調(diào)控胚胎骨骼肌的發(fā)育和分化[2-3]，其中每個(gè)MRF都可單獨(dú)作為肌生成的主要調(diào)節(jié)因子，精細(xì)調(diào)節(jié)胚胎骨骼肌的發(fā)育過程[4]。Six家族蛋白質(zhì)作為轉(zhuǎn)錄因子存在兩個(gè)特征性結(jié)構(gòu)域：一個(gè)是可與DNA結(jié)合的Six型同源結(jié)構(gòu)域，另一個(gè)是與轉(zhuǎn)錄激活或抑制因子相互作用的氨基末端結(jié)構(gòu)域。Six1（The sine oculis–related homeobox 1）和Six4通常被認(rèn)為是遺傳調(diào)節(jié)級(jí)聯(lián)的制高點(diǎn)[5-6]，能夠誘導(dǎo)肌源性祖細(xì)胞發(fā)育成為肌源性細(xì)胞系。除經(jīng)典調(diào)控因子外，Pax7能與Wdr5-Ash2L-MLL2組蛋白甲基轉(zhuǎn)移酶（HMT）復(fù)合物結(jié)合，在Myf5的作用下，致使染色質(zhì)周圍的H3K4三甲基化[7]。因此，通過Pax因子富集HMT復(fù)合物的途徑被普遍認(rèn)為是可能重塑染色質(zhì)結(jié)構(gòu)以控制譜系特異性基因表達(dá)的保守機(jī)制[8]。近年來，有關(guān)骨骼肌發(fā)育非編碼RNA的研究較多[9-11]，但受樣品、設(shè)備及分析方法等因素的限制，在全轉(zhuǎn)錄組水平上同時(shí)對(duì)多種非編碼RNA的關(guān)聯(lián)分析還未見報(bào)道。【本研究切入點(diǎn)】肌肉形成與發(fā)生機(jī)制是一個(gè)復(fù)雜的動(dòng)態(tài)過程，目前對(duì)該生理過程還缺乏系統(tǒng)的研究。本研究將使用二代測(cè)序技術(shù)，從基因轉(zhuǎn)錄組層面對(duì)肌肉發(fā)生機(jī)制進(jìn)行系統(tǒng)的、全面的特征化和量化分析[12]?！緮M解決的關(guān)鍵問題】應(yīng)用RNA sequencing技術(shù)對(duì)綿羊不同發(fā)育時(shí)期的胚胎骨骼肌進(jìn)行全轉(zhuǎn)錄組測(cè)序和生物信息學(xué)分析，對(duì)骨骼肌早期發(fā)育過程中的關(guān)鍵調(diào)控機(jī)制進(jìn)行探究，發(fā)掘重要的功能性調(diào)控非編碼RNA。

1 材料與方法

1.1 動(dòng)物道德倫理聲明

本研究方法嚴(yán)格按照中華人民共和國農(nóng)業(yè)農(nóng)村部制定的相關(guān)指導(dǎo)方針和規(guī)定進(jìn)行，所有試驗(yàn)方案均經(jīng)中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所批準(zhǔn)。

1.2 實(shí)驗(yàn)動(dòng)物和樣品制備

試驗(yàn)選擇年齡相近、體況良好、體重在 55—60 kg的成年中國美利奴羊（購自新疆巴州種畜場(chǎng)），經(jīng)埋栓同期發(fā)情處理后進(jìn)行人工授精配種，通過超數(shù)排卵獲得胚胎并移植到體況相近的受體母羊，妊娠母羊均按照國家農(nóng)業(yè)行業(yè)肉羊飼養(yǎng)標(biāo)準(zhǔn)（NY/T816- 2004）進(jìn)行飼喂。胚胎分別在妊娠第85天、105天和135天時(shí)通過外科手術(shù)取出（每組采樣個(gè)體為 3 只），采集背最長肌，并分成小塊于-80℃冰箱保存?zhèn)溆谩?/p>

1.3 方法

1.3.1 RNA 的提取及質(zhì)控 準(zhǔn)備足量樣品，按照Trizol試劑盒提取樣品總RNA，提取后使用微量核酸蛋白測(cè)定儀（scandrop100）檢測(cè)OD 260 nm /OD 280 nm 值以鑒定RNA樣品濃度，排除RNA污染可能性。Bioanalyzer 2100（Agilent, Santa Clara, CA）用于評(píng)估總 RNA 質(zhì)量，以 RIN≥7 且260/280≥1.8為閾值，RNase-free DNase I (Ambion Inc., Texas, USA)用于消除潛在的基因組 DNA 污染。分離的RNA 樣本保存在-80℃環(huán)境下。

1.3.2 cDNA 文庫構(gòu)建和 RNA 測(cè)序 RNA質(zhì)檢合格后，分別從3個(gè)階段肌肉樣品中各取出約10 μg RNA，分別構(gòu)建測(cè)序文庫，主要包括酶促RNA 片段化、cDNA 合成、測(cè)序接頭銜接及 PCR 擴(kuò)增等過程。首先，根據(jù)EpicentreRibo-Zero Gold試劑盒（Illumina，San Diego，USA）的操作說明，消耗樣品中的核糖體RNA、線性RNA并純化后，在高溫下使其片段化，隨后轉(zhuǎn)錄成cDNA文庫，最后分別利用Illumina Hiseq 4000進(jìn)行測(cè)序。

1.3.3 數(shù)據(jù)處理及轉(zhuǎn)錄水平分析 對(duì)測(cè)序產(chǎn)生的原始數(shù)據(jù)進(jìn)行Cutadapt[13]過濾，除去低質(zhì)量、含N的序列后通過Bowtie2[14]和Tophat2[15]將讀數(shù)映射到綿羊的3.1基因組，其中獲得的注釋文件包含轉(zhuǎn)錄本的分布信息。使用StringTie[16]識(shí)別新的轉(zhuǎn)錄本，并在此基礎(chǔ)上利用Ballgown[17]估計(jì)轉(zhuǎn)錄物的表達(dá)水平。然后通過FPKM算法消除偏好影響，進(jìn)行基因表達(dá)定量。在本研究中，采用edgeR包確定差異表達(dá)基因[18]，同時(shí)使用上四分位數(shù)算法校正數(shù)據(jù)，即設(shè)定| log2（兩個(gè)樣本的FPKM比率）| ≥1且≤0.05為閾值。TopHat- fusion[19]和CIRCExplorer2[20]可以在未映射序列中識(shí)別反向剪接讀數(shù)來估計(jì)circRNA的表達(dá)水平[21]，circRNA及其相應(yīng)親本基因的共表達(dá)關(guān)系通過Pearson相關(guān)系數(shù)2值估算確定。miRNA序列的處理需使用專有腳本，同時(shí)使用多個(gè)過濾器除去不可映射的測(cè)序讀數(shù)。通過比對(duì)miRBase 21.0中pre-miRNA（mir）和成熟miRNA（miR）的獨(dú)特序列或已公開數(shù)據(jù)庫進(jìn)行“映射”，剩余序列進(jìn)行BLAST。最后通過RNA折疊軟件從側(cè)翼80nt序列處預(yù)測(cè)含有發(fā)夾RNA結(jié)構(gòu)的序列。本研究通過TMeV軟件對(duì)差異表達(dá)的circRNA進(jìn)行熱圖繪制，直觀展示各樣本的表達(dá)一般規(guī)律、差別及變化，從而通過聚類預(yù)測(cè)新型circRNA的功能及作用。CircRNA-miRNA互作分析主要采用Targetscan與miRanda軟件進(jìn)行靶標(biāo)預(yù)測(cè)后取交集，將得到的數(shù)據(jù)上傳至Cytoscape繪制可視化互作網(wǎng)絡(luò)圖。

1.3.4 差異表達(dá)基因 GO 及 KEGG 通路富集分析 經(jīng)定量確定的差異表達(dá)RNA仍需深入了解其生物學(xué)功能以發(fā)掘重要的調(diào)控因子，在生物信息學(xué)分析中，常采用基因本體論（gene ontology，GO）[22]進(jìn)行分析。在GO體系中，可以依次從細(xì)胞成分、生物學(xué)過程和分子功能3個(gè)方面定位基因功能，使用其注釋富集功能將具有相似注釋的轉(zhuǎn)錄物整合歸類。在預(yù)測(cè)各類未知RNA時(shí)，除了通過RNA的宿主基因進(jìn)行功能預(yù)測(cè)外，后期還可間接通過它們?cè)谏镄盘?hào)通路中的參與情況進(jìn)行功能分析。生物通路富集分析主要依賴于京都基因和基因組百科全書（kyoto encyclopedia of genes and genomes , KEGG）[23]數(shù)據(jù)庫，其中途徑富集分析可以精確定位差異表達(dá)基因的主要代謝轉(zhuǎn)導(dǎo)和信號(hào)傳導(dǎo)途徑。本研究使用Scripts in house軟件進(jìn)行GO（http://geneontology.org）和KEGG（http://www.kegg. jp/kegg）富集分析，設(shè)置統(tǒng)計(jì)≤0.05為結(jié)果顯著閾值。

1.3.5 實(shí)時(shí)熒光定量 PCR 驗(yàn)證 為確保測(cè)序結(jié)果的準(zhǔn)確性和有效性，需要通過qRT-PCR 驗(yàn)證基因表達(dá)水平。選擇驗(yàn)證的RNA至少在一組對(duì)照組中表達(dá)差異顯著，且與骨骼肌的生長發(fā)育相關(guān)。經(jīng)過組間對(duì)比和篩選，最終分別選擇circRNA和miRNA各4個(gè)進(jìn)行驗(yàn)證。根據(jù)轉(zhuǎn)錄本數(shù)據(jù)庫設(shè)計(jì)轉(zhuǎn)錄物引物（表1），使用CFX-Connect TM實(shí)時(shí)系統(tǒng)（BioRad）上的SYBR green I master mix（Roche，GmBH，Basel，Switzerland）進(jìn)行實(shí)時(shí)PCR。試驗(yàn)結(jié)果用目的基因值減去內(nèi)參基因值得到ΔCt，再以相對(duì)參照樣品的ΔCt均值為對(duì)照組，用各組ΔCt值減去對(duì)照組得到ΔΔCt，以2-ΔΔCt平均值表達(dá)出各個(gè)基因的相對(duì)表達(dá)量，利用SPSS 22.0軟件對(duì)所得表達(dá)量進(jìn)行單因素方差分析(One-way ANOVA)，結(jié)果分為差異顯著（≤0.05）和差異極顯著（≤0.01）。

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

2 結(jié)果

2.1 RNA-seq 數(shù)據(jù)分析

采用RNA-Seq和Small RNA-Seq深度測(cè)序，每個(gè)樣本文庫內(nèi)獲得約90 000 000的原始數(shù)據(jù)，質(zhì)量控制后平均得到約 88 000 000的有效讀段。Q-score＞20的reads 在 99%以上，超過97%的reads Q-score＞30。在有效讀段中，約85% reads比對(duì)到綿羊參考基因組上，其中超過82%是唯一比對(duì)，除此之外還進(jìn)行了非剪接/剪接比對(duì)計(jì)數(shù)和正/反義鏈計(jì)數(shù)。綜上，充分證明樣本數(shù)據(jù)處理得當(dāng)，質(zhì)量可靠，符合下游數(shù)據(jù)分析要求。在深度分析樣本注釋結(jié)果，評(píng)估骨骼肌中轉(zhuǎn)錄本豐度及表達(dá)差異情況前，計(jì)算評(píng)估了各生物學(xué)重復(fù)間皮爾森積矩相關(guān)系數(shù)平方（2）均大于 0.98，均滿足統(tǒng)計(jì)學(xué)要求。

2.2 差異表達(dá)轉(zhuǎn)錄本鑒定

在3個(gè)對(duì)照組中，共鑒定37 474條circRNA，利用 FPKM 算法標(biāo)準(zhǔn)化基因轉(zhuǎn)錄水平差異，篩選出顯著差異circRNA 1 126條。在轉(zhuǎn)錄本識(shí)別過程中，circRNA表達(dá)量FPKM值結(jié)果顯示轉(zhuǎn)錄本表達(dá)水平存在明顯差異，但仍有相當(dāng)一部分沒有表達(dá)或者表達(dá)量極低（圖1）。對(duì)3組數(shù)據(jù)進(jìn)行比對(duì)，circRNA的表達(dá)情況表明各階段都存在特異性表達(dá)的非編碼RNA，但在D85 vs D135的發(fā)育階段中，特異性表達(dá)數(shù)量最多。兩組比對(duì)出374個(gè)差異表達(dá)基因，其中上調(diào)表達(dá)201個(gè)，下調(diào)表達(dá)173個(gè)，具有4倍以上差異的基因167個(gè)，約占差異基因的44.7%（圖2）。miRNA識(shí)別元件（miRNA response element, MRE）的存在可介導(dǎo)miRNA與靶標(biāo)轉(zhuǎn)錄物的部分互補(bǔ)序列結(jié)合，導(dǎo)致機(jī)體相應(yīng)靶標(biāo)基因的表達(dá)抑制。因此，通過TargetScan和miRanda預(yù)測(cè)miRNA和mRNA之間的靶向作用，共找到76 759個(gè)互作關(guān)系，根據(jù)GO term和KEGG pathway功能描述，篩選出與骨骼肌發(fā)育相關(guān)的互作關(guān)系497對(duì)。在D85 vs D105 vs D135的比較組中，將基因按照GO匹配顯著度排名，其中、和極為顯著。

2.3 差異表達(dá)circRNA的功能富集分析

基于GO數(shù)據(jù)庫的注釋結(jié)果表明，Unigenes中共有20 840條轉(zhuǎn)錄本得到歸類注釋，顯著富集于353個(gè)生物學(xué)過程條目，117個(gè)分子功能條目，144個(gè)細(xì)胞組分條目，主要富集于纖維組裝、細(xì)胞增殖與遷移、細(xì)胞分裂、骨骼肌收縮等生物學(xué)過程。3個(gè)階段胚胎背最長肌circRNA顯著富集于124個(gè)生物學(xué)過程條目，36個(gè)分子功能條目和40個(gè)細(xì)胞組分條目。3個(gè)發(fā)育階段都富集到與骨骼肌發(fā)育相關(guān)的生物學(xué)過程，但顯著高表達(dá)基因注釋結(jié)果略有差異（圖3）。在D85 vs D105中，信號(hào)轉(zhuǎn)導(dǎo)調(diào)節(jié)、轉(zhuǎn)錄活性調(diào)節(jié)和骨骼肌細(xì)胞分化相關(guān)的基因富集最為明顯；在D105 vs D135中，肌纖維發(fā)育、肌細(xì)胞增殖和胚胎發(fā)育調(diào)節(jié)相關(guān)的基因富集最為明顯；在D85 vs D135中，顯著富集并與其它階段重合的通路有細(xì)胞增殖與轉(zhuǎn)運(yùn)，肌細(xì)胞分化和信號(hào)轉(zhuǎn)導(dǎo)等。通過KEGG 通路富集分析，發(fā)現(xiàn)差異表達(dá)基因在多個(gè)信號(hào)轉(zhuǎn)導(dǎo)通路中被顯著富集，其中富集于心肌病、Wnt信號(hào)通路和Ras信號(hào)通路的基因比例較大。結(jié)合RNA的GO分析和 KEGG分析結(jié)果發(fā)現(xiàn)，綿羊胚胎骨骼肌的發(fā)育同時(shí)受到多種機(jī)制的調(diào)控，在不同階段涉及的信號(hào)通路略有差別。差異表達(dá)基因被顯著富集于69個(gè)通路，其中涵蓋了多個(gè)基因。如圖4所示，MAPK、PI3K-Akt、Ras、Regulation of actin cytoskeleton等與肌肉發(fā)育相關(guān)能量代謝及轉(zhuǎn)運(yùn)通路顯著富集。在所有被富集途徑中，PI3K-Akt信號(hào)傳導(dǎo)途徑是參與細(xì)胞凋亡、蛋白質(zhì)合成和細(xì)胞周期等細(xì)胞基本功能的較活躍途徑。在信號(hào)系統(tǒng)中常通過信使作用產(chǎn)生級(jí)聯(lián)反應(yīng)，首先生長因子與受體酪氨酸激酶（RTK）或G蛋白偶聯(lián)受體（GPCR）結(jié)合分別刺激Ia和Ib類活化，繼而催化細(xì)胞膜上磷脂酰肌醇-3,4,5-三磷酸（PIP3）的產(chǎn)生。然后作為第二個(gè)信使激活，便可以通過磷酸化調(diào)控關(guān)鍵的細(xì)胞過程?；罨a(chǎn)生的激活了、等第二信使導(dǎo)致、[24]、等在細(xì)胞增殖過程中上調(diào)表達(dá)，[25]、等下調(diào)表達(dá)，穩(wěn)定維持肌細(xì)胞的增殖合成，血管生成以及DNA損傷修復(fù)。在繪制的信號(hào)轉(zhuǎn)導(dǎo)通路中清晰可見，大部分與細(xì)胞發(fā)育增殖、蛋白合成和DNA損傷修復(fù)相關(guān)的基因表達(dá)顯著，呈現(xiàn)橙色（圖4）。此外，通過對(duì)D85 vs D105、D105 vs D135和D85 vs D135等3組差異表達(dá)circRNA進(jìn)行聚類篩選，最終挑選出37個(gè)具有骨骼肌細(xì)胞分化、骨骼肌纖維發(fā)育和骨骼肌細(xì)絲組裝等功能描述的circRNA，根據(jù)轉(zhuǎn)錄水平差異對(duì)其繪制聚類熱圖（圖5），由表達(dá)模式聚類可預(yù)測(cè)未能得到注釋的轉(zhuǎn)錄本。在熱圖中，這些circRNA都至少在一個(gè)時(shí)間節(jié)點(diǎn)顯著上調(diào)或下調(diào)。

圖1 不同樣本circRNA的FPKM值分布

圖2 不同樣本circRNA表達(dá)統(tǒng)計(jì)圖

A：D85 vs D105；B：D105 vs D135；C：D85 vs D135；Y軸：pathway名稱；X軸：富集因子；氣泡大小：富集到某個(gè)通路中差異表達(dá)的數(shù)量/背景中所有數(shù)量；顏色反映了富集在某通路差異表達(dá)數(shù)及富集顯著性

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 差異表達(dá)circRNA的KEGG通路分析

Fig. 3 KEGG pathway analysis of differentially expressed circRNAs

2.4 差異表達(dá)miRNA-circRNA互作網(wǎng)絡(luò)

在轉(zhuǎn)錄調(diào)控中，網(wǎng)絡(luò)調(diào)控分析可以鑒定出潛在調(diào)控重要性狀的關(guān)鍵候選基因[26]。有研究表明circRNA可作為參與miRNA表達(dá)調(diào)控的競(jìng)爭(zhēng)性內(nèi)源RNA（competing endogenous RNA, ceRNA），本研究采用Targetscan和miRanda軟件來分析circRNA-miRNA相互作用關(guān)系。Targetscan基于種子區(qū)域進(jìn)行miRNA目標(biāo)預(yù)測(cè)，而miRanda基于circRNA和miRNA的自由能的組合，根據(jù)軟件評(píng)分結(jié)果，確定了877個(gè)顯著差異circRNA，它們均與多個(gè)miRNA之間存在互作，共鑒定到6 340個(gè)miRNA-circRNA相互作用對(duì)，其中194個(gè)相互作用對(duì)與骨骼肌發(fā)育相關(guān)，我們將這些與骨骼肌發(fā)育相關(guān)的miRNA和circRNA構(gòu)建了相互作用網(wǎng)絡(luò)。在網(wǎng)絡(luò)圖中，某一節(jié)點(diǎn)與另一節(jié)點(diǎn)的連線稱之為度，度的大小可直觀反映基因在網(wǎng)絡(luò)中的重要程度，擁有相對(duì)較大度值的基因可能具有更重要的生物學(xué)價(jià)值。計(jì)算共表達(dá)網(wǎng)絡(luò)中所有基因的度，根據(jù)共表達(dá)網(wǎng)絡(luò)中所有基因度值排名，提取網(wǎng)絡(luò)中連接度較大的circRNA、miRNA，對(duì)其表達(dá)情況進(jìn)行驗(yàn)證。circRNA5502、circRNA5505、circRNA5561、circRNA30828和circRNA23584在互作網(wǎng)絡(luò)中度值較大，表明其宿主基因可能在骨骼肌細(xì)胞生長代謝及發(fā)育過程中起重要作用（圖6）。

白色小框是根據(jù)已有的數(shù)據(jù)繪制的具有一般參考意義的，概括詳盡的代謝圖；綠色小框?yàn)樵撐锓N特有的基因或酶，具有更詳細(xì)的信息；橙色方框表示該基因在此時(shí)差異性表達(dá)

2.5 構(gòu)建circRNA-miRNA-mRNA共表達(dá)調(diào)控網(wǎng)絡(luò)

將篩選得到的與肌肉發(fā)育及肌細(xì)胞增殖分化相關(guān)的差異表達(dá)circRNA及其靶向的miRNA和mRNA進(jìn)行轉(zhuǎn)錄調(diào)控分析和共表達(dá)網(wǎng)絡(luò)分析。在互作網(wǎng)絡(luò)中，4個(gè)特異性差異表達(dá)circRNA和其互作的miRNA與mRNA共形成208個(gè)互作關(guān)系，根據(jù)度值發(fā)現(xiàn)、、等在網(wǎng)絡(luò)中具有重要調(diào)控地位（圖7）。表達(dá)的circRNA1616受chi-miR- 30f-3p反式調(diào)控作用于和，表達(dá)的circRNA19073，circRNA2765作用于、、、等。結(jié)合三者互作的GO分析結(jié)果（圖8），發(fā)現(xiàn)這些表達(dá)性基因主要被富集在胞質(zhì)、核質(zhì)和ATP結(jié)合等能量代謝相關(guān)的生物學(xué)過程，KEGG通路中顯著富集到的通路都是參與代謝和細(xì)胞生長的途徑，如Ubiquitin mediated proteolysis、Wnt signaling pathway、Insulin signaling pathway、Focal adhesion、Cell cycle。根據(jù)以上結(jié)果可推測(cè)，在骨骼肌發(fā)育的D85至D135過程中，骨骼肌發(fā)生重點(diǎn)進(jìn)行了轉(zhuǎn)移，circRNA后期主要參與能量代謝和物質(zhì)運(yùn)輸?shù)恼{(diào)控。

橫軸：胚胎日齡，縱軸：circRNA名稱；藍(lán)色：下調(diào)，紅色：上調(diào)

2.6 實(shí)時(shí)熒光定量檢測(cè)

實(shí)時(shí)熒光定量結(jié)果表明，獲得的所有基因的溶解曲線均為單峰，說明定量試驗(yàn)設(shè)計(jì)的引物較優(yōu)，并且在定量檢測(cè)過程中不存在非特異性擴(kuò)增與引物二聚體（圖9）。通過對(duì)8條RNA的qPCR檢測(cè)結(jié)果進(jìn)行校正后，檢測(cè)差異表達(dá)水平與深度測(cè)序結(jié)果一致，說明高通量測(cè)序數(shù)據(jù)可靠。繼續(xù)對(duì)功能聚類后篩選出的關(guān)鍵RNA進(jìn)行分析，以D85 vs D105、D105 vs D135和D85 vs D135進(jìn)行組間顯著性檢驗(yàn)。設(shè)定≤0.05為顯著，≤0.01代表極顯著，差異具有統(tǒng)計(jì)學(xué)意義。結(jié)果顯示：circRNA17939、chi-miR-133a-5p和oar-miR- 150在組間差異極顯著，circRNA17939在D85 vs D105階段差異極顯著，bta-miR-365-3p在D85 vs D105階段差異顯著，在D105 vs D135差異極顯著，符合數(shù)據(jù)鑒定和功能分析的結(jié)果。

3 討論

本項(xiàng)目前期動(dòng)物試驗(yàn)表明：在美利奴綿羊胚胎期，從D75至D135，各體尺指標(biāo)（體重、體長和體高）均隨日齡增加而增長。在D75骨骼肌可見中空管狀的初級(jí)肌纖維典型結(jié)構(gòu)，且在初級(jí)肌纖維四周圍繞著大量的次級(jí)肌纖維；D105 和D135切片中只有次級(jí)肌纖維結(jié)構(gòu)而無初級(jí)肌纖維結(jié)構(gòu)特征?？梢娫贒75后，初級(jí)肌纖維與次級(jí)肌纖維發(fā)生融合，肌纖維的形成基本完成；在D105時(shí)肌纖維數(shù)量基本恒定；到D135時(shí)肌纖維的發(fā)育均表現(xiàn)為肌纖維面積不同程度的增大，肌纖維直徑顯著增加[9]。李雪嬌等在綿羊胎兒骨骼肌組織學(xué)結(jié)構(gòu)發(fā)育特征研究中表明胎兒D105是肌纖維從肌增生到肌肥大的關(guān)鍵時(shí)間[27]。基于此，對(duì)綿羊胚胎D75至D105之間進(jìn)行了進(jìn)一步細(xì)化分析，結(jié)果發(fā)現(xiàn)D75至D85肌纖維開始相互融合，D85時(shí)已無顯著特征的初級(jí)肌纖維，形成成熟的肌纖維結(jié)構(gòu)，肌纖維數(shù)量達(dá)到整個(gè)胎兒發(fā)育階段的峰值。因此選擇D85、D105 和D135時(shí)間節(jié)點(diǎn)構(gòu)建綿羊骨骼肌發(fā)育階段的全轉(zhuǎn)錄組文庫，通過 RNA-seq 測(cè)序和生信分析篩選出在轉(zhuǎn)錄水平上可能參與調(diào)控的RNA，結(jié)果共得到顯著差異表達(dá)circRNA1 126 個(gè)。對(duì)篩選結(jié)果進(jìn)行GO和KEGG富集分析，發(fā)現(xiàn)3個(gè)階段都顯著富集到大量與肌肉發(fā)育相關(guān)的基因和通路，如RNA信號(hào)轉(zhuǎn)導(dǎo)、物質(zhì)代謝和轉(zhuǎn)錄活性通路等。GO分析結(jié)果顯示，circRNA及其互作的miRNA和mRNA在骨骼肌的生長發(fā)育中都具有重要功能，主要包括肌細(xì)胞發(fā)育（GO:0055001）、骨骼肌衛(wèi)星細(xì)胞分化（GO:0014816）、肌纖維發(fā)育（GO:0048747）、骨骼肌收縮（GO:0003009）、骨骼肌纖維發(fā)育（GO:0006979）和骨骼肌組織生長的積極調(diào)節(jié)（GO:0048633）等生物學(xué)過程。在D85 vs D105比較組中，鑒定到circRNA23653、circRNA19073和circRNA6142等與肌細(xì)胞增殖和調(diào)控相關(guān)的circRNA，還發(fā)現(xiàn)circRNA7527與慢速和快速纖維之間的轉(zhuǎn)換密切相關(guān)。在D85 vs D135比較組中，發(fā)現(xiàn)與肌纖維發(fā)育相關(guān)的circRNA15個(gè)：circRNA1615、circRNA5174、circRNA3179、circRNA15129、circRNA7311、circRNA1616、circRNA6291、circRNA2324、circRNA18785、circRNA19839、circRNA1602、circRNA19702、circRNA31369、circRNA7569、circRNA489。比較分析發(fā)現(xiàn)circRNA7527未在D85 vs D135差異表達(dá)，說明circRNA7527只在D105時(shí)期特異性表達(dá)，推測(cè)其功能可能是調(diào)控肌纖維的轉(zhuǎn)換。根據(jù)轉(zhuǎn)錄水平的對(duì)比，篩選到、、、、、、、、、、、等15個(gè)重要基因，涉及到肌纖維發(fā)育的信號(hào)通路包括TGFβ、MAPK、PI3K- Akt、Ras、Akt/mTOR等。circRNA18785、circRNA19839、circRNA1602、circRNA19702、circRNA31369、circRNA7569、circRNA489。比較分析發(fā)現(xiàn)circRNA7527未在D85 vs D135差異表達(dá)，說明circRNA7527只在D105時(shí)期特異性表達(dá)，推測(cè)其功能可能是調(diào)控肌纖維的轉(zhuǎn)換。根據(jù)轉(zhuǎn)錄水平的對(duì)比，篩選到、、、、、、、、、、、等15個(gè)重要基因，涉及到肌纖維發(fā)育的信號(hào)通路包括TGFβ、MAPK、PI3K- Akt、Ras、Akt/mTOR等。

橙色橢圓表示差異表達(dá)的miRNA；綠色箭頭表示差異表達(dá)的circRNA；邊表示兩者之間的互相作用關(guān)系。顏色由淺到深代表表達(dá)差異程度，顏色越深差異越顯著，反之則顯著程度降低

Y軸：生物學(xué)過程；X軸：富集因子；氣泡大小反映了富集到某個(gè)GO_Term的差異表達(dá)的數(shù)量/背景中所有數(shù)量；顏色反映了富集在某過程差異表達(dá)數(shù)及富集顯著性

成熟的miRNA通過特定序列互補(bǔ)配對(duì)識(shí)別靶mRNA，并根據(jù)互補(bǔ)程度沉默復(fù)合體降解mRNA或阻抑mRNA的翻譯。與miRNA單一的作用機(jī)制不同，circRNA可從轉(zhuǎn)錄、轉(zhuǎn)錄后及表觀遺傳學(xué)等多個(gè)水平參與基因的調(diào)控。circRNA可以包含核糖體進(jìn)入位點(diǎn)，可以翻譯表達(dá)有效蛋白，轉(zhuǎn)錄時(shí)與mRNA前體發(fā)生剪切競(jìng)爭(zhēng)，還可與基因上游啟動(dòng)子序列結(jié)合調(diào)控下游基因表達(dá)，作為相關(guān)分子元件與特定蛋白結(jié)合調(diào)節(jié)相應(yīng)蛋白活性，因此circRNA的功能研究更為復(fù)雜。本研究中發(fā)現(xiàn)多個(gè)circRNA通過靶向調(diào)控或吸附miRNA而順式調(diào)控其宿主基因的表達(dá)。在D105時(shí)期找到關(guān)鍵性的circRNA7527，并在bta-miR-135a、bta-miR-615、chi-miR-133a-5p和mmu-mir-6240-p5_5序列中發(fā)現(xiàn)可與其結(jié)合的位點(diǎn)，此外，根據(jù)功能性靶向預(yù)測(cè)，發(fā)現(xiàn)這些miRNA都可作用于。在D105 至D135時(shí)期，bta-miR-135a、bta-miR-615和mmu-mir-6240- p5_5隨著circRNA的下調(diào)而上調(diào)，但chi-miR-133a-5p卻隨之下調(diào)，說明bta-miR-135a、bta-miR-615和mmu- mir-6240-p5_5受到的調(diào)控機(jī)制可能與chi-miR-133a- 5p不完全相同。除此之外，預(yù)測(cè)出的miRNA（bta-miR- 135a、bta-miR-615和mmu-mir-6240-p5_5）的靶基因上下調(diào)水平呈現(xiàn)相反趨勢(shì)，該表達(dá)規(guī)律符合哺乳動(dòng)物體內(nèi)miRNA與mRNA的作用機(jī)制。但chi-miR-133a-5p與的下調(diào)情況相同，說明其之間可能不存在靶向關(guān)系。有研究證明在胚胎發(fā)生過程中心肌和骨骼肌組織分化開始時(shí)，在家族其他基因開始表達(dá)后隨即表達(dá)，但單個(gè)的功能缺失會(huì)造成肌細(xì)胞的分化[28]。由此可見circRNA7527、bta-miR-135a、bta-miR-615、chi-miR- 133a-5p和mmu-mir-6240-p5_5中的一個(gè)或幾個(gè)對(duì)肌纖維的分化具有重要調(diào)控意義。在比較組中還發(fā)現(xiàn)circRNA2324可與bta-miR-219作用靶向調(diào)節(jié)，但該circRNA只在D85和D105兩個(gè)階段上調(diào)表達(dá)，在D135時(shí)顯著下調(diào)，說明在胚胎骨骼肌發(fā)育后期，circRNA2324的作用受到抑制。[29-30]編碼神經(jīng)纖維素蛋白，具有特異性GTP酶激活結(jié)構(gòu)域，可作為/絲裂原活化蛋白激酶（）信號(hào)傳導(dǎo)的負(fù)調(diào)節(jié)物。由于胞內(nèi)表達(dá)可誘導(dǎo)細(xì)胞退出終止細(xì)胞周期，導(dǎo)致肌細(xì)胞終末分化[31]。所以早期間充質(zhì)中的失活可導(dǎo)致肌纖維數(shù)量驟減，機(jī)體全部肌肉纖維化[32]。此外，的缺失會(huì)造成的減少，而是早期負(fù)向調(diào)節(jié)肌肉分化的重要因子，與核內(nèi)復(fù)合物相互作用并抑制其功能[33]。上游的和/或激酶的組成型活化可以抑制肌細(xì)胞分化[34-35]，的下游MAPK信號(hào)傳導(dǎo)參與肌原性分化的調(diào)節(jié)[36]，這與之前進(jìn)行的功能性pathway分析結(jié)果相符。

圖中為RNA-seq結(jié)果和qRT-PCR的定量結(jié)果；折線圖表示測(cè)序的FPKM值，柱狀圖表示相對(duì)定量結(jié)果；顯著性檢驗(yàn)：D105柱狀圖上標(biāo)符號(hào)為D85 vs D105比較結(jié)果；D135柱狀圖上2個(gè)標(biāo)符號(hào)，左邊為D85 vs D135比較結(jié)果，右邊為D105 vs D135比較結(jié)果；* P≤0.05，** P≤0.01

KEGG富集分析發(fā)現(xiàn)大多數(shù)差異表達(dá)circRNA主要富集在與肌細(xì)胞發(fā)育周期和分化的生理過程中，介導(dǎo)了細(xì)胞對(duì)外界的反應(yīng)。、、和5/是MAPK通路中重要亞族[37-38]，有研究證明它們是調(diào)控肌肉細(xì)胞和脂肪細(xì)胞增殖和分化的重要候選基因[39]。在繪制的信號(hào)通路圖中，發(fā)現(xiàn)由激活的第二信使調(diào)控的、、等基因在細(xì)胞增殖、存活和肌動(dòng)蛋白的組裝過程中起著重要作用。根據(jù)差異表達(dá)情況可知，這3個(gè)基因在該通路中的受調(diào)控作用顯著（圖4）。常常通過與骨骼肌和非骨骼肌型α-輔肌動(dòng)蛋白的第三血影蛋白重復(fù)序列結(jié)合，并以Ca2 +不敏感的方式與骨骼肌型α-輔肌動(dòng)蛋白的的EF-手狀基序的區(qū)域結(jié)合而直接關(guān)聯(lián)與細(xì)胞骨架網(wǎng)絡(luò)的連接，在肌動(dòng)蛋白骨架組裝中具有重要作用。蛋白通過在GDP結(jié)合狀態(tài)和GTP結(jié)合狀態(tài)間循環(huán)轉(zhuǎn)換來轉(zhuǎn)導(dǎo)胞外生長因子的信號(hào)，活化的（-GTP）主要可激活兩種蛋白激酶，即和激酶（）。激活絲裂原活化蛋白激酶（），包括細(xì)胞外信號(hào)調(diào)節(jié)激酶（）和激酶（），活化的直接促進(jìn)活化但不促進(jìn)活化，而參與活化，但僅在過表達(dá)后引起ERK活化，即兩種不同的依賴性級(jí)聯(lián)：一個(gè)由引發(fā)導(dǎo)致的活化，另一個(gè)是由引發(fā)導(dǎo)致的活化。在富集到的PI3K-Akt，Rap1和Ras信號(hào)通路中，顯著體現(xiàn)了引發(fā)的活化，促使基因表達(dá)的級(jí)聯(lián)模式。以上結(jié)果均表明MAPK等信號(hào)通路是參與調(diào)控肌細(xì)胞增殖分化的關(guān)鍵通路，直接或間接地參與綿羊胚胎肌肉的發(fā)育調(diào)控。

值得一提的是，LIU等對(duì)初生大白豬（）的肌肉進(jìn)行了轉(zhuǎn)錄組和蛋白組共分析，發(fā)現(xiàn)肌肉中差異表達(dá)基因主要被富集在分泌糖蛋白的Wnt、MAPK、TGF-β和 cAMP信號(hào)通路中，并驗(yàn)證了可通過促進(jìn)成纖維細(xì)胞增殖和神經(jīng)細(xì)胞分化而影響肌肉發(fā)育及出生后肌肉肥大[40]，這與本研究中被富集到的顯著通路和纖維生長關(guān)鍵基因基本相同。VALENTIN等對(duì)不同胚胎期的大白豬骨骼肌發(fā)育機(jī)制進(jìn)行了多組學(xué)關(guān)聯(lián)分析，結(jié)果發(fā)現(xiàn)骨骼肌發(fā)生后期至肌肉成熟階段，轉(zhuǎn)錄物和蛋白質(zhì)的功能主要集中在物質(zhì)能量代謝（氧化磷酸化等）中[41]。此外，HONG等在家禽（）中也有相似的胚胎階段性研究，雞胚的E11、E16及D1階段蛋白質(zhì)互作網(wǎng)絡(luò)分析表明所有差異表達(dá)蛋白質(zhì)主要參與蛋白質(zhì)合成、能量代謝、肌肉收縮和氧化磷酸化過程[42]。本研究結(jié)果與上述豬和家禽中的試驗(yàn)結(jié)論一致，這表明在綿羊胚胎發(fā)育中后期，體內(nèi)細(xì)胞充分動(dòng)員轉(zhuǎn)錄因子和蛋白分子參與肌纖維發(fā)育，積極調(diào)動(dòng)能量運(yùn)輸維持肌肉生成，后期需深度挖掘氧化磷酸化和泛素化途徑以全面系統(tǒng)的解釋發(fā)生機(jī)制。

4 結(jié)論

本研究首次獲得綿羊胚胎骨骼肌發(fā)育的全轉(zhuǎn)錄組圖譜，分析了circRNAs的差異表達(dá)，并找到快肌慢肌轉(zhuǎn)換的關(guān)鍵轉(zhuǎn)錄物circRNA7527以及與肌細(xì)胞增殖分化相關(guān)的核心轉(zhuǎn)錄因子circRNA17939、circRNA2324和circRNA19073等。靶基因預(yù)測(cè)和功能分析表明，在D85至D105時(shí)主要調(diào)控肌管融合和肌細(xì)胞的增殖，在D105至D135時(shí)主要涉及肌纖維的分化和類型轉(zhuǎn)換，期間相關(guān)差異表達(dá)RNA主要參與PI3K-Akt、Rap1、Ras、MAPK等信號(hào)通路，與上述RNA的靶基因預(yù)測(cè)和功能聚類分析一致。

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

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

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

[4] BUCKINGHAM M, RIGBY P W. Gene regulatory networks and transcriptional mechanisms that control myogenesis., 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., 2000, 22(7): 616-626.

[6] TESSMAR K, LOOSLI F, WITTBRODT J. A screen for co-factors of Six3., 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., 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., 2008, 10(1): 77-84.

[9] 李雪嬌, 劉晨曦, 楊開倫, 劉明軍. 德美羊與中美羊胎兒期骨骼肌組織學(xué)結(jié)構(gòu)發(fā)育特征差異性研究. 草食家畜, 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., 2017(4): 1-6.(in Chinese)

[10] 張偉. 綿羊骨骼肌高表達(dá)miRNA靶基因的高通量獲取及部分候選靶基因生物功能研究[D]. 石河子: 石河子大學(xué), 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] 張世芳, 魏彩虹, 陸健, 張小寧, 周鑫磊, 張淑珍, 王光凱, 曹家雪, 趙福平, 張莉, 杜立新 深度測(cè)序鑒定綿羊microRNA轉(zhuǎn)錄組. 中國畜牧獸醫(yī), 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.2013, 40(9): 19-22. (in Chinese)

[12] 黃萬龍, 張秀秀, 李嬡, 苗向陽. 利用RNA-seq技術(shù)篩選大白豬皮下和肌內(nèi)脂肪組織差異表達(dá)基因. 遺傳, 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., 2017, 39(6): 501-511.(in Chinese)

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

[14] LANGMEAD B, SALZBERG S L. Fast gapped-read alignment with Bowtie 2., 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., 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., 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., 2015, 33(3): 243.

[18] ANDERS S, HUBER W. Differential expression analysis for sequence count data., 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.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., 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., 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., 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., 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., 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., 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?, 2011, 146(3): 353-358.

[27] 李雪嬌, 劉晨曦, 孫亞偉, 楊開倫, 劉明軍. 德國美利奴羊胎兒期骨骼肌組織學(xué)結(jié)構(gòu)發(fā)育特征研究. 西北農(nóng)林科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 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., 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., 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., 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 progenitorsby both cAMP- and Ras-dependent mechanisms., 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., 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., 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., 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., 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., 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., 1998, 273(17): 10436-10444.

[37] BOST F, AOUADI M, CARON L, BINéTRUY B. The role of MAPKs in adipocyte differentiation and obesity., 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., 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., 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., 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., 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., 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)

【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 and0.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 thegene. 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 targetingthat involved in the MAPK signaling pathway were found, which provided reference for livestock myofiber development research and other research on non-coding RNA.

sheep; embryo; whole transcriptome; skeletal muscle; growth and development; circRNA

2019-03-01；

2019-05-30

國家自然科學(xué)基金聯(lián)合基金重點(diǎn)支持項(xiàng)目NSFC（U1503285）、 中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)專項(xiàng)（Y2017XM02）

石田培，E-mail：yangstp@foxmail.com。通信作者張莉，E-mail：zhangli07@caas.cn

（責(zé)任編輯 林鑒非）