基于RNA-Seq技術的牦牛體外受精胚胎發育轉錄組分析

字向東，羅斌，夏威，鄭玉才，熊顯榮，李鍵，鐘金城，朱江江，張正帆

?

基于RNA-Seq技術的牦牛體外受精胚胎發育轉錄組分析

字向東1，羅斌1，夏威1，鄭玉才1，熊顯榮1，李鍵1，鐘金城2，朱江江2，張正帆1

（1西南民族大學生命科學與技術學院, 成都 610041；2西南民族大學青藏高原研究院, 成都 610041）

【目的】探明不同發育階段牦牛體外受精（IVF）胚胎的轉錄組差異，了解差異表達基因（DEG）在功能、分類和代謝通路的差異，為揭示牦牛早期胚胎發育調控機制，促進牦牛胚胎體外生產技術的發展提供理論基礎。【方法】以IVF技術生產的牦牛2-細胞、4-細胞、8-細胞、桑椹胚和囊胚5個發育階段的胚胎為樣本分別提取總RNA，采用Smart-Seq2擴增技術構建5個測序文庫，應用HiSeqTM2500高通量測序技術進行轉錄組測序，對獲得的有效序列進行功能注釋及相關生物信息學分析。【結果】牦牛IVF后的卵裂率和囊胚率分別為69.3%和26.2%。2-細胞、4-細胞、8-細胞、桑椹胚和囊胚5個發育階段的牦牛胚胎Clean reads為47 355 570—50 855 888條，其中有85.65%—90.02%的reads能比對到牦牛參考基因組序列上；8-細胞的胚胎比對上的轉錄本最多（14 893），而囊胚比對上的轉錄本最少（9 827）。牦牛胚胎轉錄本主要有5種可變剪接類型，轉錄起始區域可變剪接（TSS）和轉錄結束區域可變剪接（TTS）所占比例最大；牦牛2-細胞、4-細胞、8-細胞、桑椹胚和囊胚基因組分別有116 601、234 131、196 420、70 841和94 840個位點存在單核苷酸多態性（SNP）。在4-細胞、8-細胞、桑椹胚、囊胚期開始表達的基因分別有1 221、1 116、142和564個；隨著胚胎發育的進行，、、、、和等母源基因的表達量逐漸減少，而、、、、、和等胚胎基因的表達量則逐漸增加。以|log2ratio| ≥1且Q-value ＜ 0.05為篩選標準，在2-細胞和4-細胞胚胎、4-細胞和8-細胞胚胎、8-細胞胚胎和桑椹胚以及桑椹胚和囊胚比對中分別篩選到6 922、7 601、8 071和10 555個DEGs。GO分析表明，4個發育階段的DEGs歸類注釋都涉及生物過程（BP）、細胞組分（CC）和分子功能（MF）3大類62個二級條目。通過KEGG pathway數據庫分析，2-細胞和4-細胞期胚胎的DEGs參與到308條通路中，顯著富集剪接體、RNA轉運和泛素介導的蛋白水解等11條通路；4-細胞和8-細胞胚胎的DEGs參與到310條通路，顯著富集嗅覺轉導、神經活性配體-受體互作和核苷酸切除修復等9條通路；8-細胞期胚胎和桑椹胚的DEGs參與到316條通路，顯著富集嗅覺轉導、泛素介導的蛋白水解和神經活性配體-受體互作等10條通路；桑椹胚和囊胚的DEGs參與到315條通路，顯著富集剪接體和RNA轉運2條通路。【結論】利用高通量測序技術對牦牛IVF胚胎不同發育階段的轉錄組進行測序和分析，揭示了牦牛胚胎發育不同階段差異表達基因的數量，獲得了差異表達基因的功能、分類和代謝通路。為豐富牦牛胚胎轉錄組信息，揭示牦牛胚胎發育分子調控機制及完善其胚胎體外培養技術研究奠定基礎。

牦牛；體外受精；胚胎；轉錄組；RNA-Seq

0 引言

【研究意義】牦牛（）是分布在海拔2 000―5 000 m以青藏高原及其毗鄰高山、亞高山地區的珍奇稀有牛種。目前，全世界有約1 400余萬頭牦牛，其中，中國擁有的牦牛數量占世界95%以上[1]。牦牛是我國青藏高原人民賴以生存的生產和生活資料，對高寒低氧的獨特適應能力決定了該畜種在青藏高原畜牧業中占據不可替代的特殊地位。牦牛的繁殖率低（40%―60%）[1-2]，體外受精（fertilization, IVF）的囊胚發育率也不高（10%―30%）[3-4]。因此，研究牦牛早期胚胎發育調控機理，對提高牦牛繁殖效率及完善牦牛胚胎體外生產（production, IVP）技術具有重要意義。【前人研究進展】轉錄組（transcriptome）是特定組織或細胞在某一環境或生理條件下所轉錄表達的所有 RNA 的總和，既包括編碼蛋白的 mRNA，也包括非編碼RNA，是連接基因組遺傳信息與蛋白質組生物功能的紐帶。因此，轉錄組測序技術（RNA sequencing，RNA-Seq）為研究基因表達及調控提供了重要的手段和方法，得到廣泛使用[5-7]。雖然卵母細胞、早期胚胎的 RNA含量低[8]，難以達到轉錄組測序建庫的RNA最小起始量要求，但是，隨著Smart-seq2等微量轉錄組學測序技術的建立與發展[9]，研究者已成功利用RNA-Seq技術解析了人[10]、豬[11]和牛[12-14]等物種的早期胚胎發育調控機制。雖然目前已利用RNA-Seq技術在牦牛卵巢[15]、卵母細胞成熟[16]、胚胎冷凍損傷[17]及犏牛雄性不育[18]的分子調控機制方面開展過一些初步研究，但對牦牛胚胎發育調控機制方面的研究尚未見報道。【本研究切入點】胚胎發育是眾多基因表達在時間和空間上的聯系和配合共同作用的結果[10-13,19]。然而，迄今為止，與牦牛胚胎發育調控相關的功能基因的研究非常有限，而逐一研究每個基因的表達模式與胚胎發育的關系從效率上和可行性上都不太理想。牦牛基因組測序的成果[20]及高通量測序技術的發展，為從組學水平研究牦牛胚胎發育的分子機理提供了快速、有效的方法。【擬解決的關鍵問題】本研究采用Smart-Seq2方法對牦牛IVF 2-、4-、8-細胞、桑椹胚和囊胚總RNA擴增建庫，應用RNA-Seq技術進行高通量測序分析，以期為揭示牦牛早期胚胎發育調控機制，完善牦牛IVP奠定基礎。

