鴨胚胎發(fā)育中后期胸肌發(fā)育阻滯的RNA-seq分析

劉宏祥，徐文娟，朱春紅，陶志云，宋衛(wèi)濤，章雙杰，李慧芳

?

鴨胚胎發(fā)育中后期胸肌發(fā)育阻滯的RNA-seq分析

劉宏祥，徐文娟，朱春紅，陶志云，宋衛(wèi)濤，章雙杰，李慧芳

（江蘇省家禽科學(xué)研究所，江蘇揚(yáng)州 225125）

【目的】選擇兩個中國地方品種高郵鴨和金定鴨，對鴨胚胎發(fā)育中后期胸肌進(jìn)行轉(zhuǎn)錄組分析研究，旨在探明胸肌發(fā)育阻滯的分子變化機(jī)制，為鴨骨骼肌調(diào)控機(jī)理研究打下基礎(chǔ)?！痉椒ā吭?1胚齡和27胚齡兩個時間點(diǎn)，分別解剖高郵鴨、金定鴨各3只，采集胸大肌，提取總RNA構(gòu)建文庫，利用Illumina的HiseqTM2000進(jìn)行高通量測序，并利用生物信息學(xué)方法進(jìn)行差異表達(dá)基因挖掘、基因功能注釋等分析，探討21胚齡和27胚齡兩個時間點(diǎn)之間胸肌發(fā)育阻滯的分子機(jī)制。【結(jié)果】高郵鴨、金定鴨21胚齡和27胚齡胸大肌組織RNA-seq質(zhì)量Q20均在94%以上，Q30均在89%以上，測序得到的結(jié)果可靠，可用于后續(xù)分析。RNA水平相關(guān)性檢查和基因mRNA表達(dá)量聚類圖結(jié)果都表明，21胚齡（27胚齡）高郵鴨和金定鴨之間表達(dá)模式的相關(guān)性高于高郵鴨（金定鴨）21胚齡和27胚齡之間的表達(dá)模式。不同品種內(nèi)時間點(diǎn)之間的差異基因數(shù)量（高郵鴨6 128個，金定鴨6 452個）遠(yuǎn)多于同一時間點(diǎn)不同品種間的差異基因數(shù)量（21胚齡522個，27胚齡299個）。qRT-PCR驗(yàn)證試驗(yàn)結(jié)果與RNA-seq分析結(jié)果相關(guān)性較強(qiáng)。通過GO和KEGG富集分析發(fā)現(xiàn)，高郵鴨、金定鴨胸肌在21胚齡到27胚齡階段，能量代謝相關(guān)基因（主要為輔酶Q相關(guān)基因、ATP酶合成相關(guān)基因和細(xì)胞色素C相關(guān)基因）均顯著上調(diào)，DNA復(fù)制和細(xì)胞周期相關(guān)基因（主要為微型染色體維持蛋白（MCM）相關(guān)基因、復(fù)制因子C（RFC）相關(guān)基因）均顯著下調(diào)。相關(guān)基因表達(dá)的變化可能與此階段成肌細(xì)胞增殖速度減慢，逐漸退出細(xì)胞周期開始準(zhǔn)備下一階段融合成多核肌管并形成肌纖維有關(guān)。對肌肉生長發(fā)育相關(guān)的關(guān)鍵基因分析發(fā)現(xiàn)，促進(jìn)肌肉生長的1和誘導(dǎo)成肌細(xì)胞末端分化的顯著下調(diào)，促進(jìn)肌纖維分化融合的1基因、誘導(dǎo)肌祖細(xì)胞向成肌細(xì)胞轉(zhuǎn)化的1基因顯著上調(diào)?！窘Y(jié)論】鴨胚胎中后期胸肌發(fā)育過程中大量基因差異表達(dá)。其中能量代謝相關(guān)基因的上調(diào)和DNA復(fù)制、細(xì)胞周期相關(guān)基因的下調(diào)，以及肌肉發(fā)育相關(guān)基因1顯著上調(diào)，1、等顯著下調(diào)，可能與鴨胚胎中后期胸肌發(fā)育阻滯現(xiàn)象密切相關(guān)。

鴨；胸肌發(fā)育；轉(zhuǎn)錄組測序；差異表達(dá)基因

