犏牛體外受精胚胎的發(fā)育轉(zhuǎn)錄組分析

字向東，劉 霜，夏 威，熊顯榮，黃 林，張正帆，李志雄，李 鍵，鐘金城，王 利，朱江江

(1.西南民族大學(xué)生命科學(xué)與技術(shù)學(xué)院，成都 610041；2.西南民族大學(xué)青藏高原研究院，成都 610041)

牦牛(Bosgrunniens)是分布在海拔2 000～5 000 m以青藏高原為中心，及其毗鄰高山、亞高山地區(qū)的特有牛種。對(duì)高寒低氧的獨(dú)特適應(yīng)能力決定了它在青藏高原畜牧業(yè)中不可替代的地位[1]。但是，相對(duì)于普通牛品種而言，牦牛生產(chǎn)生長速度慢，產(chǎn)肉性能和產(chǎn)奶性能都很低。利用普通牛(Bostaurus)優(yōu)良奶牛品種雜交改良牦牛繁殖的后代(犏牛)產(chǎn)奶性能成倍提高，深受牧民歡迎。但是，牦牛與普通牛種間雜交受胎率低[1-2]，體外受精(Invitrofertilization, IVF)的雜種胚胎發(fā)育率也不高[3-4]，其原因還不清楚。因此，研究犏牛胚胎發(fā)育調(diào)控機(jī)制，對(duì)提高牦牛種間雜交效率，完善犏牛胚胎體外生產(chǎn)技術(shù)具有十分重要的意義。

胚胎發(fā)育是眾多基因表達(dá)在時(shí)間和空間上的聯(lián)系及配合共同作用的結(jié)果[5-8]，逐一研究每個(gè)候選基因的表達(dá)模式與胚胎發(fā)育的關(guān)系從效率上和可行性上都受到限制。目前，犏牛胚胎發(fā)育調(diào)控機(jī)制方面的研究尚屬空白。隨著RNA高通量測序技術(shù)的發(fā)展[9-10]和牦?；蚪M測序的完成[11]，為從組學(xué)水平研究犏牛胚胎發(fā)育的分子機(jī)制和保存提供了快速、有效的方法。因此，本研究擬用娟珊牛精子體外受精牦牛卵母細(xì)胞，通過體外培養(yǎng)獲得雜種胚胎，然后利用微量RNA高通量測序技術(shù)分析不同發(fā)育階段的犏牛胚胎轉(zhuǎn)錄組，揭示犏牛早期胚胎發(fā)育機(jī)制, 為提高犏牛胚胎體外生產(chǎn)效率提供理論基礎(chǔ)。

1 材料與方法

1.1 犏牛胚胎的生產(chǎn)

在屠宰場采集的牦牛卵巢放入盛有DPBS液(29～33 ℃)的保溫瓶中, 2 h內(nèi)運(yùn)回到實(shí)驗(yàn)室。參照X. Xiao等[4]的方法開展牦牛卵母細(xì)胞體外成熟(Invitromaturation, IVM)：將從直徑為2～8 mm的卵泡中抽取的卵丘-卵母細(xì)胞復(fù)合體(COC)放入含10% FCS、5 μg·mL-1FSH、50 IU·mL-1LH和1 μg·mL-117β-E2的TCM199液中，在38.5 ℃、5%CO2、飽和濕度的CO2培養(yǎng)箱中成熟培養(yǎng)24 h。采用娟珊牛冷凍精液進(jìn)行體外受精(IVF)。參照石仙等[12]的方法：用Vitrolife(瑞典)體外胚胎生產(chǎn)系列試劑在38.5 ℃、5% CO2、5% O2、90% N2、飽和濕度的三氣培養(yǎng)箱中進(jìn)行精子體外獲能、IVF和胚胎的體外培養(yǎng)(Invitroculture, IVC)。在受精48 h的2-細(xì)胞、68～72 h的4-細(xì)胞、88～96 h的8-細(xì)胞，每個(gè)階段收集10枚發(fā)育良好的胚胎；受精第6 天的桑椹胚和第7天囊胚各收集1枚，提取RNA。

1.2 胚胎總RNA提取、cDNA文庫構(gòu)建及測序