1 材料與方法

1.1 試劑與耗材

M199（10×）、17β-雌二醇（17β-E2）和丙酮酸鈉（Sigma，美國）；FSH（Folltropin?-V）和LH（Lutropin?-V）（Bioniche，加拿大）；胎牛血清（FCS）（Gibco）；SpermRinseTM、G-IVFTM、G-1TM、G-2TM和透明質酸酶（HYASE-10×）（Vitrolife，瑞典）；細胞裂解液和RNase Inhibitor（安諾優達基因科技有限公司，北京）；細胞裂解液（Sigma-Aldrich，美國）；SMARTer PCR cDNA synthesis Kit（Takara）；Agilent 2100 High Sensitivity DNA Assay Kit（Agilent Technologies，美國）；Gel Extraction Kit（CWBIO）。

1.2 牦牛胚胎的體外生產（IVP）

參照Xiao等[4]建立的牦牛IVP技術生產牦牛胚胎。①牦牛卵母細胞體外成熟（IVM）：在10―12月份從成都郊區屠宰場采集牦牛卵巢，2 h內帶回到西南民族大學動物科學國家民委重點實驗室，用無菌生理鹽水清洗3次。從卵巢表面直徑為2―8 mm的卵泡中抽取卵丘-卵母細胞復合體（COC），放入含10% （v/v）FCS、5 μg·mL-1FSH、50 IU·mL-1LH和1 μg·mL-117β-E2的TCM199液中，在38.6℃、5%CO2、飽和濕度的CO2培養箱中成熟培養24 h。②精子獲能：牦牛細管凍精在35℃水浴解凍后，200 μL精液加入到含1 mL精子SpermRinseTM獲能液的1.5 mL無菌EP管底部，獲能培養50 min。取500 μL上清液，500 g離心5 min，棄上清液，留底部50 μL用于IVF。③IVF與胚胎體外培養（culture, IVC）：將7―10 μL獲能精液加入含50―80個成熟COCs的100 μL G-IVFTM受精液微滴中，精卵共培養（受精）22 h即完成IVF。將受精卵用0.2%透明質酸酶去除卵丘細胞后，先后移入G-1TM和G-2TM胚胎培養液微滴IVC。精子獲能、IVF和IVC均在5% CO2、5% O2、90% N2、飽和濕度的三氣培養箱中完成。收集5個階段的胚胎開展牦牛早期胚胎發育轉錄組分析，樣本包括: 10個受精42―46 h的2-細胞胚胎、10個68―72 h的4-細胞胚胎、8個88―96 h的8-細胞胚胎、5個5―6d的桑椹胚、3個7d的囊胚。

1.3 RNA提取、測序文庫構建及測序

測序文庫構建及測序由安諾優達基因科技（北京）有限公司完成。首先采用細胞裂解液分別裂解2-、4-、8-細胞、桑椹胚和囊胚等5個發育階段的牦牛胚胎，釋放總RNA，然后，使用Smart-Seq2方法[7]進行擴增: 加入反應buffer、反轉錄酶、含公共序列的oligo-dT引物和TSO引物，反應得到一鏈擴增產物；不轉管，直接加入含共同序列的ISPCR引物和PCR擴增試劑，反應得到長度約1―2 kb的二鏈擴增產物cDNA。采用Agilent 2100 High Sensitivity DNA Assay Kit檢測富集擴增獲得的cDNA產物片段分布情況，根據檢測結果判定擴增產物cDNA質量，確定后續文庫構建。

每個樣本各選取20 ng擴增產物cDNA作為起始原料進行文庫構建。使用Bioruptor?Sonication System（Diagenode Inc.）進行樣本cDNA片段化處理，使之斷裂為200 bp左右的小片段，然后用2%（w/v）瓊脂糖凝膠電泳檢測片段化處理效果。經超聲打斷后，進行樣本cDNA末端修復、加堿基A、測序接頭拼接等各步反應，每步反應后使用Beckman Ampure XP磁珠進行純化。取接頭產物進行PCR擴增，樣本分別引入不同的Index標簽，便于上機測序時區分。最后，用2%（w/v）瓊脂糖凝膠電泳檢測PCR擴增產物，切取DNA片段的凝膠塊，使用CWBIO Gel Extraction Kit快速瓊脂糖凝膠DNA回收試劑盒回收DNA，再溶于EB緩沖液中，即為最終的文庫。最后采用測序策略為PE125模式，使用HiSeqTM2500測序平臺對構建好的文庫進行高通量測序。

1.4 轉錄組數據分析

為保證測序數據質量，對HiSeqTM2500測序所得的Raw reads進行數據過濾，去除接頭序列、空讀序列以及低質量序列（Phred quality ＜5）后，得到Clean reads。采用TopHat v2.0.12軟件[21]將Clean reads與牦牛參考基因組（https://www.ncbi.nlm.nih.gov /genome/?term=yak）[20]進行比對分析。利用RPKM法（Reads per kilobase transcriptome per million mapped reads）計算基因表達量[22]，無生物學重復樣品的DEGseq法進行兩個連續發育階段之間基因差異表達分析[23]，以|log2ratio|≥1和Q-value＜0.05作為基因差異表達的閾值篩選差異表達基因（differentially expressed genes，DEGs）。將獲得的DEGs向GO數據庫（http://www.geneontology.org/）各個條目進行映射，計算其數目，用Benjamini 法[24]對P值進行校正后，Q＜0.05的GO條目即為在DEGs中顯著富集的GO條目。通過與KEGG （Kyoto encyclopedia of genes and genomes，京都基因與基因組百科全書）數據庫（http://wego.genomics.org.cn）進行比對，對基因涉及的信號通路或代謝途徑（pathway）進行分析。

1.5 可變剪接與單核苷酸多態性的分析及新轉錄本的預測