0 引言

【研究意義】動物肌肉的生長發(fā)育大體可以分為兩個階段：出雛（生）前的胚胎期階段和出雛（生）后階段。動物肌纖維數(shù)量在出雛（生）前就已經(jīng)固定[1-3]，后期肌肉量的增加主要由于肌纖維的肥大，因此禽類上市日齡時的產(chǎn)肉量取決于胚胎中后期肌纖維的數(shù)量。胸肌是禽類肌肉的重要組成部分，因此研究禽類胚胎期胸肌發(fā)育的調(diào)控機(jī)制，挖掘影響肌纖維發(fā)育的關(guān)鍵基因，對育種實(shí)踐中禽類產(chǎn)肉量的選育具有重要指導(dǎo)意義?！厩叭搜芯窟M(jìn)展】動物在胚胎發(fā)育期間，成肌細(xì)胞增殖、分化成為多核的肌管，最終形成成熟的肌纖維，此過程涉及肌纖維的增生和肥大[2]。哺乳動物胚胎期胸肌重量持續(xù)增長[4]，而禽類胚胎期胸肌發(fā)育過程較為復(fù)雜，表現(xiàn)為胚胎前期發(fā)育迅速而中后期發(fā)育遲緩甚至萎縮。MOORE等[5]報(bào)道火雞胸肌肌纖維橫截面積在胚胎后期減小。CHEN等[6-7]對鴨的研究發(fā)現(xiàn)，鴨胸肌重量和肌纖維橫截面積從22胚齡到出雛期呈減小趨勢。本課題組前期研究發(fā)現(xiàn)高郵鴨和金定鴨胚胎期21胚齡后胸肌重量沒有增加[8]。這可能與此時胸肌肌纖維主要表現(xiàn)為融合，而非肥大有關(guān)[7]，但具體分子機(jī)制還未完全了解。近年來，隨著鴨基因組測序的完成[9]（基因組草圖見http://www.ensembl.org/Anas_platyrhynchos/Info/ Index），以及高通量測序、分析技術(shù)的發(fā)展，RNA-seq分析鴨不同發(fā)育時期基因mRNA表達(dá)水平變化情況及其調(diào)控機(jī)制成為可能。【本研究切入點(diǎn)】鴨胚胎中后期胸肌發(fā)育阻滯現(xiàn)象已有多方報(bào)道，但具體分子機(jī)制還未完全明晰?！緮M解決的關(guān)鍵問題】高郵鴨和金定鴨分別為中國地方體型差異較大的品種，其肌肉生長具有明顯的表型差異，這為研究肌肉生長發(fā)育的分子機(jī)制提供了良好的素材。本研究擬通過RNA-seq技術(shù)方法，探討胚胎發(fā)育中后期高郵鴨、金定鴨的胸肌轉(zhuǎn)錄組表達(dá)差異情況，為深入了解此階段胸肌生長發(fā)育的分子機(jī)制奠定基礎(chǔ)。

1 材料與方法

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

本試驗(yàn)研究對象為中國地方品種鴨——蛋肉兼用型品種高郵鴨和蛋用型品種金定鴨，于2014年飼養(yǎng)于江蘇省高郵鴨集團(tuán)。選擇正常的高郵鴨和金定鴨種蛋（重量接近）各10枚于微電腦全自動孵化器（山東德州，科裕孵化機(jī)）中孵化，21胚齡和27胚齡時隨機(jī)挑選高郵鴨、金定鴨胚蛋各3枚。解剖胚胎，采集右側(cè)胸大肌樣品，迅速置于液氮中凍存。

1.2 試驗(yàn)方法

1.2.1 總RNA提取與質(zhì)量鑒定 使用Trizol試劑（Invitrogen）從約30 mg胸大肌樣品中提取總RNA，使用1%甲醛凝膠對RNA樣品進(jìn)行電泳，鑒定RNA完整性和降解情況；使用Qubit RNA試劑盒（Invitrogen）測定RNA濃度；使用Nano Photometer分光光度計(jì)（德國Implen）檢測RNA純度。

