犏牛體外受精胚胎的發育轉錄組分析

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

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

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

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

1 材料與方法

1.1 犏牛胚胎的生產

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

1.2 胚胎總RNA提取、cDNA文庫構建及測序

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

1.3 轉錄組數據分析

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

1.4 差異表達基因熒光定量 PCR驗證

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

表1驗證RNA-Seq的qRT-PCR基因及其引物

Table1GenesandqPCRprimersusedforRNA-Seqvalidation

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

2 結 果

2.1 體外受精效果與測序質量分析

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

表2測序和比對結果概述

Table2TheresultsofRNA-seqandmappingtothereferencegenome

比對數據Mappeddata2-細胞2-cell4-細胞4-cell8-細胞8-cell桑椹胚Morula囊胚Blastocyst干凈數據序列數CleanReads6088044663332216477918505492401849296682比對上基因組的序列數MappedReads5547920556339924430873724899907139435615比對上基因組的Reads比例/%Mappingrate91.1388.9690.1689.2180.00未比對到基因組的Reads數UnmappedReads54012416992292470447859249479861067比對到基因組多位點的Reads數Multi-mapReads25187072643767221959321191722583620比對到基因組多位點的Reads比例/%Multi-maprate4.144.174.643.865.24

2.2 犏牛胚胎的基因表達的分析

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

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

表3犏牛胚胎在不同發育階段的基因表達分析統計表

Table3Geneexpressionatthedifferentdevelopmentalstagesofcrossbredembryosoftheyak

項目Item2-細胞2-cell4-細胞4-cell8-細胞8-cell桑椹胚Morula囊胚Blastocyst表達的基因數Expressedgenes144911540013792960411860新轉錄本Noveltranscripts496564747441519首次表達的基因數1sttimeexpressedgenes-1606390148535差異表達基因數Diferetiallyexpressedgenes36906332896510298上調基因數Up-regulatedgenes2861205829085747下調基因數Down-regulatedgenes829427460574551

2.3 差異表達基因的 qRT-PCR驗證

為了進一步驗證轉錄組分析結果的可靠性，選出4個基因，采用qRT-PCR 技術進行表達驗證，結果表明，RNA-seq測序結果與qRT-PCR的結果基本一致(圖 2)。說明本研究利用RNA-seq技術分析基因的結果是可靠的。

圖2 差異表達基因的qRT-PCR與RNA-seq結果的對比Fig.2 Comparison of the differentially expressed genes between qRT-PCR and RNA-seq

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

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

2.5 差異表達基因的代謝路徑分析

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

3 討 論

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

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

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

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

Table4Top5ofenrichedKEGGpathwaysofDEGsduringdevelopmentofcrossbredembryosoftheyak

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

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

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

4 結 論

本研究首次利用RNA-Seq技術從轉錄組學揭示犏牛早期胚胎發育機制, 獲得了眾多差異基因和有關通路的富集，為完善犏牛胚胎體外生產技術提供理論基礎。

[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] 鄭 杰, 蒲思穎, 楊遠瀟,等.基于高通量測序的犏牛囊胚玻璃化冷凍損傷機制研究[J].畜牧獸醫學報，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] 石 仙, 熊顯榮, 蘭道亮, 等. 糖對牦牛卵母細胞體外成熟及其發育能力的影響[J]. 中國畜牧獸醫, 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] 羅 斌. 牦牛體外受精胚胎的發育轉錄組研究[D]. 成都: 西南民族大學, 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.