利用ASprofile軟件的RPKM工具分析獲得可變剪接事件結構詳情和表達量；用Samtools-0.1.19進行單核苷酸多態性（single nucleotide polymorphysm，SNP）分析[25]。利用 Cufflinks軟件v2.2.1（http://cufflinks. cbcb.umd.edu/ howitworks.html）將比對上基因組的測序序列進行組裝拼接[26]。經過濾掉低質量序列（長度≤180 bp，Q 值≤10）后，將組裝的轉錄本序列與牦牛基因組上的基因注釋信息進行比較，如組裝的轉錄本序列未與現有基因比對上，而是位于現有基因之間的基因組上，同時滿足以下條件: 距離現有的注釋基因 200 bp 以上、長度不短于180 bp、測序深度不小于2，則這些序列確定為潛在的新轉錄本及新基因[27]。

1.6 實時熒光定量PCR驗證

為了進一步驗證測序結果的準確性，以為內參基因[28]，挑選5個基因，采用實時熒光定量PCR （qRT-PCR）方法驗證基因表達量。用Primer Premier 5軟件設計基因定量引物，基因名稱及引物信息見表1。以1.3中合成的cDNA第一鏈為模板，擴增目的基因。qRT-PCR反應體系為10 μL：上、下游引物各0.8 μL，Sso AdvancedTM SYBR?Green Super mix 5 μL，ddH2O 2.9 μL，cDNA模板0.5 μL。RT-qPCR擴增程序: 95℃預變性3 min；95℃變性10 s，60℃（根據實際引物退火溫度進行調整）退火20 s，72℃延伸60 s，30個循環；最后72℃延伸5 min；4 ℃保存。采用2-ΔΔCt計算得到基因的相對表達量。

表1 實時熒光定量PCR引物信息

2 結果

2.1 測序質量評估及基本數據分析

牦牛IVF后的卵裂率和囊胚率分別為69.3%和26.2%。經檢測，本研究中構建的牦牛2-細胞、4-細胞、8-細胞、桑椹胚和囊胚等5個發育階段的胚胎微量文庫滿足轉錄組測序要求。利用Agilent High Sensitivity DNA Kit試劑盒經Agilent 2100 Bioanalyzer檢測，5個樣本峰圖顯示在1―2 kb片段長度范圍存在明顯的目的產物主峰，有部分1 kb以下小片段產物但比例較小，表明原始樣本完整性較好，判定合格，符合建庫要求。5個樣本Q30的百分比為90.27% ― 93.18%，說明測序質量和文庫構建質量高，測序數據準確可靠，可滿足后續分析。測序結果中5個樣本堿基A-T、C-G含量都基本對應重合，說明堿基組成穩定平衡，測序質量高。從堿基質量分布圖可知，5個樣本堿基質量穩定在30%―40%之間，低質量堿基比例小，說明測序質量較好。根據飽和量分析圖可知，5個樣本均得到了較高的基因表達量。

2.2 測序結果比對分析

對HiSeqTM2500測序所得的原始測序序列進行數據過濾后，得到2-細胞、4-細胞、8-細胞、桑椹胚和囊胚等5個發育階段的牦牛胚胎過濾后測序序列為47355570―50 855 888條。采用TopHat軟件將獲得的Clean reads與參考基因組進行比對分析[20,21]，結果表明，每個階段有85.65%―90.02% Clean reads比對上牦牛參考基因，比對到基因組多個位置的序列比例（multi_map rate）為3.98%―4.93%，符合要求（表2）。各發育階段比對上的轉錄本和預測新轉錄本的數量見表3，其中，8-細胞期比對上的轉錄本最多（14 893個），而囊胚比對上的轉錄本數量最少（9 827）。利用ASprofile軟件可變剪接分析表明，牦牛胚胎主要有5大類剪接事件：（1）外顯子跳躍（skipped exon, SKIP）和盒式外顯子跳躍（multi-SKIP, MSKIP），（2）內含子滯留（intron retention, IR）和多重內含子滯留（multi-IR, MIR），（3）可變5′端或3′端剪接（alternative exon, AE），（4）轉錄起始區域可變剪接（transcription start site, TSS），（5）轉錄結束區域可變剪接（transcription terminal site, TTS），其中，TSS和TTS所占比例最大（圖1）。利用Samtools-0.1.19進行SNP分析顯示，從2-細胞發育到囊胚，牦牛胚胎在每個發育階段有70841―234131個位點存在單核苷酸多態性（圖2）。

表2 牦牛早期胚胎序列與牦牛基因組比對統計表

SKIP：單外顯子跳躍；MSKIP：多外顯子跳躍；IR：單內含子保留；MIR：多內含子保留；AE：可變5’或3’端剪切；TSS：轉錄起始區域可變剪切；TTS：轉錄結束區域可變剪切；XSKIP：邊界模糊型單外顯子跳躍；XMSKIP：邊界模糊型多外顯子跳躍；XIR：邊界模糊型單內含子保留；XMIR：邊界模糊型多內含子保留；XAE：邊界模糊型5′或3′端可變剪切

2.3 早期胚胎的基因表達特征及差異表達基因分析

在基因表達水平的reads值的基礎上，利用RPKM法[22]計算得到基因在胚胎不同發育階段的RPKM值，4-細胞開始表達的基因有1 221個，8-細胞開始表達的基因有1 116個，而在桑椹胚和囊胚開始表達的基因分別只有142個和564個。隨著胚胎發育的進行，、、、、和等母源基因的表達量逐漸減少（圖3），而、、、、、和等基因在4-細胞期的表達量開始增加。和分別在8-細胞和囊胚才開始大量表達。

表3 檢測到的轉錄本及新轉錄本的數量

2：2-細胞；4：4-細胞；8：8-細胞；M：桑椹胚；B：囊胚。下同

圖3 牦牛早期胚胎發育過程中部分母源基因表達量變化圖

采用DEGseq法分析牦牛胚胎發育階段DEG的結果顯示，從2-細胞到4-細胞、4-細胞到8-細胞、8-細胞到桑椹胚及桑椹胚到囊胚4個連續發育階段的DEGs數分別為6 922、7 601、8 071和10 555個。除8-細胞發育到桑椹胚的上調DEGs只有2 349個以外，其余3個發育階段的上調DEGs數都在4 100個以上，而從2-細胞發育到囊胚期的每個階段下調DEGs數都在不斷增加（圖4）。