采用RNeasy Micro Kit分別提取胚胎總RNA后，用Smart-Seq2方法[9]進(jìn)行擴(kuò)增：加入反應(yīng)buffer、反轉(zhuǎn)錄酶、含公共序列的Oligo-dT引物和TSO引物，反應(yīng)得到第一條cDNA鏈；不轉(zhuǎn)管，直接加入含共同序列的ISPCR引物和PCR擴(kuò)增試劑，反應(yīng)得到長度1～2 kb的第二條cDNA鏈。對(duì)富集擴(kuò)增的cDNA采用Agilent 2100 High Sensitivity DNA Assay Kit檢測擴(kuò)增產(chǎn)物cDNA樣品片段分布情況。峰圖顯示在1～2 kb片段長度存在明顯的目的產(chǎn)物主峰，表明原始樣本完整性較好，可以進(jìn)行文庫構(gòu)建。每個(gè)樣本各選取20 ng擴(kuò)增產(chǎn)物cDNA作為起始樣本構(gòu)建文庫。使用Bioruptor?Sonication System(Diagenode Inc.)進(jìn)行樣本cDNA片段化，使之?dāng)嗔褳?00 bp左右的小片段。經(jīng)超聲打斷后，進(jìn)行樣本cDNA末端修復(fù)、加 poly(A)并連接測序接頭，每步反應(yīng)后使用Beckman Ampure XP磁珠進(jìn)行純化。取接頭產(chǎn)物進(jìn)行PCR擴(kuò)增，樣本分別引入不同的Index標(biāo)簽。最后用2%(w/v)瓊脂糖凝膠電泳檢測PCR擴(kuò)增產(chǎn)物后,切取DNA片段的凝膠塊，使用CWBIO Gel Extraction Kit快速瓊脂糖凝膠DNA回收試劑盒回收DNA，再溶于EB緩沖液中，獲得測序文庫。用HiSeqTM2500進(jìn)行高通量測序。

1.3 轉(zhuǎn)錄組數(shù)據(jù)分析