1.2.2 測序文庫的構(gòu)建 RNA檢測合格后，將mRNA打斷成短片段并反轉(zhuǎn)錄成cDNA，選擇合適大小的片段進(jìn)行PCR富集得到最終的cDNA文庫。使用qRT-PCR方法對文庫的有效濃度進(jìn)行準(zhǔn)確定量（文庫有效濃度＞2 nM），以保證文庫質(zhì)量。

1.2.3 Illumina測序 庫檢合格后，交由北京諾禾致源公司使用HiseqTM2000測序儀（Illumina）進(jìn)行轉(zhuǎn)錄組測序。

1.2.4 qRT-PCR驗(yàn)證 根據(jù)轉(zhuǎn)錄組測序結(jié)果，挑選16個差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證。使用Primer 3在線工具設(shè)計(jì)引物（http://fokker.wi.mit.edu/primer3/ input.htm）（所選基因和設(shè)計(jì)的引物見表1）。

表1 用于qRT-PCR驗(yàn)證的基因及其引物

*看家基因Housekeeping genes

1.3 測序數(shù)據(jù)的處理分析

1.3.1 測序數(shù)據(jù)的處理 為保證信息分析質(zhì)量，對原始序列（Raw reads）進(jìn)行過濾去雜（去除測序接頭序列、重復(fù)冗余序列、低質(zhì)量序列）得到高質(zhì)量的序列數(shù)據(jù)（clean reads）。采用TopHat2方法[10]將clean reads與鴨參考基因組進(jìn)行比對。

1.3.2 差異基因的篩選 本試驗(yàn)每個胚齡點(diǎn)都有3個生物學(xué)重復(fù)。通過DESeq軟件包，對前步得到的clean reads進(jìn)行標(biāo)準(zhǔn)化，并評估read count的分布，進(jìn)行BH檢驗(yàn)（Benjamini和Hochberg于1995年提出的控制假陽性結(jié)果錯誤率的多重比較方法[11]），最終得到不同組別之間的差異基因比較結(jié)果。

1.3.3 mRNA差異表達(dá)基因GO分類和KEGG富集分析 使用R語言中的goseq包[12]將mRNA差異表達(dá)基因序列與GO（gene ontology）數(shù)據(jù)庫進(jìn)行比對分析，獲得GO功能注釋；將mRNA差異表達(dá)基因與KEGG（kyoto encyclopedia of genes and genomes）數(shù)據(jù)庫進(jìn)行BLASTX比對獲得mRNA差異表達(dá)基因相對應(yīng)的Pathway注釋信息。

2 結(jié)果

2.1 RNA-seq數(shù)據(jù)總體評價

本研究對高郵鴨、金定鴨21胚齡和27胚齡胸大肌組織轉(zhuǎn)錄組測序，高郵鴨和金定鴨兩個時間點(diǎn)均獲得了43+ M條的reads，總堿基數(shù)5+ Gb，Q20接近95%，Q30接近90%（表2），說明RNA-seq測序結(jié)果可靠，可用于后續(xù)分析。高郵鴨、金定鴨21胚齡和27胚齡樣品比對到參考基因組上的reads均達(dá)到總reads數(shù)的60%左右（表3），比對率均較高。

表2 測序數(shù)據(jù)評估統(tǒng)計(jì)

Q20、Q30分別指質(zhì)量值大于等于20或30的堿基所占百分比。GC為測序結(jié)果中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水平相關(guān)性檢查

使用皮爾遜相關(guān)系數(shù)對樣品進(jìn)行相關(guān)性分析。同一時間點(diǎn)、品種內(nèi)的3個生物學(xué)重復(fù)的相關(guān)系數(shù)平方（2）均大于0.95，同一時間點(diǎn)不同品種之間的相似性高于同一品種不同時間點(diǎn)之間的相似性（圖1）。

2.3 基因mRNA表達(dá)聚類圖

根據(jù)高郵鴨、金定鴨21胚齡、27胚齡各基因RPKM進(jìn)行層次聚類分析。同一時間點(diǎn)不同品種之間的表達(dá)模式相近，而同一品種不同時間點(diǎn)的表達(dá)模式差異較大（圖2）。

2.4 mRNA差異表達(dá)基因的篩選

表4顯示了21胚齡和27胚齡樣品的RPKM值。21胚齡高郵鴨、金定鴨胸大肌組織高表達(dá)基因較多，RPKM值大于3的基因均占到60%左右。