為了驗證RNA-Seq所得結果可信性，從高通量測序結果中選取、、、和等5個DEGs，采用qRT-PCR技術，隨機在2-、4-、8-細胞、桑椹胚和囊胚等5個階檢測其表達情況，結果表明: qRT-PCR結果與RNA-Seq結果基本一致（圖5）。

2.4 差異表達基因的GO富集與KEGG通路分析

GO分析顯示，從2-細胞到4-細胞期、4-細胞到8-細胞期、8-細胞期到桑椹胚及桑椹胚到囊胚4個發育階段分別有6 473、7 133、7 577和9 930個DEGs得到歸類注釋，都涉及生物過程（Biological process，BP）、細胞組分（cellular component，CC）和分子功能（molecular function，MF）3大類62個二級條目（圖6）。在BP分類中有23個二級條目中占比例最大依次為細胞過程（cellular process）、單有機體過程（single-organism process）、生物調節（biological regulation）和代謝過程（metabolic process）。在CC分類中有18個二級條目占比例最高依次是細胞部分（cell part）占比例最多，其次為細胞器（organelle）及細胞器部件（organelle part）。在MF分類中有21個二級條目，占比例最大的是綁定分子（binding），其次催化活性（catalytic activity）及分子傳感器活性（molecular transducer activity）。不同發育階段的GO二級條目排列順序有一定的特異性，例如，從桑椹胚發育到囊胚的過程中MF類別的成形素（morphogen）條目上調。

2 vs 4：2-細胞vs 4-細胞；4 vs 8：4-細胞vs 8-細胞；8 vs M：8-細胞vs桑椹胚；M vs B：桑椹胚vs囊胚

X-軸表示基因表達量，Y-軸表示胚胎發育階段X-axis denotes gene expression level, Y-axis denotes development stages

BP1：細胞過程；BP2：單有機體過程；BP3：生物調節；BP4：代謝過程；BP5：細胞組成的組織或合成；BP6：發育過程；BP7：定位；BP8：刺激反應；BP9：多細胞生物過程；BP10：免疫系統過程；BP11：繁殖過程；BP12：移動；BP13：生物粘附；BP14：多生物體過程；BP15：行為；BP16：信號；BP17：生長；BP18：節律過程；BP19：激素分泌；BP20：生物相；BP21：細胞聚集；BP22：細胞殺傷；BP23：繁殖；CC1：細胞部分；CC2：細胞器；CC3：細胞器部分；CC4：膜；CC5：膜部分；CC6：高分子復合物；CC7：細胞外區域部分；CC8：細胞連接；CC9：胞外區；CC10：膜包圍腔；CC11：突出部分；CC12：突觸； CC13：細胞外基質；CC14：細胞Cell；CC15：膠原三聚體；CC16：細胞外基質部分；CC17：細胞核；CC18：病毒體部分；MF1：結合；MF2：催化；MF3：分子功能的調節；MF4：核酸結合的轉錄因子；MF5：分子傳感器；MF6：運輸；MF7：酶的調節；MF8：蛋白結合的轉錄因子；MF9：結構分子；MF10：鳥嘌呤核苷酸交換因子；MF11：通道調節；MF12：電子載體；MF13：抗氧化；MF14：翻譯調控；MF15：受體調控；MF16：化學引誘物；MF17：化學排斥物；MF18：蛋白標簽；MF19：金屬伴侶；MF20：成形素；MF21：營養庫

牦牛胚胎發育過程中差異表達基因KEGG分析表明，2-細胞到4-細胞期涉及到308條通路，4-細胞到8-細胞期有310條通路，8-細胞期到桑椹胚有316條通路，桑椹胚到囊胚有315條通路，各發育階段的富集前5條通路如表4所示，每個階段的通路種類及富集性都有一定差異，2-細胞到4-細胞階段及桑椹胚到囊胚階段剪接體通路（spliceosome）最為富集，而4-細胞到8-細胞及8-細胞到桑椹胚階段嗅覺轉導通路（olfactory transduction）最為富集。2-細胞到4-細胞、4-細胞到8-細胞、8-細胞到桑椹胚和桑椹胚到囊胚4個階段分別有11、9、10和2個通路顯著富集。

3 討論

新一代高通量測序技術的不斷發展，已徹底改變了轉錄組學的研究，使RNA-Seq無需預先設計探針即可對特定條件下任何生物生長發育階段整體轉錄活動進行測序，準確探測到各種條件下的基因表達情況，發現了許多未知的分子調控機制[5-7]。本研究首次利用RNA-Seq技術從轉錄組學角度揭示牦牛早期胚胎發育機制，為完善牦牛胚胎體外生產提供新思路，同時也為進一步完善牦牛基因結構信息和胚胎發育相關的新基因提供理論基礎。由于哺乳動物2-細胞到囊胚期每個胚胎的總RNA只有200―2 000 pg[8]，這樣微量RNA不能滿足轉錄組測序文庫的構建及高通量測序的基本要求，因此，分別提取2-細胞、4-細胞、8-細胞、桑椹胚和囊胚等5個發育階段的牦牛胚胎總RNA，使用先進Smart-Seq2擴增技術對樣本進行富集并構建測序文庫[9]，再應用RNA高通量測序技術對其進行高通量測序分析。測序質量評估、數據分析及qPCR驗證結果都表明，測序質量和文庫構建質量高，測序數據準確可靠。

表4 牦牛胚胎發育過程中差異表達基因富集前5（Top 5）KEGG通路表

3.1 牦牛早期胚胎發育錄組特點

對HiSeqTM2500測序所得的原始測序序列進行數據過濾后，得到2-細胞、4-細胞、8-細胞、桑椹胚和囊胚5個發育階段的牦牛早期胚胎過濾后測序序列為47 355 570―50 855 888條，比對分析顯示，每個階段有85.65%―90.02%測序序列比對上參考基因的序列（表2）。測序序列數與牦牛基因組比對顯示，8-細胞期比對上的轉錄本最多（14 893），而囊胚期比對上的轉錄本數量最少（9 827）（表3），可能與牦牛IVF囊胚率低（26.2%）有關，說明牦牛胚胎的體外培養系統有待改進。其余階段的轉錄本數與體內受精的普通牛胚胎的轉錄本數基本一致[29-30]。

