奶牛胚胎早期死亡的转录组初探

吕小青，刘林，薛建华，王彦平，李艳华，麻柱

（1.北京奶牛中心，北京 100192；2.农业农村部奶牛遗传育种重点实验室，北京 100192;3.奶牛遗传育种与繁殖北京市重点实验室，北京 100192）

早期胚胎丢失是奶牛繁殖失败的主要原因之一，主要发生在胚胎附植前后，有多个因素影响胚胎丢失，如遗传（SNPs、CNV等）、环境、疾病、病原微生物等都能引起胚胎死亡[1,2]。受精失败引起的胚胎损失达17%，卵裂阶段的胚胎损失达37%[3,4]。本研究筛选了两个不同卵裂阶段，对奶牛早期死亡胚胎进行了转录组分析。

转录组测序现在已经被广泛应用于人类、小鼠、猪、牛等动物研究中[5～7]，能快速地获得样本的序列信息和表达信息，具有高通量、全面等优点，为研究基因表达调控的重要方法。本研究的目的是挖掘引起奶牛胚胎早期死亡的关键基因或突变，目前国内还没有相关的报道，具体包括采用转录组测序技术对奶牛8-细胞期和桑葚胚期死亡和正常的胚胎样本进行基因表达变化的分析，本研究分别针对奶牛胚胎的两种不同时期进行细胞组学测序，联合分析不同时期胚胎之间的差异表达基因，获得与奶牛胚胎发育相关的生物标志物，为奶牛群体筛选奠定基础。

1 材料与方法

1.1 体外胚胎的生产

屠宰场收集的卵巢通过保温瓶运回实验室后，用10mL注射器采集卵母细胞。将鉴定可用的卵母细胞洗涤4遍后，放入成熟液滴中，成熟培养22～24h，并利用鲜精体外受精后，精卵共培养6～8h。之后将受精后的合子转移到胚胎发育液滴中进行培养，培养条件为38.5℃、5% CO2、饱和湿度，期间每隔48h半量换液。

在受精88～96h后收集8-细胞阶段正常胚胎和死亡胚胎，并在受精第6天时收集正常桑葚胚和死亡桑葚胚，提取RNA。胚胎发育阶段和胚胎分级评定标准参考国际胚胎移植协会标准。

1.2 cDNA文库构建及测序

利用链霉蛋白酶将胚胎的透明带消化掉，并利用玻针将单个卵裂球转移至裂解液中，立即冻存于-80℃冰箱。之后对单个卵裂球进行裂解，并反转录成1st cDNA。通过PCR扩增富集、核酸纯化、构建文库等之后，利用Illumina平台进行测序[8]，测序策略为PE 150。

1.3 单细胞转录组测序及数据处理

测序完成后，用 Casava1.8.2 软件将图像格式的测序数据转化为fastq格式，得到的原始数据进行过滤获得高质量Clean Reads，再进行后续分析，后续分析都基于Clean Reads[9,10]。本研究对获得的测序数据与牛基因组参考序列（UMD3.1）进行比对分析，采取HiSAT2软件将各样品过滤后的RNA-seq数据同基因组进行比对分析，比对率达92%。每个组设置两个生物学重复，采用DESeq2进行基因差异表达分析，并选取|log2Ratio|≥1和q＜0.05的基因作为显著差异表达基因，从而获得上下调基因个数。

2 结果与分析

2.1 8-细胞期胚胎和桑葚胚表达谱变化

通过对奶牛8-细胞期和桑葚胚时期死亡和正常的胚胎样本进行基因转录组测序，联合分析不同时期胚胎之间的差异表达基因，挖掘奶牛胚胎早期死亡的关键基因。

本研究设置了8-细胞期死亡组和正常组、桑葚胚死亡组和正常组，共4个组别，每组两个生物学重复，因此共建立了8个cDNA文库，并且利用了Illumina平台进行测序。8个测序样本共计获得52、257、166条Reads，之后用DESeq2进行差异基因分析，并选取|log2Ratio|≥1和q＜0.05的为差异显著。

通过8-细胞期正常组与死亡组，桑葚胚期正常组与死亡组差异基因比较分析，获得上下调基因个数。分析结果显示，在8-细胞期获得2倍以上差异表达基因157个，其中上调基因6个，下调基因151个；在桑葚胚期获得2倍以上差异表达基因269个，其中上调基因10个，下调基因259个。

2.2 差异表达基因的GO统计

基因功能注释（Gene Ontology，GO）可分为三个部分，分别描述基因的分子功能（Molecular Function，MF）、细胞组分（Cellular Component，CC）、生物过程（Biological Process，BP）。本研究通过差异基因分析在CC分类中有14个条目表达差异，其中差异基因占比最大的是细胞部分，其次为细胞器；在BP分类中有18个条目表达差异，其中差异基因占比最大的是繁殖过程；在MF分类中有10个条目表达差异，其中差异基因占比最大的是结合，其次是催化活性。GO详细统计结果如图1。

图1 差异基因的GO统计柱状图

2.3 差异基因KEGG通路分析

本研究对8-细胞期正常胚胎和死亡胚胎的差异表达基因，进行了KEGG通路分析，同时对桑葚胚期也进行了KEGG信号分析，与胚胎早期死亡相关性较高的通路为细胞凋亡信号通路和内质网蛋白合成通路。

3 讨论

Shinji Sasaki等对日本和牛人工授精30～60d期间死亡的791个胚胎进行了以CNV为变异类型的GWAS分析，发现定位到8号染色体上378 127到412 061bp位置的ANXA 10基因第2～6外显子区域的拷贝数变化（CNV） 与胚胎死亡有关。这项研究确定了一种缺失型CNV，在奶牛中包含ANXA10，与人工授精后30～60d的胚胎死亡率相关[11]。

2011年，VanRaden等人利用北美58 453头荷斯坦牛、5 288头娟姗牛及1 991头瑞士褐牛的50K Illumina SNPs基因型数据[12]，通过软件（Findhap.f90 version2）将其转化为单倍型数据，得到长度小于75个标记的单倍型，通过不同世代个体比对，发现11种群体频率很高的单倍型从来没有以纯合的形式出现，然后进行系谱分析，观察其遗传方式，并通过和具体繁殖性状（SCR、妊娠率、死胎率）的关联分析，推测出5种单倍型，当其纯合时，会导致早期流产或死胎。

本研究利用高通量测序对8-细胞期、桑葚胚期的正常和死亡胚胎进行了单细胞转录组测序，通过对差异表达基因的通路分析，发现细胞凋亡信号通路差异极显著。但是转录组测序数据量大，下一步将利用定量PCR技术对获得的差异显著基因进行验证，以进一步获得引起胚胎早期死亡的重要标志基因。

本研究初步对早期死亡胚胎进行了转录组分析，为进一步研究胚胎发育潜能和体外胚胎生产提供了基础，对提高胚胎生产效率具有重要意义。