0 引言
【研究意义】牦牛（Bos grunniens）是分布在海拔2 000―5 000 m以青藏高原及其毗邻高山、亚高山地区的珍奇稀有牛种。目前,全世界有约1 400余万头牦牛,其中,中国拥有的牦牛数量占世界95%以上[1]。牦牛是我国青藏高原人民赖以生存的生产和生活资料,对高寒低氧的独特适应能力决定了该畜种在青藏高原畜牧业中占据不可替代的特殊地位。牦牛的繁殖率低（40%―60%）[1,2],体外受精（in vitro fertilization, IVF）的囊胚发育率也不高（10%―30%）[3,4]。因此,研究牦牛早期胚胎发育调控机理,对提高牦牛繁殖效率及完善牦牛胚胎体外生产（in vitro production, IVP）技术具有重要意义。【前人研究进展】转录组（transcriptome）是特定组织或细胞在某一环境或生理条件下所转录表达的所有 RNA 的总和,既包括编码蛋白的 mRNA,也包括非编码RNA,是连接基因组遗传信息与蛋白质组生物功能的纽带。因此,转录组测序技术（RNA sequencing,RNA-Seq）为研究基因表达及调控提供了重要的手段和方法,得到广泛使用[5,6,7]。虽然卵母细胞、早期胚胎的 RNA含量低[8],难以达到转录组测序建库的RNA最小起始量要求,但是,随着Smart-seq2等微量转录组学测序技术的建立与发展[9],研究者已成功利用RNA-Seq技术解析了人[10]、猪[11]和牛[12,13,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-1 FSH、50 IU·mL-1 LH和1 μg·mL-1 17β-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与胚胎体外培养（in vitro 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验证
为了进一步验证测序结果的准确性,以H2A为内参基因[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计算得到基因的相对表达量。Table 1
表1
表1实时荧光定量PCR引物信息
Table 1The qRT-PCR primers
| 基因名称Gene | 上游引物Forward primer (5′-3′) | 下游引物（5′→3′）Reverse primer (5′-3′) | 产物长度Product size（bp） |
|---|---|---|---|
| H2A | GCGTATTACCCCTCGTCACT | CTTCTGTTGTCCTTTCTTTCC | 139 |
| CXCL-9 | AGCCTGCCTAAACCCAGATT | AGTGGTTGGTGAAGTGGGAG | 187 |
| NANOG | GCTGTGTGTGCTCAATGACA | GCCTGGGTATTCTGCCATTG | 211 |
| CDC-20 | ATATCCTGTCCAGTGGCTCG | TTGATGCTGGGTGAAGGTCT | 230 |
| MELK | CAGAGGCAGATGTTTGGAGC | CCAGGGGTGGTTCAAGAGAT | 224 |
| PCNA | AGTAGACTTTTGGGAGGGGC | ACAGTGAGTTCTACCACCGG | 161 |
新窗口打开
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个发育阶段的牦牛胚胎过滤后测序序列为47 355 570―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-细胞发育到囊胚,牦牛胚胎在每个发育阶段有70 841―234 131个位点存在单核苷酸多态性（图2）。Table 2
表2
表2牦牛早期胚胎序列与牦牛基因组比对统计表
Table 2RNA-Seq and mapping to the reference genome
| 2-细胞2-cell | 4-细胞4-cell | 8-细胞8-cell | 桑椹胚Morula | 囊胚Blastocyst | |
|---|---|---|---|---|---|
| 干净数据序列数Clean reads | 50855888 | 47978848 | 48831636 | 47355570 | 48455410 |
| 比对上基因组的序列数 Mapped reads | 45777962 | 42998819 | 43541753 | 42557017 | 41500742 |
| 比对上基因组的序列比例Percentage of mapped reads | 90.02 | 89.62 | 89.17 | 89.87 | 85.65 |
| 未比对到基因组的Reads数Unmapped reads | 5077926 | 4980029 | 5289883 | 4798553 | 6954668 |
| 比对到基因组多位点的Reads数Multi-map reads | 2125900 | 1909991 | 2021485 | 1961978 | 2389740 |
| 比对到基因组多位点的Reads比例Percentage of multi-mapped reads | 4.18 | 3.98 | 4.14 | 4.14 | 4.93 |
新窗口打开
Table 3
表3
表3检测到的转录本及新转录本的数量
Table 3Detected transcripts and novel transcripts
| 胚胎 Embryos | 比对上的转录本 No. of transcripts | 新转录本 No. of novel transcripts |
|---|---|---|
| 2-细胞2-cell | 14326 | 583 |
| 4-细胞4-cell | 13791 | 864 |
| 8-细胞8-cell | 14893 | 665 |
| 桑椹胚Morula | 11469 | 595 |
| 囊胚Blastocyst | 9827 | 759 |
新窗口打开
显示原图|下载原图ZIP|生成PPT图1牦牛早期胚胎发育过程中可变剪接事件分析统计图
SKIP：单外显子跳跃;MSKIP：多外显子跳跃;IR：单内含子保留;MIR：多内含子保留;AE：可变5’或3’端剪切;TSS：转录起始区域可变剪切;TTS：转录结束区域可变剪切;XSKIP：边界模糊型单外显子跳跃;XMSKIP：边界模糊型多外显子跳跃;XIR：边界模糊型单内含子保留;XMIR：边界模糊型多内含子保留;XAE：边界模糊型5′或3′端可变剪切
-->Fig. 1Alternative splicing events detected during embryonic development
SKIP: Skipped exon; MSKIP: Multi-SKIP; IR: Intron retention; MIR: Multi-IR; AE: Alternative exon; TSS: Transcription start site; TTS: Transcription terminal site; XSKIP: Approximate SKIP; XMSKIP: Approximate MSKIP; XIR: Approximate IR; XMIR: Approximate MIR; XAE: Approximate AE(5′, 3′, or both)
-->
显示原图|下载原图ZIP|生成PPT图2牦牛早期胚胎发育过程中单核苷酸多态性分析统计图
2：2-细胞;4：4-细胞;8：8-细胞;M：桑椹胚;B：囊胚。下同
-->Fig. 2Putative SNPs detected during embryonic development
2: 2-cell; 4: 4-cell; 8: 8-cell; M: Morula; B: Blastocyst. The same as below
-->
2.3 早期胚胎的基因表达特征及差异表达基因分析
在基因表达水平的reads值的基础上,利用RPKM法[22]计算得到基因在胚胎不同发育阶段的RPKM值,4-细胞开始表达的基因有1 221个,8-细胞开始表达的基因有1 116个,而在桑椹胚和囊胚开始表达的基因分别只有142个和564个。随着胚胎发育的进行,BMP15、KIT、GDF9、STAT3、ZP3和ZP4等母源基因的表达量逐渐减少（图3）,而SARS、IL18、ACO2、TXN2、ATP5B、PCGF4、UBE3A、MAPK13、SNURF和JUP等基因在4-细胞期的表达量开始增加。NANOG和CLDN4分别在8-细胞和囊胚才开始大量表达。
显示原图|下载原图ZIP|生成PPT图3牦牛早期胚胎发育过程中部分母源基因表达量变化图
-->Fig. 3The expression tendency of some maternal genes during embryonic development
-->
采用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）。
显示原图|下载原图ZIP|生成PPT图4牦牛早期胚胎发育过程中的差异表达基因统计图
2 vs 4：2-细胞vs 4-细胞;4 vs 8：4-细胞vs 8-细胞;8 vs M：8-细胞vs桑椹胚;M vs B：桑椹胚vs囊胚
-->Fig. 4DEGs statistic during early embryonic development of the yak
2 vs 4: 2-cell vs 4-cell; 4 vs 8: 4-cell vs 8-cell; 8 vs M: 8-cell vs Morula; M vs B: Morula vs Blastocyst
-->
为了验证RNA-Seq所得结果可信性,从高通量测序结果中选取CDC-20、PCNA、MELK、NANOG和CXCL-9等5个DEGs,采用qRT-PCR技术,随机在2-、4-、8-细胞、桑椹胚和囊胚等5个阶检测其表达情况,结果表明: qRT-PCR结果与RNA-Seq结果基本一致（图5）。
显示原图|下载原图ZIP|生成PPT图5RNA-Seq差异表达基因的qRT-PCR验证
X-轴表示基因表达量,Y-轴表示胚胎发育阶段
-->Fig. 5Validation of sequencing results by qRT-PCR
X-axis denotes gene expression level, Y-axis denotes development stages
-->
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）条目上调。
显示原图|下载原图ZIP|生成PPT图6牦牛早期胚胎发育过程中的差异表达基因GO分类注释图
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：营养库
-->Fig. 6Gene Ontology classification of the DEGs during embryonic development
BP1: Cellular process; BP2: Single-organism process; BP3: Biological regulation; BP4: Metabolic process; BP5: Cellular component organization or biogenesis; BP6: Developmental process; BP7: Localization; BP8: Response to stimulus; BP9: Multicellular organismal process; BP10: Immune system process; BP11: Reproductive process; BP12: Locomotion; BP13: Biological adhesion; BP14: Multi-organism process; BP15: Behavior; BP16: Signaling; BP17: Growth; BP18: Rhythmic process; BP19: Hormone secretion; BP20: Biological phase; BP21: Cell aggregation; BP22: Cell killing; BP23: Reproduction; CC1: Cell part; CC2: Organelle; CC3: Organelle part; CC4: Membrane; CC5: Membrane part; CC6: Macromolecular complex; CC7: Extracellular region part; CC8: Cell junction; CC9: Extracellular region; CC10: Membrane-enclosed lumen; CC11: Synapse part; CC12: Synapse; CC13: Extracellular matrix; CC14: Cell; CC15: Collagen trimer; CC16: Extracellular matrix part; CC17:Nucleoid; CC18: Virion part; MF1: Binding; MF2: Catalytic; MF3: Molecular function regulator; MF4: Nucleic acid binding transcription factor; MF5: Molecular transducer; MF6: Transporter; MF7: Enzyme regulator; MF8: Protein binding transcription factor; MF9: Structural molecule; MF10: Guanyl-nucleotide exchange factor; MF11: Channel regulator; MF12: Electron carrier; MF13: Antioxidant; MF14: translation regulator; MF15: Receptor regulator; MF16: Chemoattractant; MF17: Chemorepellent; MF18: Protein tag; MF19: Metallochaperone; MF20: Morphogen; MF21: Nutrient reservoir
-->
牦牛胚胎发育过程中差异表达基因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个通路显著富集。
Table 4
表4
表4牦牛胚胎发育过程中差异表达基因富集前5（Top 5）KEGG通路表
Table 4Top 5 of enriched KEGG pathways of DEGs during embryonic development
| KEGG代谢通路 KEGG pathways | 上调 Up-regulated genes | 下调 Down-regulated genes | Q-值 Q-value |
|---|---|---|---|
| 2-细胞到4-细胞2-cell to 4-cell | |||
| 剪接体Spliceosome | 12 | 12 | 0.000 |
| RNA转运RNA transport | 18 | 23 | 0.000 |
| 泛素介导的蛋白水解Ubiquitin mediated proteolysis | 23 | 16 | 0.001 |
| 内吞作用Endocytosis | 40 | 15 | 0.006 |
| 氧化磷酸化作用Oxidative phosphorylation | 81 | 3 | 0.012 |
| 4-细胞到8-细胞4-cell to 8-cell | |||
| 嗅觉转导Olfactory transduction | 55 | 34 | 0.000 |
| 神经活性配体-受体互作Neuroactive ligand-receptor interaction | 60 | 29 | 0.000 |
| 核苷酸切除修复Nucleotide excision repair | 6 | 2 | 0.001 |
| 剪接体Spliceosome | 36 | 3 | 0.001 |
| 细胞因子及其受体互作Cytokine-cytokine receptor interaction | 35 | 52 | 0.009 |
| 8-细胞到桑椹胚8-cell to morula | |||
| 嗅觉转导Olfactory transduction | 17 | 64 | 0.000 |
| 泛素介导的蛋白水解Ubiquitin mediated proteolysis | 17 | 25 | 0.000 |
| 神经活性配体-受体互作Neuroactive ligand-receptor interaction | 19 | 73 | 0.000 |
| 核糖体Ribosome | 3 | 142 | 0.001 |
| RNA转运RNA transport | 11 | 44 | 0.001 |
| 桑椹胚到囊胚 Morula to blastocyst | |||
| 剪接体Spliceosome | 46 | 36 | 0.003 |
| RNA转运RNA transport | 52 | 45 | 0.020 |
| 嗅觉转导Olfactory transduction | 61 | 23 | 0.120 |
| 细胞周期-酵母Cell cycle - yeast | 19 | 28 | 0.122 |
| 蛋白转运Protein export | 6 | 5 | 0.122 |
新窗口打开
3 讨论
新一代高通量测序技术的不断发展,已彻底改变了转录组学的研究,使RNA-Seq无需预先设计探针即可对特定条件下任何生物生长发育阶段整体转录活动进行测序,准确探测到各种条件下的基因表达情况,发现了许多未知的分子调控机制[5,6,7]。本研究首次利用RNA-Seq技术从转录组学角度揭示牦牛早期胚胎发育机制,为完善牦牛胚胎体外生产提供新思路,同时也为进一步完善牦牛基因结构信息和胚胎发育相关的新基因提供理论基础。由于哺乳动物2-细胞到囊胚期每个胚胎的总RNA只有200―2 000 pg[8],这样微量RNA不能满足转录组测序文库的构建及高通量测序的基本要求,因此,分别提取2-细胞、4-细胞、8-细胞、桑椹胚和囊胚等5个发育阶段的牦牛胚胎总RNA,使用先进Smart-Seq2扩增技术对样本进行富集并构建测序文库[9],再应用RNA高通量测序技术对其进行高通量测序分析。测序质量评估、数据分析及qPCR验证结果都表明,测序质量和文库构建质量高,测序数据准确可靠。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主要涉及参与胞外配体信号（BMP4、BRINP3、FGF2、STAT5A、TBC1D24 和WBP1）、内吞作用和胞吐作用（MON1B、PMM2、SLC18A2和TTLL3）、凋亡调控（PARM1）、细胞应激保护（HSPA1L）、能量代谢（PCCB）、蛋白互作（C1QB和WBP1）和转录调控（IRF9）等途径影响胚胎发育能力[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]。本研究发现牦牛胚胎的NANOG在8-细胞期开始表达。CLDN4是Claudin蛋白家族的成员之一,在缝隙连接（gap junction）形成中起不可或缺的作用。研究证明,妇女黄体期子宫内膜CLDN4 mRNA的表达量与妊娠相关[33]。在小鼠胚胎培养液中添加CLDN4抑制因子会导致胚胎不能形成正常的囊胚腔[34]。在乳腺癌细胞活性研究中还发现CLDN4具有促进细胞的组织侵入作用[35]。因此,与在人的囊胚中检测结果类似[36],本研究发现在牦牛囊胚中CLDN4的表达量显著提高,说明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]。BMP15、KIT、GDF9、STAT3、ZP3和ZP4等基因被证明为标志性的母源基因[29-30,39],随着胚胎发育的进行,牦牛的这些母源基因表达量逐渐减少（图3）。在牦牛胚胎4-和8-细胞期开始表达的基因分别有1 221个和1 116个,而在桑椹胚和囊胚期开始表达的基因分别有142个和564个。另外,牦牛胚胎在2-细胞到4-细胞期以及4-细胞到8-细胞期即出现大量差异表达基因,但是EGA的标志性基因NANOG和CLDN4[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技术对牦牛体外受精胚胎发育不同阶段转录组进行了分析,获得了众多差异基因和有关通路的富集。不同阶段差异表达基因在数量、功能、分类和代谢通路都有各自的特异性。对于完善牦牛胚胎体外生产技术,及进一步完善牦牛基因结构信息和胚胎发育相关的新基因提供了理论基础。The authors have declared that no competing interests exist.