牛胚胎的SNP主要涉及參與胞外配體信號（和）、內吞作用和胞吐作用（）、凋亡調控（）、細胞應激保護（）、能量代謝（）、蛋白互作（）和轉錄調控（）等途徑影響胚胎發育能力[31]。本研究SNP分析顯示，牦牛2-細胞、4-細胞、8-細胞、桑椹胚和囊胚基因組分別有116 601、234 131、196 420、70 841和94 840個位點存在SNP，說明基因組SNP在牦牛胚胎發育中起重要調控作用。可變剪接是提高轉錄組和蛋白組復雜性的一個重要機制。與Yan等[10]利用單細胞RNA-Seq技術對人的單卵裂球轉錄組測序分析結果一樣，本研究發現牦牛2-細胞、4-細胞、8-細胞、桑椹胚和囊胚中存在大量的可變剪接，說明可變剪接在真核生物中普遍存在,并且在生物體胚胎發育的不同階段中，基因的剪接方式是不斷變化的，以此調控細胞的增殖、分化、遷移和凋亡。

有些基因在胚胎發育的特異階段表達，在特異階段發揮關鍵的調控作用[13,19,29-30]。NANOG能調控合子基因組的激活，缺乏NANOG的胚胎合子基因組激活率低，發育受阻[32]。本研究發現牦牛胚胎的在8-細胞期開始表達。CLDN4是Claudin蛋白家族的成員之一，在縫隙連接（gap junction）形成中起不可或缺的作用。研究證明，婦女黃體期子宮內膜mRNA的表達量與妊娠相關[33]。在小鼠胚胎培養液中添加CLDN4抑制因子會導致胚胎不能形成正常的囊胚腔[34]。在乳腺癌細胞活性研究中還發現CLDN4具有促進細胞的組織侵入作用[35]。因此，與在人的囊胚中檢測結果類似[36]，本研究發現在牦牛囊胚中的表達量顯著提高，說明CLDN4可能在牦牛囊胚發育和附植過程中胚胎滋養層細胞逐漸侵入子宮上皮及基質層都發揮重要作用。

3.2 胚胎基因組的激活

胚胎發育是遺傳信息按一定的時間、空間和次序表達的結果，即按照發育的遺傳程序（genetic program）展開的結果。早期胚胎的發育屬母型調控，即由來自卵母細胞發生及成熟期間合成的大量mRNAs和蛋白質來調控，母源mRNA 在胚胎發育早期起著重要生理作用[37]。隨著發育的進行，母源mRNA和蛋白質逐漸降解，而胚胎基因組激活（embryonic genome activation，EGA）啟動，發育從母型調控向胚胎型調控的過渡（maternal- to-embryonic transition，MET）。不同的物種MET發生的時期不同，小鼠MET發生在2-細胞期[38]，人[39]和豬[38]的發生在4-細胞到8-細胞期，牛胚胎基因組的主要轉錄開始于8-細胞到16-細胞期[12,38]。、、、、和等基因被證明為標志性的母源基因[29-30,39]，隨著胚胎發育的進行，牦牛的這些母源基因表達量逐漸減少（圖3）。在牦牛胚胎4-和8-細胞期開始表達的基因分別有1 221個和1 116個，而在桑椹胚和囊胚期開始表達的基因分別有142個和564個。另外，牦牛胚胎在2-細胞到4-細胞期以及4-細胞到8-細胞期即出現大量差異表達基因，但是EGA的標志性基因和[29-30,33]分別在8-細胞期和囊胚期才開始大量表達。因此，綜上分析筆者認為牦牛EGA可能發生在4-細胞到8-細胞期。

3.3 胚胎差異表達基因的GO注釋及KEGG分析

從2-細胞到4-細胞、4-細胞到8-細胞、8-細胞到桑椹胚及桑椹胚到囊胚4個發育階段的DEG數分別為6 922、7 601、8 071和10 555個（圖4），說明不同發育階段的牦牛胚胎的發育調控機制存在明顯的時序性差異。GO分析表明，從2-細胞到4-細胞、4-細胞到8-細胞、8-細胞到桑椹胚及桑椹胚到囊胚4個發育階段的DEGs歸類注釋都涉及生物過程（BP）、細胞組分（CC）和分子功能（MF）3大類62個二級條目（圖6）。4-細胞到8-細胞期的BP類別的發育程序（developmental process）條目上調和桑椹胚到囊胚的MF類別的成形素條目上調。普通牛在早期胚胎發育過程中，GO條目也在4-細胞到8-細胞期發生較大改變[29,30]。多細胞動物形態形成中，形成素具有給予細胞位置信號的作用，從而決定胚胎細胞形成不同組織、器官和構成有序空間結構的圖式形成（pattern formation）以及胚胎發育的反應速度和擴散速度[40]。牦牛桑椹胚到囊胚階段的成形素條目上調，可能導致具有內細胞團和滋養層的囊胚形成，從而決定胚胎細胞的發育結果。