表4 不同胚齡胸大肌RPKM值

圖1 樣品間RNA-seq水平相關(guān)性檢查

圖2 基因mRNA差異表達(dá)聚類圖

高郵鴨和金定鴨不同時間點(diǎn)之間的mRNA差異表達(dá)基因均較多。高郵鴨21胚齡和27胚齡之間、金定鴨21胚齡和27胚齡之間的mRNA顯著差異表達(dá)基因數(shù)量均超過6 000個。21胚齡金定鴨和高郵鴨mRNA顯著差異表達(dá)基因、27胚齡金定鴨和高郵鴨mRNA顯著差異表達(dá)基因分別為522個和299個（圖3），明顯少于品種內(nèi)不同時間點(diǎn)之間的mRNA差異表達(dá)基因數(shù)量。

縱坐標(biāo)1.3以上表示padj<0.05 Numbers higher than 1.3 in Y-axis indicate padj＜0.05

2.5 qRT-PCR驗(yàn)證

從差異表達(dá)分析結(jié)果中挑選出16個基因進(jìn)行qRT-PCR驗(yàn)證。結(jié)果顯示，所選擇的16個基因qRT-PCR結(jié)果（??Ct值）與RNA-seq結(jié)果（RPKM比值）具有較強(qiáng)的相關(guān)性（> 0.70）（表5），表明RNA-seq結(jié)果可靠。

2.6 骨骼肌相關(guān)基因差異表達(dá)情況

挑選出與骨骼肌生長發(fā)育相關(guān)的關(guān)鍵基因1、1、1、和，發(fā)現(xiàn)在21胚齡和27胚齡兩個時間點(diǎn)，高郵鴨和金定鴨之間1、1、、1 mRNA表達(dá)水平均沒有顯著變化。兩個品種21胚齡到27胚齡1、均顯著下調(diào)，1、1均顯著上調(diào)，均沒有顯著變化（表6）。

2.7 mRNA差異表達(dá)基因GO富集分析

分別對高郵鴨、金定鴨21胚齡和27胚齡胸大肌組織的上調(diào)差異基因、下調(diào)差異基因進(jìn)行GO（Gene Ontology）分析，高郵鴨和金定鴨分別富集到24個（上調(diào)和下調(diào)分別2個和22個）和26個顯著GO條目（上調(diào)和下調(diào)分別1個和25個（圖4）。

B1：細(xì)胞成分運(yùn)動（cellular component movement），B2：基于微管的過程（microtubule-based process），B3：基于微管的運(yùn)動（microtubule-based movement），B4：細(xì)胞周期（cell cycle），B5：分解代謝過程的正向調(diào)控（positive regulation of catabolic process），B6：自噬調(diào)控（regulation of autophagy），B7：自噬正調(diào)控（positive regulation of autophagy），B8：細(xì)胞膜組織（cellular membrane organization），B9：（evasion or tolerance of host defenses by virus），B10：細(xì)胞分解過程的正向調(diào)控（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：其他共生互作生物的防御響應(yīng)（response to defenses of other organism involved in symbiotic interaction），B16：對宿主防御的響應(yīng)（response to host defenses），B17：對宿主的響應(yīng)（response to host）

C1：線粒體（mitochondrion），C2：核糖體（ribosome），C3：核糖核蛋白復(fù)合物（ribonucleoprotein complex），C4：細(xì)胞外基質(zhì)（extracellular matrix），C5：大分子復(fù)合物（macromolecular complex），C6：非膜結(jié)合細(xì)胞器（non-membrane-bounded organelle），C7：胞內(nèi)非膜結(jié)合細(xì)胞器（intracellular non-membrane-bounded organelle）

M1：核糖體結(jié)構(gòu)成分（structural constituent of ribosome），M2：肌動活性（motor activity），M3：微管肌動活性（microtubule motor activity），M4：結(jié)構(gòu)分子活性（structural molecule activity），M5：微管結(jié)合（microtubule binding），M6：細(xì)胞支架蛋白結(jié)合（cytoskeletal protein binding），M7：DNA依賴的ATP酶活性（DNA-dependent ATPase activity），M8：微管蛋白結(jié)合（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差異表達(dá)基因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）為上調(diào)基因所富集，與線粒體能量代謝有關(guān)；其余均由下調(diào)基因所富集，主要與細(xì)胞周期、DNA復(fù)制有關(guān)。每個共有GO條目下的共有基因數(shù)量見表7。

