基于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