llumina HisSeqTM2500測序產(chǎn)生的原始圖像數(shù)據(jù)經(jīng)Base calling轉(zhuǎn)化為序列數(shù)據(jù)(Raw Reads)，結(jié)果以FASTQ文件格式存儲(chǔ)，包含Reads的序列以及堿基的測序質(zhì)量。Raw Reads經(jīng)公共序列污染及酶切位點(diǎn)、Adapter污染、含N過多(>5%)、低質(zhì)量(Q≤19)以及比對(duì)到rRNA庫上的Reads過濾后，獲得純凈 Reads(Clean Reads)。采用TopHat v2.0.12軟件[13]將測序過濾后所得的Clean Reads與牦牛基因組(http://me.lzu.edu.cn/yak)[11]進(jìn)行比對(duì)分析，獲得唯一比對(duì)上參考基因的 Reads(Unique reads)。采用RPKM法(每百萬 Reads 中來自于某基因每千堿基長度的 Reads 數(shù)，Reads per kilobase transcriptome per million mapped Reads)計(jì)算基因表達(dá)量[14]，參照無生物學(xué)重復(fù)樣品的DEGseq法基因差異表達(dá)分析[15]，以|log2Ratio|≥1和Q<0.05為閾值, 篩選出不同發(fā)育階段間的差異表達(dá)基因(Differentially expressed genes, DEGs)。將DEGs與參考基因比較，得出DEGs顯著富集的Gene Ontology(GO)功能條目，并進(jìn)一步挖掘出與DEGs顯著相關(guān)的生物學(xué)功能。然后將DEGs向GO數(shù)據(jù)庫(http://www.geneontology.org/)各個(gè)條目進(jìn)行映射，計(jì)算其數(shù)目，用Y. Benjamini和Y. Hochberg方法[16]對(duì)P-value進(jìn)行校正后，以Q-value<0.05為閾值定義DEGs中顯著富集的GO條目。通過與京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes，KEGG)數(shù)據(jù)庫(http://wego.genomics.org.cn)進(jìn)行比對(duì)，對(duì)基因涉及的信號(hào)通路或代謝途徑(Pathway)進(jìn)行分析。

1.4 差異表達(dá)基因熒光定量 PCR驗(yàn)證

隨機(jī)選擇與發(fā)育相關(guān)的SKP1、CD63、ZAR1和H3等4個(gè)DEGs，采用Primer5.0軟件進(jìn)行基因的定量引物設(shè)計(jì)(表1)，以1.2中反轉(zhuǎn)錄獲得的2-、4-、8-細(xì)胞、桑椹胚和囊胚的cDNA等體積混合，然后倍比稀釋，將稀釋的樣本作為模板，構(gòu)建其熒光定量標(biāo)準(zhǔn)曲線。qRT-PCR反應(yīng)體系為10 μL：上下游引物各0.8 μL，Sso AdvancedTMSYBR?Green Super mix 5 μL，ddH2O 2.9 μL，cDNA模板0.5 μL。擴(kuò)增條件：95 ℃預(yù)變性3 min，95 ℃ 10 s、59.5 ℃ 20 s，35個(gè)循環(huán)。以H2A為內(nèi)參，采用2-ΔΔCt計(jì)算基因的差異倍數(shù)。每個(gè)樣品設(shè)3次重復(fù)。

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

Table1GenesandqPCRprimersusedforRNA-Seqvalidation

基因Gene正向引物(5'-3')Forwardprimer反向引物(5'-3')Reverseprimer產(chǎn)物長度/bpProductssizeH2AGCGTATTACCCCTCGTCACTCTTCTGTTGTCCTTTCTTTCC138SKP1ATCCAGTCCCCTTGCCAAATAGAGTGTCCCTTGGTCAACT173CD63AAGATTTTGAGTGCTGCGGGCCGCAATCTTCTCCACACAG164ZAR1GCTGGGAAAGTGCCTATGTGGAAGTTTCACTGGGCAGGAG158H3GCTCTGGAATGGGAGGTTCTTAGTACATCCCACGCCATCC219

2 結(jié) 果

2.1 體外受精效果與測序質(zhì)量分析

通過3個(gè)批次的重復(fù)試驗(yàn)，IVF/IVC后雜種胚胎的平均卵裂率為78.4%，囊胚率為36.3%，說明采用的IVF/IVC體系效果好，測序結(jié)果能代表正常胚胎轉(zhuǎn)錄組水平。通過分析堿基的組成和質(zhì)量值分布可以控制原始數(shù)據(jù)的質(zhì)量。本研究中2-、4-、8-細(xì)胞、桑椹胚和囊胚 5 個(gè)發(fā)育時(shí)期的犏牛胚胎 RNA 樣本測序所得原始數(shù)據(jù)的堿基質(zhì)量分布和堿基含量分布良好，且Q30>85%，說明測序質(zhì)量較好。原始數(shù)據(jù)經(jīng)過濾后獲得47 791 850(8-細(xì)胞)～63 332 216(4-細(xì)胞)條 Clean Reads。將 Clean Reads與牦牛基因組比對(duì)，有80.00%(囊胚)～91.13%(2-細(xì)胞)的序列可被牦?；蚪M注釋，比對(duì)到基因組多位的Reads比例為3.86%～5.24%(表2)。

表2測序和比對(duì)結(jié)果概述

Table2TheresultsofRNA-seqandmappingtothereferencegenome

比對(duì)數(shù)據(jù)Mappeddata2-細(xì)胞2-cell4-細(xì)胞4-cell8-細(xì)胞8-cell桑椹胚M(jìn)orula囊胚Blastocyst干凈數(shù)據(jù)序列數(shù)CleanReads6088044663332216477918505492401849296682比對(duì)上基因組的序列數(shù)MappedReads5547920556339924430873724899907139435615比對(duì)上基因組的Reads比例/%Mappingrate91.1388.9690.1689.2180.00未比對(duì)到基因組的Reads數(shù)UnmappedReads54012416992292470447859249479861067比對(duì)到基因組多位點(diǎn)的Reads數(shù)Multi-mapReads25187072643767221959321191722583620比對(duì)到基因組多位點(diǎn)的Reads比例/%Multi-maprate4.144.174.643.865.24

2.2 犏牛胚胎的基因表達(dá)的分析

采用RPKM法[14]計(jì)算基因表達(dá)量的結(jié)果表明，犏牛胚胎從2-細(xì)胞期發(fā)育到囊胚期每個(gè)階段有9 604(桑椹胚)～15 400(4-細(xì)胞)個(gè)表達(dá)基因；共有19 072個(gè)基因表達(dá)，其中7 785個(gè)為共表達(dá)基因，11 287個(gè)為階段特異表達(dá)基因。4-細(xì)胞開始表達(dá)的基因數(shù)目最多(1 606個(gè))，而桑椹胚開始表達(dá)的基因數(shù)目最少(148個(gè))。隨著胚胎發(fā)育的進(jìn)行，BMP15、GDF9和ZP4等母源基因的表達(dá)量降低，而ATP5B、UBEA3和SNURF等胚胎基因的開始表達(dá)且表達(dá)量不斷增加(圖1)。以|log2Ratio|≥1和Q<0.05為篩選條件篩選不同發(fā)育階段間的DEGs，結(jié)果表明,桑椹胚～囊胚的DEGs最多(10 298個(gè))，而2-～4-細(xì)胞胚的DEGs最少(3 690個(gè))(表3)。

圖1 母源基因(A)和胚胎基因(B)的表達(dá)變化趨勢Fig.1 The expression tendency of maternal genes (A) and embryonic genes (B)

表3犏牛胚胎在不同發(fā)育階段的基因表達(dá)分析統(tǒng)計(jì)表

Table3Geneexpressionatthedifferentdevelopmentalstagesofcrossbredembryosoftheyak

項(xiàng)目Item2-細(xì)胞2-cell4-細(xì)胞4-cell8-細(xì)胞8-cell桑椹胚M(jìn)orula囊胚Blastocyst表達(dá)的基因數(shù)Expressedgenes144911540013792960411860新轉(zhuǎn)錄本Noveltranscripts496564747441519首次表達(dá)的基因數(shù)1sttimeexpressedgenes-1606390148535差異表達(dá)基因數(shù)Diferetiallyexpressedgenes36906332896510298上調(diào)基因數(shù)Up-regulatedgenes2861205829085747下調(diào)基因數(shù)Down-regulatedgenes829427460574551

2.3 差異表達(dá)基因的 qRT-PCR驗(yàn)證

為了進(jìn)一步驗(yàn)證轉(zhuǎn)錄組分析結(jié)果的可靠性，選出4個(gè)基因，采用qRT-PCR 技術(shù)進(jìn)行表達(dá)驗(yàn)證，結(jié)果表明，RNA-seq測序結(jié)果與qRT-PCR的結(jié)果基本一致(圖 2)。說明本研究利用RNA-seq技術(shù)分析基因的結(jié)果是可靠的。

圖2 差異表達(dá)基因的qRT-PCR與RNA-seq結(jié)果的對(duì)比Fig.2 Comparison of the differentially expressed genes between qRT-PCR and RNA-seq

2.4 差異表達(dá)基因的 GO 分類和顯著性富集分析

GO分析顯示(圖3)，從2-～4-細(xì)胞、4-～8-細(xì)胞、8-細(xì)胞～桑椹胚及桑椹胚～囊胚4個(gè)發(fā)育階段分別有3 461(占93.8%)、5 972(占94.3%)、8 438(占94.1%)和9 749(占94.7%)個(gè)DEGs得到歸類注釋，都涉及生物過程(Biological process，BP)、細(xì)胞組分(Cellular component，CC)和分子功能(Molecular function，MF)3大類67個(gè)二級(jí)條目。在BP分類中有23個(gè)二級(jí)條目，其中占比例最大的4個(gè)二級(jí)條目在4個(gè)發(fā)育階段相同，從高到低依次為細(xì)胞過程(Cellular process)、單有機(jī)體過程(Single-organism process)、生物調(diào)節(jié)(Biological regulation)和代謝過程(Metabolic process)。但是，各發(fā)育階段第5個(gè)二級(jí)條目開始出現(xiàn)差異。在CC分類中有22個(gè)二級(jí)條目，二級(jí)條目中細(xì)胞部分(Cell part)占比例最多，其次為細(xì)胞器 (Organelle)及細(xì)胞器部件(Organelle part)，前10個(gè)二級(jí)條目的排序在各發(fā)育階段相同，但從第11個(gè)二級(jí)條目開始排序在不同發(fā)育階段出現(xiàn)差異。在MF分類中有22個(gè)二級(jí)條目，占比例最大的2個(gè)二級(jí)條目在4個(gè)發(fā)育階段完全相同，依次是綁定分子(Binding)所占比例最多和催化活性(Catalytic activity)，但從第3個(gè)二級(jí)條目開始排序在不同發(fā)育階段出現(xiàn)差異。

2.5 差異表達(dá)基因的代謝路徑分析

犏牛胚胎發(fā)育過程中差異表達(dá)基因KEGG分析表明，2-～4-細(xì)胞涉及到300條通路(Pathways)，4-～8-細(xì)胞有302條通路，8-細(xì)胞～桑椹胚有314條通路，桑椹胚～囊胚有314條通路，各發(fā)育階段的富集前5條通路(Top 5 pathways)如表4所示，每個(gè)階段的通路種類及富集性都有差異。2-～4-細(xì)胞、4-～8-細(xì)胞、8-細(xì)胞～桑椹胚和桑椹胚～囊胚4個(gè)階段分別有15、0、18和1條通路顯著富集(Q<0.05)。

3 討 論

本研究利用RNA-seq技術(shù)從轉(zhuǎn)錄組學(xué)的角度揭示了犏牛早期胚胎發(fā)育機(jī)制, 為完善犏牛胚胎體外生產(chǎn)，提高犏牛胚胎體外生產(chǎn)效率具有重要價(jià)值。普通牛從2-細(xì)胞～囊胚，每個(gè)胚胎的總RNA只有200～2 000 pg[17]。犏牛胚胎應(yīng)該與此相近，這樣微量RNA不能滿足轉(zhuǎn)錄組測序文庫的構(gòu)建及高通量測序的基本要求，因此，分別提取每個(gè)發(fā)育階段的犏牛早期胚胎總RNA后，使用Smart-Seq2擴(kuò)增技術(shù)對(duì)樣本進(jìn)行富集并構(gòu)建測序文庫[9]，再應(yīng)用RNA高通量測序技術(shù)對(duì)其進(jìn)行高通量測序分析。通過從擴(kuò)增產(chǎn)物cDNA樣品片段分布、堿基的組成和質(zhì)量值分布、Q30、飽和量、qRT-PCR驗(yàn)證等全面分析的結(jié)果表明，以犏牛胚胎微量RNA作為模板，采用Smart-Seq2擴(kuò)增技術(shù)與RNA-seq技術(shù)相結(jié)合開展犏牛胚胎轉(zhuǎn)錄組高通量測序分析，結(jié)果可靠。

犏牛胚胎從2-細(xì)胞發(fā)育到囊胚每個(gè)階段有9 604(桑椹胚)～15 400(4-細(xì)胞胚)個(gè)表達(dá)基因。與普通牛[18-19]和牦牛的胚胎[20]比較，犏牛桑椹胚的轉(zhuǎn)錄本少。這可能是由于犏牛胚胎為種間雜種胚胎，胚胎的遺傳信息決定了在這一階段不可避免地出現(xiàn)“發(fā)育阻滯”現(xiàn)象，也有可能是胚胎體外培養(yǎng)條件的不完善引起“體外發(fā)育阻滯”[21]。可以通過對(duì)體內(nèi)發(fā)育的犏牛胚胎的轉(zhuǎn)錄組進(jìn)一步深入研究，探明原因。如果確認(rèn)為“體外發(fā)育阻滯”，則可以通過優(yōu)化胚胎培養(yǎng)條件得到改善。犏牛胚胎從2-細(xì)胞發(fā)育到囊胚共有19 072個(gè)基因表達(dá)，其中7 785個(gè)為共表達(dá)基因，11 287個(gè)為階段特異表達(dá)基因。4-細(xì)胞開始表達(dá)的基因數(shù)目最多(1 606個(gè))，而桑椹胚開始表達(dá)的基因數(shù)目最少(148個(gè))，顯示胚胎發(fā)育是眾多基因表達(dá)在時(shí)間和空間上的聯(lián)系和配合共同作用的結(jié)果，不是單個(gè)基因控制的。但是，有些基因在胚胎發(fā)育的特異階段發(fā)揮關(guān)鍵的調(diào)控作用，例如，NANOG能調(diào)控合子基因組的激活, 缺乏NANOG的胚胎合子基因組激活率低，發(fā)育受阻[22]。本研究發(fā)現(xiàn)，與普通牛[7,18-19]和牦牛胚胎[20]一樣，犏牛胚胎NANOG基因在8-細(xì)胞開始表達(dá)。CLDN4是Claudin蛋白家族的成員之一，婦女黃體期子宮內(nèi)膜CLDN4 mRNA的表達(dá)量與妊娠相關(guān)[23]，囊胚中基因的表達(dá)與囊胚的發(fā)育與附植相關(guān)[24]。推測CLDN4基因在犏牛囊胚中開始表達(dá)，可能在囊胚發(fā)育和妊娠的建立中發(fā)揮重要作用。早期胚胎的發(fā)育受到由來自卵母細(xì)胞發(fā)生及成熟期間合成的大量轉(zhuǎn)錄本和蛋白質(zhì)的調(diào)控[7,18-19,25]，但隨著發(fā)育的進(jìn)行，母源轉(zhuǎn)錄本和蛋白質(zhì)降解，而胚胎基因組激活(Embryonic genome activation, EGA)，發(fā)育從母型調(diào)控向胚胎型調(diào)控的過渡(Maternal-to-embryonic transition, MET)。本研究發(fā)現(xiàn)，隨著胚胎發(fā)育的進(jìn)行，犏牛胚胎的BMP15、GDF9、ZP4和ZP3等母源基因表達(dá)量逐漸減少，而ATP5B、UBEA3和SNURF等胚胎基因的開始表達(dá)且表達(dá)量不斷增加(圖2)，這與在其它物種的研究結(jié)果一致[7,18-19,26]。

A.2-～4-細(xì)胞；B.4-～8-細(xì)胞；C.8-細(xì)胞～桑椹胚；D.桑椹胚～囊胚A.2-cell-4-cell; B.4-cell-8-cell; C.8-cell-morula; D.Morula-blastocyst圖3 犏牛早期胚胎發(fā)育過程中的差異表達(dá)基因GO分類注釋圖Fig.3 Gene Ontology classification of the DEGs of 4 developmental stages of crossbred embryos of the yak

表4牦牛胚胎發(fā)育過程中差異表達(dá)基因富集前5(Top5)KEGG通路

Table4Top5ofenrichedKEGGpathwaysofDEGsduringdevelopmentofcrossbredembryosoftheyak

發(fā)育階段DevelopmentstageKEGG通路KEGGpathway上調(diào)基因Up-regulatedgenes下調(diào)基因Down-regulatedgenesQ-值Q-value2-細(xì)胞～4-細(xì)胞2-cell-4-cell泛素介導(dǎo)調(diào)節(jié)Ubiquitinmediatedproteolysis1020.0005內(nèi)質(zhì)網(wǎng)加工蛋白Proteinprocessinginendoplasmicreticulum1920.0081細(xì)胞周期Cellcycle1410.0152剪切體Spliceosome1410.0176RNA轉(zhuǎn)運(yùn)RNAtransport1630.02364-細(xì)胞～8-細(xì)胞4-cell-8-cell基底細(xì)胞癌Basalcellcarcinoma7190.0572RNA轉(zhuǎn)運(yùn)RNAtransport10300.0572系統(tǒng)性紅斑狼瘡Systemiclupuserythematosus33220.0798神經(jīng)活性配體-受體互作Neuroactiveligand-receptorinteraction20480.0798剪切體Spliceosome15190.07988-細(xì)胞～桑椹胚8-cell-morula核糖體Ribosome13561.63E-13氧化磷酸化機(jī)制Oxidativephosphorylation21311.82E-06亨廷頓氏病Huntington'sdisease31531.88E-06RNA轉(zhuǎn)運(yùn)RNAtransport18461.52E-05神經(jīng)活性配體-受體互作Neuroactiveligand-receptorinteraction18538.03E-05桑椹胚～囊胚M(jìn)orula-blastocystRNA轉(zhuǎn)運(yùn)RNAtransport61340.0387核糖體Ribosome17340.1013內(nèi)質(zhì)網(wǎng)加工蛋白Proteinprocessinginendoplasmicreticulum72390.1799核苷酸切除修復(fù)Nucleotideexcisionrepair17100.1799氨基酸合成Biosynthesisofaminoacids29100.1799

以|log2Ratio|≥1和Q<0.05為篩選條件篩選出桑椹胚～囊胚的DEGs最多(10 298個(gè))，而2-～4-細(xì)胞胚的DEGs最少(3 690個(gè))(表3)。GO分析顯示，從2-～4-細(xì)胞、4-～8-細(xì)胞、8-細(xì)胞～桑椹胚及桑椹胚～囊胚4個(gè)發(fā)育階段分別有3 461、5 972、8 438和9 749個(gè)DEGs(占DEGs數(shù)目的94%左右)得到歸類注釋，都涉及生物過程(BP)、細(xì)胞組分(CC)和分子功能(MF)3大類67個(gè)二級(jí)條目(圖3)。在BP分類中，比例最大的4個(gè)二級(jí)條目在4個(gè)發(fā)育階段都相同，但各發(fā)育階段第5個(gè)二級(jí)條目開始出現(xiàn)差異；在CC分類中前10個(gè)二級(jí)條目的排序在4個(gè)發(fā)育階段都相同，但從第11個(gè)二級(jí)條目開始在不同發(fā)育階段出現(xiàn)差異；在MF分類中，只有占比例最大的2個(gè)二級(jí)條目在4個(gè)發(fā)育階段都相同，但從第3個(gè)二級(jí)條目開始排序在不同發(fā)育階段出現(xiàn)差異。說明在犏牛胚胎發(fā)育過程中，胚胎細(xì)胞組分相對(duì)穩(wěn)定，但細(xì)胞內(nèi)分子過程和分子功能在不斷地發(fā)生變化。例如，分子功能調(diào)節(jié)器(Molecular function regulator)具有調(diào)節(jié)基因產(chǎn)物活性的功能，它從2-～4-細(xì)胞和4-～8-細(xì)胞位列BF分類中的第5位上調(diào)到8-細(xì)胞～桑椹胚中的第4位和桑椹胚～囊胚中的第3位。

本研究發(fā)現(xiàn), 犏牛早期胚胎發(fā)育過程中DEGs富集的主要通路有剪接體(Spliceosome)、神經(jīng)活性配體-受體互作(Neuroactive ligand-receptor interaction)、細(xì)胞因子及其受體相互作用(Cytokine-cytokine receptor interactions)、泛素介導(dǎo)的蛋白水解(Ubiquitin mediated proteolysis)、RNA轉(zhuǎn)運(yùn)(RNA transport)、核糖體(Ribosome)等(表4)，這與在其它動(dòng)物的發(fā)現(xiàn)基本一致[22,27-28]，說明這些通路在犏牛胚胎發(fā)育過程中發(fā)揮重要調(diào)控作用。

4 結(jié) 論

本研究首次利用RNA-Seq技術(shù)從轉(zhuǎn)錄組學(xué)揭示犏牛早期胚胎發(fā)育機(jī)制, 獲得了眾多差異基因和有關(guān)通路的富集，為完善犏牛胚胎體外生產(chǎn)技術(shù)提供理論基礎(chǔ)。

[1] WIENER G, HAN J L, LONG R J. The yak[M]. 2nd ed. Bangkok: Regional Office for Asia and the Pacific of the Food and Agriculture Organization of the United Nations, 2003: 1-13.

[2] ZI X D, HE S M, LU H, et al. Induction of estrus in suckled female yaks (Bosgrunniens) and synchronization of ovulation in the non-sucklers for timed artificial insemination using progesterone treatments and Co-Synch regimens[J].AnimReprodSci, 2006, 92(1-2): 183-192.

[3] ZI X D, YIN R H, CHEN S W, et al. Developmental competence of embryos derived from reciprocalinvitrofertilization between yak (Bosgrunniens) and cattle (Bostaurus)[J].JReprodDev, 2009, 55(5): 480-483.

[4] XIAO X, ZI X D, NIU H R, et al. Effect of addition of FSH, LH and proteasome inhibitor MG132 toinvitromaturation medium on the developmental competence of yak (Bosgrunniens) oocytes[J].ReprodBiolEndocrinol, 2014, 12: 30.

[5] YAN L Y, YANG M Y, GUO H S, et al. Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells[J].NatStructMolBiol, 2013, 20(9): 1131-1139.

[6] CAO S Y, HAN J Y, WU J, et al. Specific gene-regulation networks during the pre-implantation development of the pig embryo as revealed by deep sequencing[J].BMCGenomics, 2014, 15: 4.

[7] GRAF A, KREBS S, HEININEN-BROWN M, et al. Genome activation in bovine embryos: review of the literature and new insights from RNA sequencing experiments[J].AnimReprodSci, 2014, 149(1-2): 46-58.

[8] JIANG Z L, SUN J W, DONG H, et al. Transcriptional profiles of bovineinvivopre-implantation development[J].BMCGenomics, 2014, 15: 756.

[9] PICELLI S, BJ?RKLUND A K, FARIDANI O R, et al. Smart-Seq2 for sensitive full-length transcriptome profiling in single cells[J].NatMethods, 2013, 10(11): 1096-1098.

[10] 鄭 杰, 蒲思穎, 楊遠(yuǎn)瀟,等.基于高通量測序的犏牛囊胚玻璃化冷凍損傷機(jī)制研究[J].畜牧獸醫(yī)學(xué)報(bào)，2017，48(10)：1871-1881.

ZHENG Jie, PU Si-ying, YANG Yuan-xiao, et al.Exploring Mechanism for Vitrification Damage of the Cross-bred Blastocysts of the Yak via High-throughput Sequencing[J].ActaVeterinariaetZootechnicaSinica, 2017, 48(10): 1871-1881. (in Chinese)

[11] QIU Q, ZHANG G J, MA T, et al. The yak genome and adaptation to life at high altitude[J].NatGenet, 2012, 44(8): 946-949.

[12] 石 仙, 熊顯榮, 蘭道亮, 等. 糖對(duì)牦牛卵母細(xì)胞體外成熟及其發(fā)育能力的影響[J]. 中國畜牧獸醫(yī), 2017, 44(1): 155-160.

SHI X, XIONG X R, LAN D L, et al. Effect of sugar oninvitromaturation and developmental competence of yak oocytes[J].ChinaAnimalHusbandry&VeterinaryMedicine, 2017, 44(1): 155-160. (in Chinese)

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

[14] WANGER G P, KIN K, LYNCH V J. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples[J].TheoryBiosci, 2012, 131(4): 281-285.

[15] WANG L K, FENG Z X, WANG X, et al. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data[J].Bioinformatics, 2010, 26(1): 136-138.

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

[17] GILBERT I, SCANTLAND S, SYLVESTRE E L, et al. The dynamics of gene products fluctuation during bovine pre-hatching development[J].MolReprodDev, 2009, 76(8): 762-772.

[18] KUES W A, SUDHEER S, HERRMANN D, et al. Genome-wide expression profiling reveals distinct clusters of transcriptional regulation during bovine preimplantation developmentinvivo[J].ProcNatlAcadSciUSA, 2008, 105(50): 19768-19773.

[19] GRAF A, KREBS S, ZAKHARTCHENKO V, et al. Fine mapping of genome activation in bovine embryos by RNA sequencing[J].ProcNatlAcadSciUSA, 2014, 111(11): 4139-4144.

[20] 羅 斌. 牦牛體外受精胚胎的發(fā)育轉(zhuǎn)錄組研究[D]. 成都: 西南民族大學(xué), 2016.

LUO B. Transcriptomic analysis of the regulation of yak embryonic developmentinvitro[D]. Chengdu: Southwest University for Nationalities, 2016. (in Chinese)

[21] GARCIA S M, MARINHO L S R, LUNARDELLI P A, et al. Developmental block and programmed cell death inBosindicusembryos: effects of protein supplementation source and developmental kinetics[J].PLoSOne, 2015, 10(3): e0119463.

[22] LEE M T, BONNEAU A R, TAKACS C M, et al. Nanog, Pou5f1 and SoxB1 activate zygotic gene expression during the maternal-to-zygotic transition[J].Nature, 2013, 503(7476): 360-364.

[23] SERAFINI P C, SILVA I D, SMITH G D, et al. Endometrial claudin-4 and leukemia inhibitory factor are associated with assisted reproduction outcome[J].ReprodBiolEndocrinol, 2009, 7: 30.

[24] MUNCH E M, SPARKS A E, GONZALEZ BOSQUET J, et al. Differentially expressed genes in preimplantation human embryos: potential candidate genes for blastocyst formation and implantation[J].JAssistReprodGenet, 2016, 33(8): 1017-1025.

[25] TADROS W, LIPSHITZ H D. The maternal-to-zygotic transition: a play in two acts[J].Development, 2009, 136(18): 3033-3042.

[26] BRAUDE P, BOLTON V, MOORE S. Human gene expression first occurs between the four- and eight-cell stages of preimplantation development[J].Nature, 1988, 332(6163): 459-461.

[27] ZUO Y C, SU G H, WANG S S, et al. Exploring timing activation of functional pathway based on differential co-expression analysis in preimplantation embryogenesis[J].Oncotarget, 2016, 7(45): 74120-74131.

[28] RAUWERDA H, PAGANO J F B, DE LEEUW W C, et al. Transcriptome dynamics in early zebrafish embryogenesis determined by high-resolution time course analysis of 180 successive, individual zebrafish embryos[J].BMCGenomics, 2017, 18: 287.