2.8 mRNA差異表達(dá)基因KEGG富集分析

本試驗(yàn)分別對高郵鴨、金定鴨21胚齡和27胚齡胸大肌組織的上調(diào)差異基因、下調(diào)差異基因進(jìn)行KEGG分析，高郵鴨和金定鴨分別富集到8個和7個顯著KEGG通路（圖5）。高郵鴨上調(diào)基因富集到4個顯著KEGG通路，下調(diào)基因富集到4個顯著KEGG通路；金定鴨上調(diào)基因富集到4個顯著KEGG通路，下調(diào)基因富集到3個顯著KEGG通路。

表5 qRT-PCR實(shí)驗(yàn)得到的基因表達(dá)??Ct值與RPKM比值之間的相關(guān)

表6 骨骼肌相關(guān)基因mRNA表達(dá)水平變化

表7 高郵鴨、金定鴨21胚齡和27胚齡之間顯著富集的共有GO分類

Term2：檸檬酸循環(huán)（三羧酸循環(huán)）（Citrate cycle (TCA cycle)），Term6：氧化磷酸化（Oxidative phosphorylation），Term7：纈氨酸、亮氨酸和異亮氨酸降解（Valine, leucine and isoleucine degradation），Term9：碳新陳代謝（Carbon metabolism），Term11：核糖體（Ribosome），Term12：DNA復(fù)制（DNA replication），Term14：PPAR信號途徑（PPAR signaling pathway），Term15：細(xì)胞周期（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）由上調(diào)基因所富集，主要與能量代謝有關(guān)；Term12（apla03030，DNA replication）、Term15（apla04110，Cell cycle）由下調(diào)基因所富集，主要與DNA復(fù)制和細(xì)胞周期有關(guān)。每個共有KEGG通路下的共有基因數(shù)量見表8。

表8 高郵鴨、金定鴨21胚齡和27胚齡之間顯著富集的共有KEGG通路

2.9 GO、KEGG聯(lián)合分析

對共有GO條目的基因和共有KEGG通路的基因取交集，篩選出16個上調(diào)基因和15個下調(diào)基因（表9）。上調(diào)基因大多為輔酶Q相關(guān)基因、ATP酶合成相關(guān)基因、細(xì)胞色素C相關(guān)基因，下調(diào)基因大多為微型染色體維持蛋白（MCM）相關(guān)基因、復(fù)制因子C（RFC）相關(guān)基因。

3 討論

肌肉發(fā)育是一個持續(xù)不斷累積的過程，本課題組前期對高郵鴨、金定鴨的研究[8,13]以及其他課題組的多個試驗(yàn)均報(bào)道雞[14]、鴨[6-7]胚胎發(fā)育中后期，胸肌生長遲緩，甚至出現(xiàn)萎縮現(xiàn)象?；痣u胚胎后期胸肌橫截面積減小、肌肉萎縮的同時，肌肉衛(wèi)星細(xì)胞的有絲分裂活性也有降低[5]。為了解這一現(xiàn)象背后的分子機(jī)制，本試驗(yàn)以高郵鴨、金定鴨21胚齡、27胚齡胚胎胸肌為試驗(yàn)素材，通過RNA-seq比較品種之間以及21胚齡、27胚齡兩個時間點(diǎn)之間的mRNA差異表達(dá)基因。挑選部分顯著差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證，結(jié)果表明RNA-seq結(jié)果與熒光定量PCR結(jié)果相似，高通量測序具有較高的可信度。