本研究發現，牦牛早期胚胎發育過程中DEGs富集的主要通路也有時序特異性，例如，在2-細胞到4-細胞富集前3的通路由高到低依次是剪接體、RNA轉運和泛素介導的蛋白水解，而在4-細胞到8-細胞富集前3的通路則是嗅覺轉導、神經活性配體-受體互作和核苷酸切除修復（表4）。但是就總體而言，剪接體、神經活性配體-受體互作、細胞因子及其受體相互作用、泛素介導的蛋白水解、RNA轉運等通路在早期胚胎發育的各個階段都是富集通路，這與在其它動物的發現基本一致[29,41-42]。在牦牛4-細胞胚胎的剪接體通路和內吞作用通路即被激活。研究證明，在早期胚胎發育過程中剪接體通路和內吞作用通路與胚胎母源RNA降解和基因組的啟動有關[41,43-44]。細胞因子及其受體相互作用通路在早期胚胎的發育和附植過程中發揮重要作用[41,45]，胚胎附植反復失敗的婦女子宮內膜中有許多mRNA表達異常，其中表達下調的多數mRNA都涉及細胞因子及其受體相互作用通路[46]。泛素是一種在細胞內廣泛分布的高度保守的小蛋白，在DNA修復、蛋白質降解的標記（泛素化）、蛋白質的合成與轉運、基因轉錄調控及信號轉導等各個生命活動中發揮著重要的作用[47]。有研究提示，泛素介導的蛋白水解通路和ErbB信號通路的正常調節可能是保證胚胎正常發育的基本條件[48]。妊娠期的飲酒惡習會導致胎兒發育的缺陷，目前的研究認為這種胚胎發育缺陷很可能是由于酒精影響妊娠期JAK-STAT信號通路、神經活性配體-受體互作、Toll樣受體（TLR）信號通路、細胞因子及其受體相互作用等通路引起[49]。因此，在牦牛早期胚胎發育過程中，上述這些富集通路保證了胚胎的正常發育。嗅覺轉導通路4-細胞到8-細胞以及8-細胞到桑椹胚都為顯著富集的通路，但其對胚胎發育的調控作用尚不清楚。

4 結論

本研究首次利用RNA-Seq技術對牦牛體外受精胚胎發育不同階段轉錄組進行了分析，獲得了眾多差異基因和有關通路的富集。不同階段差異表達基因在數量、功能、分類和代謝通路都有各自的特異性。對于完善牦牛胚胎體外生產技術，及進一步完善牦牛基因結構信息和胚胎發育相關的新基因提供了理論基礎。

[1] Wiener G, Han J L, Long R J.. 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. Reproduction in female yaks () and opportunities for improvement., 2003, 59: 1303-1312.

[3] Zi X D, Yin R H, Chen S W, Liang G N, Zhang D W, Guo C H. Developmental competence of embryos derived from reciprocalfertilization between yak () and cattle ()., 2009, 55(5): 480-483.

[4] Xiao X, Zi X D, Niu H R, Xiong X R, Zhong J C, Li J, Wang L, Wang Y. Effect of addition of FSH, LH and proteasome inhibitor MG132 tomaturation medium on the developmental competence of yak () oocytes., 2014, 12: 30.

[5] Mortazavi A, Williams B A, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq., 2008, 5(7): 621-628.

[6] 朱志明, 陳紅萍, 林如龍, 繆中緯, 辛清武, 李麗, 張丹青, 鄭嫩珠. 山麻鴨開產期和產蛋高峰期卵巢組織轉錄組分析. 中國農業科學, 2016, 49(5): 998-1007.

ZHU Z M, CHEN H P, LIN R L, MIAO Z W, XIN Q W, LI L, ZHANG D Q, ZHENG N Z. Transcriptome analysis of ovary tissue in early laying period and egg laying peak period of Shanma ducks., 2016, 49(5): 998-1007. (in Chinese)

[7] 王文龍, 馮陳晨, 紅梅, 岳建偉, 呼和巴特爾, 劉春霞. 不同發育階段斯氏副柔線蟲比較轉錄組學分析. 中國農業科學, 2017, 50(23): 4644-4655.

WANG W L, FENG C C, HONG M, YUE J W, HUHE B T E, LIU C X. The comparative transcriptome analysis ofat different developmental stages., 2017, 50(23): 4644-4655. (in Chinese)

[8] Gilbert I, Scantland S, Sylvestre E L, Gravel C, Laflamme I, Sirard MA, Robert C. The dynamics of gene products fluctuation during bovine pre-hatching development., 2009, 76(8): 762-772.

[9] Picelli S, Bjorklund A K, Faridani O R, Sagasser S, Winberg G, Sandberg R. Smart-seq2 for sensitive full-length transcriptome profiling in single cells., 2013, 10(11): 1096?1098.

[10] Yan L Y, Yang M Y, Guo H S, Yang L, Wu J, Li R, Liu P, Lian Y, Zheng X, Yan J, Huang J, Li M, Wu X, Wen L, Lao K, Li R, Qiao J, Tang F. Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells., 2013, 20(9): 1131-1139.

[11] Cao S, Han J, Wu J, Li Q, Liu S, Zhang W, Pei Y, Ruan X, Liu Z, Wang X, Lim B, Li N. Specific gene-regulation networks during the pre-implantation development of the pig embryo as revealed by deep sequencing., 2014, 15:4.

[12] Graf A, Krebs S, Heininen-Brown M, Zakhartchenko V, Blum H, Wolf E. Genome activation in bovine embryos: review of the literature and new insights from RNA sequencing experiments., 2014, 149(1-2): 46-58.

[13] Jiang Z, Sun J, Dong H, Luo O, Zheng X, Obergfell C, Tang Y, Bi J, O'Neill R, Ruan Y, Chen J, Tian XC. Transcriptional profiles of bovinepreimplantation development., 2014, 15: 756.

[14] 于學穎, 郭芹芹, 郝海生, 孫尉峻, 趙學明, 朱化彬, 楊凌, 杜衛華. 谷胱甘肽促進牛體外受精胚胎發育的轉錄組初探. 畜牧獸醫學報, 2016, 47 (7): 1363-1372.

YU X Y, GUO Q Q, HAO H S, SUN W J, ZHAO X M, ZHU H B, YANG L, DU W H. Transcriptome of bovine IVF Embryos treated with glutathione., 2016, 47 (7): 1363-1372. (in Chinese)

[15] 蘭道亮, 熊顯榮, 位艷麗, 徐通, 鐘金城, 字向東, 王永, 李鍵. 基于RNA-Seq 高通量測序技術的牦牛卵巢轉錄組研究: 進一步完善牦牛基因結構及挖掘與繁殖相關新基因. 中國科學: 生命科學, 2014, 44(3): 307-317.

LAN D L, XIONG XR, WEI Y L, XU T, ZHONG J C, ZI X D, WANG Y, LI J. RNA-Seq analysis of yak ovary: improving yak gene structure information and mining reproduction-related genes., 2014, 44(3): 307-317. (in Chinese)