RNA-seq水平相關(guān)性檢查發(fā)現(xiàn)，同一時間點(diǎn)、品種內(nèi)的3個生物學(xué)重復(fù)之間的表達(dá)模式相似度較高，能夠保證后續(xù)差異基因的可靠分析。同一時間點(diǎn)不同品種之間的相關(guān)性高于同一品種不同時間點(diǎn)之間的相關(guān)性。高郵鴨21胚齡、27胚齡的差異基因和金定鴨21胚齡、27胚齡的差異基因數(shù)量遠(yuǎn)多于21胚齡金定鴨、高郵鴨之間的差異基因數(shù)量和27胚齡金定鴨、高郵鴨之間的差異基因數(shù)量。對高郵鴨、金定鴨兩個時間點(diǎn)的基因按照其RPKM值進(jìn)行層次聚類分析，21胚齡高郵鴨、金定鴨聚為一類，27胚齡高郵鴨、金定鴨聚為一類。以上結(jié)果均表明，在21胚齡到27胚齡的發(fā)育過程中，高郵鴨與金定鴨的胸肌表達(dá)模式具有相似性，在胸肌發(fā)育過程中時間因素的影響要大于品種之間遺傳因素的影響。

在GO和KEGG富集分析中，上調(diào)基因富集到的GO條目、KEGG通路主要與線粒體能量代謝有關(guān)；下調(diào)基因富集到的GO條目、KEGG通路主要與細(xì)胞周期、DNA復(fù)制有關(guān)。篩選出GO條目和KEGG通路中的共有基因，發(fā)現(xiàn)與能量代謝有關(guān)的上調(diào)基因主要為輔酶Q相關(guān)基因、ATP酶合成相關(guān)基因以及細(xì)胞色素c相關(guān)基因，與DNA復(fù)制和細(xì)胞周期有關(guān)的下調(diào)基因主要為微型染色體維持蛋白（MCM）相關(guān)基因、復(fù)制因子C（RFC）相關(guān)基因。輔酶Q對體內(nèi)呼吸鏈中的質(zhì)子移位及電子傳遞起重要作用，是細(xì)胞代謝和細(xì)胞呼吸的激活劑[15-16]。細(xì)胞色素c也是呼吸鏈中可流動的遞氫體或遞電子體。上調(diào)基因中，輔酶Q相關(guān)基因（NDUF*等）分別參與了NADH還原酶1α和1β等復(fù)合物的編碼合成；ATP合成酶相關(guān)基因（ATP5F1、ATP5H和ATP5J）參與了ATP合成酶復(fù)合物的編碼合成。21胚齡到27胚齡階段，鴨胸肌組織中這些基因的上調(diào)，增強(qiáng)了能量的轉(zhuǎn)化水平和能量代謝頻率。真核細(xì)胞DNA的復(fù)制和有絲分裂需要微型染色體維持蛋白（MCM，Minichromosome maintenance protein complex）[17-19]和復(fù)制因子C[20-22]的參與。21胚齡到27胚齡鴨胸肌DNA復(fù)制和細(xì)胞周期相關(guān)基因下調(diào)，說明在此階段胸肌細(xì)胞的增殖速度開始減緩。細(xì)胞周期阻滯是后續(xù)細(xì)胞融合和分化的必要條件，細(xì)胞的融合、分化也將使組織器官獲取相應(yīng)的功能和合適的尺寸[23]。

表9 GO與KEGG聯(lián)合分析后的共有基因

動物機(jī)體肌肉的形成可分為多個步驟。在胚胎期，成肌細(xì)胞增殖、分化形成多核的肌管，肌管再分化形成肌纖維，最后組裝成肌肉組織[2,24-25]。肌肉發(fā)育過程需要有許多調(diào)控因子的參與，其中1、1、1、和等起著重要作用，因此本試驗(yàn)篩選這5個基因進(jìn)行詳細(xì)分析。1是由肝臟主要分泌的調(diào)控生長的一種生長因子，近來的研究表明，1也可在肌肉組織局部表達(dá)，以旁分泌/自分泌的形式，通過IGF1-Akt/PKB途徑調(diào)節(jié)骨骼肌生長[26]，還可通過IGF1-CaN-NFATc3途徑促進(jìn)成肌細(xì)胞分化和肌纖維類型的轉(zhuǎn)換[27]。MUSTN1基因在肌纖維分化、融合過程中具有關(guān)鍵作用，并調(diào)節(jié)下游靶基因1和[28]。北京鴨上的研究發(fā)現(xiàn)，1 mRNA表達(dá)水平與鴨胸腿肌的相對生長率相關(guān)[29]。1和為MRFs（生肌調(diào)控因子）的兩個重要成員。在動物胚胎期1可以誘導(dǎo)肌祖細(xì)胞向生肌細(xì)胞系轉(zhuǎn)變[30,31]，體外試驗(yàn)也表明，1可將其他類型細(xì)胞轉(zhuǎn)變?yōu)槌杉〖?xì)胞[32-34]。在成肌細(xì)胞脫離細(xì)胞周期并從增殖過程轉(zhuǎn)變到分化過程中起著重要作用，但不影響成肌細(xì)胞的增殖[35-37]。是一種重要的負(fù)調(diào)控骨骼肌生長的關(guān)鍵基因。在功能缺失的牛[38-40]、鼠[41-43]，其肌肉量顯著增加。本研究發(fā)現(xiàn)，在21胚齡到27胚齡胸肌發(fā)育阻滯的同時，1顯著上調(diào)（高郵鴨和金定鴨基因表達(dá)水平分別提高6.5倍和4.7倍），表明1在鴨胸肌發(fā)育中起著重要作用。1與1的mRNA表達(dá)模式相似，21胚齡到27胚齡也表現(xiàn)顯著上調(diào)，而1和均顯著下調(diào)，均沒有顯著變化。本試驗(yàn)前期研究中對21胚齡、27胚齡高郵鴨、金定鴨胸肌肌肉生長相關(guān)基因進(jìn)行了qRT-PCR定量分析，發(fā)現(xiàn)兩品種21胚齡到27胚齡期間1均顯著下調(diào)，均略微上調(diào)，表達(dá)在21胚齡和27胚齡兩個時間點(diǎn)均沒有品種差異[13]。在對北京鴨的研究中，Gu等發(fā)現(xiàn)21胚齡到27胚齡階段，表達(dá)水平持續(xù)下降，而表達(dá)水平持續(xù)上升[7]。除了基因，1、和1結(jié)果都與本試驗(yàn)RNA-seq結(jié)果相似。1的上調(diào)將促使更多的肌祖細(xì)胞向成肌細(xì)胞轉(zhuǎn)化，這為出雛后肌纖維的快速生長提供了更多的材料來源。細(xì)胞周期與DNA復(fù)制相關(guān)基因mRNA表達(dá)水平顯著下降，這可能導(dǎo)致肌肉衛(wèi)星細(xì)胞的有絲分裂活性降低；另外1、的下調(diào)，最終使得鴨胚在接近出雛時胸肌發(fā)育阻滯。

4 結(jié)論

本試驗(yàn)利用高通量測序技術(shù)對高郵鴨、金定鴨胚胎中后期胸肌組織的轉(zhuǎn)錄組進(jìn)行了測序分析，結(jié)果表明在此階段胸肌組織中與能量代謝相關(guān)的基因顯著上調(diào)，與DNA復(fù)制、細(xì)胞周期相關(guān)的基因顯著下調(diào)，一些肌肉發(fā)育相關(guān)的關(guān)鍵基因，如1、顯著下調(diào)，最終導(dǎo)致胸肌肌肉細(xì)胞增殖減緩，胸肌發(fā)育阻滯。該發(fā)現(xiàn)將為下一步深入探索鴨胸肌發(fā)育阻滯機(jī)制奠定基礎(chǔ)。

[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] 胡艷, 劉宏祥, 單艷菊, 姬改革, 束婧婷, 徐文娟, 朱春紅, 陶志云, 李慧芳. 鴨發(fā)育早期骨骼肌異步發(fā)育和IGF-1/MSTN-A表達(dá)的相關(guān)性. 中國農(nóng)業(yè)科學(xué), 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.

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

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

國家自然科學(xué)基金（31172194）、江蘇現(xiàn)代農(nóng)業(yè)（水禽）產(chǎn)業(yè)技術(shù)體系建設(shè)專項(xiàng)資金（JATSI2018J225）、江蘇現(xiàn)代農(nóng)業(yè)重點(diǎn)項(xiàng)目（BE2017349）

劉宏祥，E-mail：lhxatyz@foxmail.com。徐文娟，E-mail：xuwj1980@126.com。劉宏祥和徐文娟為同等貢獻(xiàn)作者。

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