[16] 蘭道亮, 熊顯榮, 陳亞冰, 澤讓東科, 艾鷖, 李鍵. 牦牛體外成熟卵母細胞差異轉錄組學研究. 中國科學(生命科學), 2017, 47(10): 1099-1112.

LAN D L, XIONG X R, CHEN Y B, ZERANG D K, AI Y, LI J. Differential transcriptome analysis of yak oocytesmaturation., 2017, 47(10): 1099-1112. (in Chinese)

[17] 鄭杰, 蒲思穎, 楊遠瀟, 王琴, 楊繞芬, 字向東. 基于高通量測序的犏牛囊胚玻璃化冷凍損傷機制研究. 2017, 48 (10): 1871-1881.

ZHENG J, PU S Y, YANG Y X, WANG Q, YANG R F, ZI X D. Exploring mechanism for vitrification damage of the cross-bred blastocysts of the yak via high-throughput sequencing., 2017, 48 (10): 1871-1881. (in Chinese)

[18] 曾賢彬, 柴志欣, 王永, 馬志杰, 楊琴, 宋喬喬, 鐘金城. 犏牛精子發生阻滯的比較轉錄組研究. 中國科學(生命科學), 2014, 44(6): 584-601.

ZENG X B, CHAI Z X, WANG Y, MA Z J, YANG Q, SONG Q Q, ZHONG J C. Comparative transcriptome analysis of spermatogenesis arrest in cattle-yak., 2014, 44(6): 584-601. (in Chinese)

[19] 張春強, 趙德超, 韓志玲, 張文廣, 張家新. 父源性印記基因在牛早期胚胎中的表達. 中國農業科學, 2013, 46(18): 3887-3893.

ZHANG C Q, ZHAO D C, HAN Z L, ZHANG W G, ZHANG J X. The expression of paternal impringting genes in bovine early stage embryo., 2013, 46(18): 3887-3893. (in Chinese)

[20] Qiu Q, Zhang G, Ma T, Qian W, Wang J, Ye Z, Cao C, Hu Q, Kim J, Larkin DM, Auvil L, Capitanu B, Ma J, Lewin H A, Qian X, Lang Y, Zhou R, Wang L, Wang K, Xia J, Liao S, Pan S, Lu X, Hou H, Wang Y, Zang X, Yin Y, Ma H, Zhang J, Wang Z, Zhang Y, Zhang D, Yonezawa T, Hasegawa M, Zhong Y, Liu W, Zhang Y, Huang Z, Zhang S, Long R, Yang H, Wang J, Lenstra J A, Cooper D N, Wu Y, Wang J, Shi P, Wang J, Liu J. The yak genome and adaptation to life at high altitude., 2012, 44(8): 946-949.

[21] Trapnell C, Pachter L, Salzberg S L. TopHat: discovering splice junctions with RNA-seq., 2009, 25(9): 1105?1111.

[22] Wanger G P, Kin K, Lynch V J. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples., 2012, 131(4): 281?285.

[23] Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data., 2010, 26(1): 136?138.

[24] Benjamini Y, Hochberg Y. Controlling the false discovery rate - a practical and powerful approach to multiple testing., 1995, 57(1): 289-300.

[25] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and SAMtools., 2009, 25(16): 2078-2079.

[26] Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, van Baren M J, Salzberg S L, Wold B J, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation., 2010, 28(5): 511?515.

[27] Roberts A, Pimentel H, Trapnell C, Pachter L. Identification of novel transcripts in annotated genomes using RNA-Seq., 2011, 27(17): 2325-2329.

[28] Robert C, McGraw S, Massicotte L, Pravetoni M, Gandolfi F, Sirard MA. Quantification of housekeeping transcript levels during the development of bovine preimplantation embryos., 2002, 67(5): 1465-1472.

[29] Kues W A, Sudheer S, Herrmann D, Carnwath J W, Havlicek V, Besenfelder U, Lehrach H, Adjaye J, Niemann H. Genome-wide expression profiling reveals distinct clusters of transcriptional regulation during bovine preimplantation development., 2008, 105(50): 19768-19773.

[30] Graf A, Krebs S, Zakhartchenko V, Schwalb B, Blum H, Wolf E. Fine mapping of genome activation in bovine embryos by RNA sequencing., 2014, 111(11): 4139-4144.

[31] Ortega M S, Kurian J J, McKenna R, Hansen P J. Characteristics of candidate genes associated with embryonic development in the cow: evidence for a role for WBP1 in development to the blastocyst stage., 2017, 12(5): e0178041.

[32] Lee M T, Bonneau A R, Takacs C M, Bazzini A A, DiVito K R, Fleming E S, Giraldez A J. Nanog, Pou5f1 and SoxB1 activate zygotic gene expression during the maternal-to-zygotic transition., 2013, 503(7476): 360-364.

[33] Serafini P C, Silva I D, Smith G D, Motta E L, Rocha A M, Baracat E C. Endometrial claudin-4 and leukemia inhibitory factor are associated with assisted reproduction outcome., 2009, 7: 30.

[34] Moriwaki K, Tsukita S, Furuse M. Tight junctions containing claudin 4 and 6 are essential for blastocyst formation in preimplantation mouse embryos., 2007, 312(2): 509-522.

[35] Webb P G, Spillman M A, Baumgartner H K. Claudins play a role in normal and tumor cell motility.,2013, 14: 19.

[36] Munch E M, Sparks A E, Gonzalez Bosquet J, Christenson L K, Devor E J, Van Voorhis B J. Differentially expressed genes in preimplantation human embryos: potential candidate genes for blastocyst formation and implantation., 2016, 33(8): 1017-1025.

[37] Tadros W, Lipshitz H D. The maternal-to-zygotic transition: a play in two acts., 2009, 136(18): 3033-3042.

[38] Sirard M A. Factors affecting oocyte and embryo transcriptomes., 2012, 47(Suppl. 4): 148-155.

[39] Braude P, Bolton V, Moore S. Human gene expression first occurs between the four- and eight-cell stages of preimplantation development., 1988, 332(6163): 459-461.

[40] Briscoe J, Small S. Morphogen rules: design principles of gradient-mediated embryo patterning., 2015, 142(23): 3996-4009.

[41] Zuo Y, Su G, Wang S, Yang L, Liao M, Wei Z, Bai C, Li G. Exploring timing activation of functional pathway based on differential co-expression analysis in preimplantation embryogenesis., 2016, 7(45): 74120-74131.

[42] Rauwerda H, Pagano J F B, de Leeuw W C, Ensink W, Nehrdich U, de Jong M, Jonker M, Spaink H P, Breit T M. Transcriptome dynamics in early zebrafish embryogenesis determined by high-resolution time course analysis of 180 successive, individual zebrafish embryos., 2017, 18(1): 287.

[43] Sato M, Sato K. Dynamic regulation of autophagy and endocytosis for cell remodeling during early development., 2013, 14(5): 479-486.

[44] Abada A, Elazar Z. Getting ready for building: signaling and autophagosome biogenesis., 2014, 15(8): 839-852.

[45] Altm?e S, Reimand J, Hovatta O, Zhang P, Kere J, Laisk T, Saare M, Peters M, Vilo J, Stavreus-Evers A, Salumets A. Research resource: interactome of human embryo implantation: identification of gene expression pathways, regulation, and integrated regulatory networks., 2012, 26(1): 203-217.

[46] Shi C, Han H J, Fan L J, Guan J, Zheng X B, Chen X, Liang R, Zhang X W, Sun K K, Cui Q H, Shen H. Diverse endometrial mRNA signatures during the window of implantation in patients with repeated implantation failure.,Doi: 10.1080/14647273.2017.1324180.

[47] Hospenthal M K, Freund S M V, Komander D. Assembly, analysis and architecture of atypical ubiquitin chains., 2013, 20(5): 555-565.

[48] Isom S C, Stevens J R, Li R, Spollen W G, Cox L, Spate L D, Murphy C N, Prather R S. Transcriptional profiling by RNA-Seq of peri-attachment porcine embryos generated by a variety of assisted reproductive technologies., 2013, 45(14): 577-589.

[49] Kim Y Y, Roubal I, Lee Y S, Kim J S, Hoang M, Mathiyakom N, Kim Y. Alcohol-induced molecular dysregulation in human embryonic stem cell-derived neural precursor cells., 2016, 11(9): e0163812.

Transcriptomic Analysis of IVF Embryonic Development in the Yak () Via RNA-Seq

ZI XiangDong1, LUO Bin1, XIA Wei1, ZHENG YuCai1, XIONG XianRong1, LI Jian1, ZHONG JinCheng2, ZHU JiangJiang2, ZHANG ZhengFan1

(1College of Life Science and Technology, Southwest Minzu University, Chengdu 610041;2Institute of Qinghai-Tibetan Plateau, Southwest Minzu University, Chengdu 610041)

【Objective】The objectives of this study were to investigate the transcriptome differences and identify function, classification and metabolic pathways of the differentially expressed genes (DEG) at different developmental stages of yak embryos derived fromfertilization (IVF), which are necessary to better understand the mechanism that regulates embryonic development and provide theoretical basis for improvingembryo production in the yak (). 【Method】Total RNA was extracted from IVF derived yak embryos at 2-cell, 4-cell, 8-cell, morula and blastocyst stages and amplified via the Smart-seq2 method, and the constructed RNA libraries were sequenced using the HiSeqTM2500 high-throughput sequencing method. 【Result】After IVF, the average cleavage rates and blastocyst rates were 69.3% and 26.2%, respectively.A total of 47 355 570 to 50 855 888 clean reads were obtained from 2-cell, 4-cell, 8-cell, morula and blastocyst stages, respectively, of which, 85.65% to 90.02% were covered in the yak reference genome. In total, the number of transcripts mapped to yak genomewas highest for 8-cell (14 893) and lowest for blastocysts (9 827). The transcripts mainly had five patterns of alternative splice, of which, the two largest proportions were transcription start site (TSS) and transcription terminal site (TTS). The SNP numbers of the five stages of yak embryonic transcripts were 116 601, 234 131, 196 420, 70 841 and 94 840, respectively. A total of 1 221, 1 116, 142 and 564 transcripts were first detected at the 4-cell, 8-cell, morula and blastocyst stages, respectively. As embryo development proceeded, maternally derived transcripts such as,,,,andetc.were decreased, whereas embryonic transcripts such as,,,,,andetc. were increased at specific stages. When |log2ratio| ≥1 and Q-value＜0.05 were set as thresholds for identifying DEGs, a total of 6 922, 7 601, 8 071 and 10 555 DEGs were identified from 2-cell vs. 4 cell, 4-cell vs. 8-cell, 8-cell vs. morula, and morula vs. blastocyst, respectively. The GO distributions of the DEGs were classified into three categories: biological processes (BP), cellular components (CC), molecular functions (MF) with a total of 62 subcategories of two successive stages. KEGG enrichment analysis of DEGs showed that DEGs of 2-cell vs. 4-cell participated in 308 pathways, and significantly enriched in 11 pathways such as spliceosome, RNA transport and ubiquitin mediated proteolysis etc. DEGs of 4-cell vs. 8-cell participated in 310 pathways, and significantly enriched in 9 pathways such as olfactory transduction, neuroactive ligand-receptor interaction and nucleotide excision repair etc. DEGs of 8-cell vs. morula participated in 316 pathways, and significantly enriched in 10 pathways such as olfactory transduction, ubiquitin mediated proteolysis and neuroactive ligand-receptor interaction etc. DEGs of morula vs. blastocyst participated in 315 pathways, significantly enriched in 2 pathways i.e., spliceosome and RNA transport.【Conclusion】This is the first study for analyzing the transcriptomes of IVF derived yak-embryos at different stages using high-throughput sequencing. A number of DEGs and their function, classification and metabolic pathways were discovered, which enriched transcriptome information for yak embryos. In addition, the results provided a foundation and reference to uncoverthe mechanism that regulates embryonic development and improvesembryo production of the yak species.

yak;fertilization; embryo; transcriptome; RNA-Seq

（責任編輯 林鑒非）

10.3864/j.issn.0578-1752.2018.08.015

2017-07-04；

2018-02-02

國家“973”前期專項（2007CB116204）、中央高校基本科研業務費專項資金（2015NZYTD02）

字向東，Tel：028-85528868；E-mail：zixd2000@sina.com