谷胱甘肽促进牛体外受精胚胎发育的转录组初探

于学颖，郭芹芹，郝海生，孙尉峻，赵学明，朱化彬，杨　凌，杜卫华*

(1．河北工程大学，邯郸 056000； 2.中国农业科学院北京畜牧兽医研究所，北京 100193)



谷胱甘肽促进牛体外受精胚胎发育的转录组初探

于学颖1，2，郭芹芹2，郝海生2，孙尉峻2，赵学明2，朱化彬2，杨凌1*，杜卫华2*

(1．河北工程大学，邯郸 056000； 2.中国农业科学院北京畜牧兽医研究所，北京 100193)

摘要：旨在探究牛体外受精胚胎的培养液中添加3 mmol·L-1谷胱甘肽(Glutathione，GSH)后，8～16-细胞期和桑椹期的牛体外胚胎在转录组水平上的变化。收集8～16-细胞期和桑椹期的体外受精胚胎，构建文库，通过Illumina平台进行测序；筛选差异基因，并进行功能注释和代谢途径分析；采用定量PCR对10个基因的表达水平进行检测，以验证测序结果的准确性；在获得的新转录本中挑选3组进行RNA和蛋白结构的预测。通过OOSP1、THAP9等10个基因的定量检测、测序结果的质量评价(碱基错误率、Q20、Q30、GC含量)、reads与牛基因组的比对及其在基因组的分布统计分析，均说明该测序结果准确、可靠。测序获得差异表达基因4 100个(其中，处理组8～16-细胞期胚胎和桑椹胚中3 952个，对照组中884个)，已知基因2 225个，未知基因1 875个。这些差异基因主要富集到与ATP合成、呼吸链、DNA合成、翻译过程和物质代谢相关的84条GO terms上，参与细胞黏附、柠檬酸循环、谷胱甘肽代谢、溶酶体以及氨基酸代谢等27条通路。另外，与已知的转录本相比，5个新转录本中均发现不同长度的新外显子；预测的蛋白结构中除包含各自的保守结构域(丝氨酸/苏氨酸蛋白激酶的催化结构域、MAD同源结构域)外，还发现2个新的蛋白结构域(低复杂区、dwarfins蛋白家族B结构域)。综上表明，在培养液中添加GSH能导致胚胎转录谱的变化，显著提高胚胎的体外发育能力。

关键词：牛体外胚胎；转录组测序；定量PCR；新转录本；功能注释；代谢通路

胚胎发育是一个精密且高度协调的过程，为确保早期胚胎发育周期的有序进行，维持细胞内氧自由基(Reactive oxygen species，ROS)产生和清除的平衡至关重要。在正常的细胞代谢过程中，各个系统都会产生ROS并发挥一定的生物学作用，包括信号转导及功能调控等[1]。但ROS过量积累会对细胞产生有害影响，如DNA损伤、脂质中多不饱和脂肪酸氧化、氨基酸氧化、或氧化辅助因子对特定酶的灭活[2]。F.Calzi等[3]研究表明，体外培养时低氧环境下培养的胚胎发育能力要显著优于高氧环境；输卵管内的氧气浓度低于子宫内，约为大气压的40%[4-5]。

谷胱甘肽(Glutamine，GSH)是一种广泛存在于动物体细胞和配子中的三肽，主要有两种存在形式：氧化型(GSSG)和还原型(GSH)。还原型GSH会在谷胱甘肽过氧化物酶的作用下被氧化为GSSG，GSSG可通过谷胱甘肽还原酶再次被还原成GSH[6]。GSH能清除细胞内的ROS，并维持胞内的氧化还原平衡[7]；对于卵母细胞，GSH含量随其成熟的进程而逐渐升高[8]；对胚胎，胞内GSH和ROS含量能够决定合子的发育潜力[7]。因此，GSH常添加于卵母细胞成熟液和胚胎培养液中以提高胚胎的体外发育能力[9-10]，尤其是促进胚胎的致密化[11]。另外，GSH可通过表观遗传重编程调控胚胎的发育潜能[12]，降低胞内GSH含量，影响DNA的甲基化修饰[13]。关于GSH促进胚胎体外发育的研究仍停留在氧化还原平衡的层面，对胚胎内部发生的分子生物学事件、GSH对胚胎整个基因组在RNA水平的影响还没有得到完全地揭示。笔者前期的研究表明，体外培养液中添加GSH(3 mmol·L-1)能显著提高牛体外受精胚胎的囊胚率，但对GSH合成、代谢相关基因和胚胎发育相关基因的表达水平没有显著影响[10]。因此，对GSH处理胚胎的基因组转录水平进行整体性、系统性研究显得至关重要。随着高通量测序的不断发展和单细胞测序技术的诞生[14-15]，胚胎转录组研究成为可能。本研究以GSH处理的牛体外受精胚胎为对象，未处理组胚胎为对照，采用单细胞测序技术进行转录组测序，对测序结果进行验证；筛选差异表达基因，并进行基因功能注释、信号通路及网络分析；同时对新转录本的功能进行生物信息学预测。

1材料与方法

1.1主要试剂

本研究所用试剂，如无特殊说明均购自西格玛奥德里奇(上海)贸易有限公司。荧光定量PCR试剂盒Power SYBR®Green PCR Master Mix和Power SYBR®Green Cells-to-CtTMKit购自美国Life Technologies公司，细胞裂解试剂盒SMARTer®UltraTMLow RNA Kit for Illumina®Sequencing和cDNA全长扩增试剂盒Advantage®2 PCR Kit购自美国Clontech公司，纯化试剂盒AMPure XP beads购自美国Beckman Coulter公司。精液购自山东奥克斯生物技术有限公司。

1.2卵母细胞的采集、体外受精及胚胎培养

从当地屠宰场采集牛卵巢，置于保温瓶中，3 h内运回实验室，温度控制在35～37 ℃。使用37 ℃含有100 IU·mL-1青霉素和100 mg·mL-1硫酸链霉素的生理盐水将卵巢洗净，将18号针头与真空蠕动泵连接，抽取3～6 mm的窦状卵泡，体视显微镜下挑选卵丘卵母细胞复合体(Cumulus-oocyte complex，COC)，在成熟液(TCM199 medium、10% FBS、0.01 IU·L-1FSH、10 IU·L-1LH、1 μg·mL-1雌二醇、100 ng·mL-1IGF、50 ng·mL-1EGF、100 IU·mL-1青霉素、100 μg·mL-1链霉素)中洗3次，放入含750 μL成熟液的四孔板中，每孔放置50枚COC，在38.5 ℃、5% CO2、100%湿度的CO2培养箱中成熟培养22 h。

精液从液氮中取出后，空气中平衡10 s，投入37.5 ℃的水浴锅中进行解冻；用5 mL洗精液(m-BO medium[16]，10 mmol·L-1咖啡因、4 mg·mL-1BSA)进行稀释，然后500 g离心8 min，弃上清，再次离心洗涤；最后用受精液(m-BO medium，4 mg·mL-1BSA)将精子密度稀释到2×107个·mL-1。在含有15枚COC的50 μL受精滴中加入50 μL精液，于38.5 ℃、5% CO2、100%湿度的CO2培养箱中共孵育8 h。用吸卵针脱除外周的颗粒细胞后，受精卵移入前期培养液(NaCl 109.50 mol·L-1、KCl 3.10 mol·L-1、NaHCO326.20 mol·L-1、MgCl2·6H2O 0.80 mol·L-1、KH2PO31.19 mol·L-1、丙酮酸钠0.4 mol·L-1、葡萄糖1.5 mol·L-1、半乳糖酸钙5 mol·L-1、BSA 6 mg·mL-1、L-谷氨酰胺1 mol·L-1、2%必需氨基酸、1%非必需氨基酸)。48 h后，将4～8-细胞胚胎移入后期培养液(NaCl 109.50 mol·L-1、KCl 3.10 mol·L-1、NaHCO326.20 mol·L-1、MgCl2·6H2O 0.80 mol·L-1、KH2PO31.19 mol·L-1、丙酮酸钠0.4 mol·L-1、葡萄糖1.5 mol·L-1、半乳糖酸钙5 mol·L-1、10% FBS、L-谷氨酰胺1 mol·L-1、2%必需氨基酸、1%非必需氨基酸)继续培养，每隔48 h半量换液。培养过程中，分别于受精后48和132 h收集8～16-细胞胚胎和桑椹胚用于测序和荧光定量PCR检测。处理组胚胎的培养液中添加3 mmol·L-1GSH，对照组则不添加GSH。

1.3构建cDNA文库并测序

以3 mmol·L-1GSH处理的8～16-细胞期(C8C16_T)和桑椹期(ML_T)的胚胎为处理组，GSH未处理的同一阶段胚胎(C8C16_T和 ML_C)为对照组进行测序。每组20个胚胎，两个生物学重复，共8个cDNA文库进行测序。

胚胎细胞裂解后进行第一链cDNA合成，纯化后通过LD-PCR扩增得到双链cDNA。定量检测合格后，使用Covaris系统超声打断双链cDNA，获得的短片段进行末端修复、加A尾并连接测序接头，选择大小在200 bp左右的片段PCR富集得到最终的cDNA文库。建库后，使用Qubit2.0进行初步定量，再使用Agilent 2100检测文库的insert size，最后使用qPCR准确定量文库的有效浓度(>2 nmol·L-1)，保证文库质量。文库检验合格后，把不同文库按照有效浓度及目标上机数据量的需求混池后进行HiSeq/MiSeq测序。得到的原始序列中去除带接头的、低质量、带引物的序列，得到有效序列。选择牛参考基因组(版本号：UMD3.1)，通过TOPHAT(v2.0.6)软件将所有序列比对到参考基因组上。

1.4差异表达基因的GO功能注释和富集分析

Gene Ontology(GO)是基因功能国际标准分类体系，根据试验目的筛选差异基因后，研究差异基因在Gene Ontology中的分布状况。采用GOseq[17]分析方法，准确地计算出GO term被差异基因富集的概率，阐明样本差异在基因功能上的体现。

1.5差异表达基因的代谢通路富集分析

KEGG(Kyoto encyclopedia of genes and genomes)是有关Pathway的主要公共数据库[18]。通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。使用KOBAS(2.0)，以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。

1.6新转录本的功能预测

根据测序结果，选取BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1，NovelIso00290.2)和SMAD3(NovelIso00499.1，NovelIso00499.2)等与早期胚胎发育相关基因的新转录本进行简单的功能预测分析。

通过NCBI的ORF Finder(http：//www.ncbi.nlm.nih.gov/gorf/gorf.html)程序先对3个基因的新转录本进行开放阅读框(Open reading frame，ORF)分析，获得ORF的长度、起始和终止位置；再输入新转录本的序列信息，获得其编码的蛋白序列；采用SMART(http：//smart.embl-heidelberg.de)程序对蛋白序列进行功能结构域预测，与PFAM蛋白质结构域数据库进行多重序列比对。

1.7荧光定量PCR

根据测序结果中基因表达值进行初步筛选，结合 KEGG Pathway 和GO分析结果，确定与GSH代谢、氧化应激反应和细胞凋亡相关的10个基因进行荧光定量PCR反应，验证测序结果，其中GAPDH为内参基因。参考NCBI基因序列，利用primer 5.0设计引物，引物序列见列表1。

收集8～16-细胞(受精后48 h)期和桑椹胚(受精后132 h)期的牛胚胎，采用7900HT system系统 (Applied Biosystems，美国)进行检测。10 μL PCR反应体系：2×Master Mix 5 μL，上下游引物(10 μmol·L-1)各0.2 μL，cDNA模板1 μL，Rnase-Free Water 3.6 μL，混匀样品。PCR反应条件为95 ℃预变性30 s；95 ℃变性5 s，60 ℃退火15 s，72 ℃延伸15 s，40个循环；阴性对照用0.2 μL Rnase-Free Water代替模板。试验对每个样品检测进行3个技术重复。通过定量PCR获得Ct值，采用2-ΔΔCt计算得到基因的相对表达量。

表1　基因引物信息

1.8测序结果的验证

比较荧光定量PCR方法获得的8～16-细胞期与桑椹期胚胎的P-value和FC(Fold Change)，与转录组测序获得的相应结果进行比较。其中，P-value是通过比较两个时期的基因相对表达水平计算得到的，显著性分析按双侧P<0.05定义，FC=log2(桑椹胚时期基因的相对表达水平/8～16-细胞时期基因的相对表达水平)。全部数据采用SAS 9.2软件进行单因素方差分析。

2结果

2.1测序数据的质量评估

根据表2，以ML_T1_1为例，测序得到原始序列(Raw reads)42 426 942条，有效序列(Clean reads)36 950 922条，总的数据量(Clean bases)达3.65G。测序错误率(Error rate)为0.03%，低于0.5%，在可接受范围内。测序错误率小于1%的碱基(Q20)在整个数据中大于96%，错误率小于0.1%的碱基(Q30)占整个数据92%以上，即准确测序的碱基比例较高。GC含量39.36%，处于公认的40%～60%。总之，8个测序样品中错误率均控制在0.04%以下，Q20、Q30的比例均在88%以上，GC含量也在40%左右，因此该测序数据准确可靠，可用于后续的分析。

2.2有效序列与牛基因组比对结果

测序获得的有效序列与牛基因组AC_000160.1比对结果见表3。以处理组桑椹胚为例，其中约58 000 000条序列比对到参考基因组上，比对率约92%，其中88% 序列唯一比对到参考基因组上。在允许前6～7个碱基错配的前提下，8个样品中至少有86%的序列唯一比对到基因组上，其中部分序列与可变剪接有关，如处理组桑椹胚中约19%序列比对在剪接结合区。所以，8个样品中均有91%以上的序列比对到牛基因组上。

2.3有效序列在牛基因组上的分布情况

将比对到基因组上的有效序列按定位区域外显子(Exon)、 内含子(Intron)和基因间区(Intergenic)进行统计，见图1。以处理组桑椹胚为例，其中67.0%的序列分布在Exon上，10.9%的序列分布在Intron上，剩余22.1%分布在Intergenic。

图1　序列在基因组不同区域的分布Fig.1　Distribution of reads in different regions of bovine genome

表2　测序数据质量指标汇总

ML_T、ML_C、C8C16_T、C8C16_C分别为GSH处理组桑椹胚、对照组桑椹胚、GSH处理组8～16-细胞期胚胎和对照组8～16-细胞期胚胎。Raw reads为原始序列。Clean reads为有效数据。Error rate是测序错误率，Q20、Q30分别指测序错误率分别≤1%、≤0.1%的碱基数目比例ML_T、ML_C、C8C16_T、C8C16_C are groups of morulae with GSH treatment，morulae without treatments，8-16-cell stages embryos with or without GSH treated，respectively.Raw reads are obtained by conversion from raw sequence reads data.Clean reads refers to the valid data which discard the raw reads of low quality sequence reads and fitting contamination reads.Error rate refers to the sequencing error rate.Q20，Q30 are the proportion of bases with sequencing error rate ≤1% or ≤0.1% respectively

表3　测序数据与牛基因组的比对结果

ML_T、ML_C、C8C16_T、C8C16_C分别为GSH处理组桑椹胚、对照组桑椹胚、GSH处理组8～16-细胞期胚胎和对照组8～16-细胞期胚胎

ML_T、ML_C、C8C16_T、C8C16_C are groups of morulae with GSH treatment，morulae without treatments，8-16-cell stages embryos with or without GSH treated，respectively

2.4差异表达基因的筛选

来自4个组别、2次重复的8个样品的测序结果共获得314 135 947条序列，各个重复间的皮尔逊相关系数保证样品制备和测序技术的重复性。分别在8～16-细胞期胚胎和桑椹胚中检测到13 489和11 527个基因，检测到的基因约占牛基因组中基因(22 000个)的50%。

根据测序结果，比较处理组的8～16-细胞期胚胎和桑椹胚，共发现3 952个差异基因(2 156个已知基因和1 787个未知基因)；而比较对照组的8～16-细胞期胚胎和桑椹胚，则发现884个差异基因(363个已知基因和521个未知基因)，其中736个基因在两个组中均存在(图2)。在 4 100个差异基因中，已知基因2 225个，未知基因1 875个。已知差异基因中，60个仅在对照组中表达，1 862个仅在处理组表达，303个在两组中均表达。处理组中，有下调基因193个，上调基因170个；对照组中，1 148个基因下调，1 017个基因上调。

2.5差异表达基因的功能注释

通过GO分析发现，获得的差异基因主要富集到与ATP合成、呼吸链、DNA合成、翻译过程和物质代谢相关的84条GO terms上，如质子转运ATP合成酶复合物(Proton-transporting ATP synthase complex)、NADH脱氢酶复合物(NADH dehydrogenase complex)、呼吸链复合体I(Respiratory chain complex I)、DNA聚合酶活化(DNA polymerase activity)、线粒体功能(Mitochondrial part)、核糖体代谢过程(Ribosome metabolic process)等，而对照组中的差异基因只富集到与呼吸链和ATP相关的16条terms上，如呼吸链(Respiratory chain)，ATP生物合成过程(ATP biosynthetic process)，ATP合成与质子转运(ATP synthesis coupled proton transport)等。

左边区域为处理组中特异存在的差异基因，右边区域为对照组中特异存在的差异基因，重叠区域为在两个组中均存在的差异基因The left area represents the differentially expression genes in treated group，the right area represents the differentially expressed genes in control group，and the overlapping area represents common genes of differentially expressed in two groups图2　差异基因维恩图Fig.2　Venn diagram of differentially expression genes

2.6差异表达基因的代谢通路和信号转导途径分析

通过KEGG Pathway信号通路分析，共发现27条通路，对照组中有矿物质吸收、细胞黏附、叶酸合成等途径，而处理组则发现了柠檬酸循环(TCA循环)、谷胱甘肽代谢、溶酶体以及氨基酸代谢等多条通路(表4)。

表4　差异基因KEGG富集结果

2.7新转录本的生物信息学分析

转录组测序共获得325 585条新转录本，根据新转录本信息和文献查阅，确定3组与早期胚胎发育相关基因的新转录本进行功能预测分析，包括BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1，NovelIso00290.2)和SMAD3(NovelIso00499.1，NovelIso00499.2)。

首先对3个基因的新转录本进行ORF预测，BCL6(NovelIso00109.1)ORF长246 bp(图3a)，起始密码子位于 1 bp处，终止密码子位于246 bp处，编码81个氨基酸残基。GSK3B(NovelIso00290.1)ORF长735 bp(图3b)，起始密码子位于2 bp处，终止密码子位于736 bp 处，编码244个氨基酸残基。GSK3B(NovelIso00290.2)ORF的起始密码子位于2 394 bp 处，终止密码子位于2 870 bp 处，长477 bp 的 ORF(图 3c)，编码158个氨基酸残基。SMAD3(NovelIso00499.1)ORF开始于1 806 bp 处，终止于2 768 bp 处，长963 bp，编码320个氨基酸残基(图 3d)。SMAD3(NovelIso00499.2)也找出一条同样长度的ORF(图 3e)，但起始位点不同，起始于662 bp处，终止于1 624 bp处，编码320个氨基酸残基。

a～e.BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1)、GSK3B(NovelIso00290.2)、SMAD3(NovelIso00499.1)和SMAD3(NovelIso00499.2)序列的ORF分析结果a-e.The ORF analysis of BCL6 (NovelIso00109.1)，GSK3B (NovelIso00290.1)，GSK3B (NovelIso00290.2)，SMAD3 (NovelIso00499.1)，SMAD3(NovelIso00499.2) sequences,respectively图3　基因序列的 ORF 分析Fig.3　ORF analysis of genes sequences

采用SMART 程序对基因编码序列进行功能结构域的预测，发现在GSK3B(NovelIso00290.1)基因的氨基酸序列中29～243位发现丝氨酸/苏氨酸蛋白激酶的催化结构域(Serine/Threonine protein kinases catalytic domain，S-TKc)(图4a)。在GSK3B(NovelIso00290.2)基因氨基酸序列29～46位发现低复杂区(Low-complexity region，LCR)(图4b)。在SMAD3(NovelIso00499.1)基因的氨基酸序列中分别发现了1～26位的MAD 同源(MAD homology 1，MH1)结构域、28～36位LCR结构域 和 125～296位的dwarfins蛋白家族的B结构域(Domain B in dwarfin family proteins，DWB)(图4c)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同样位点发现同样的结构域。BCL6(NovelIso00109.1)编码的蛋白序列未发现功能结构域。

a.GSK3B(NovelIso00290.1)基因中预测到的丝氨酸/苏氨酸蛋白激酶的催化结构域(S-TKc)；b.GSK3B(NovelIso00290.2)基因中预测到的低复杂区(LCR，粉色框)；c.SMAD3(NovelIso00499.1)基因中发现的MAD同源(MH1)结构域、LCR结构域和dwarfins蛋白家族B结构域(DWB)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同样的结构域。BCL6(NovelIso00109.1)基因编码的蛋白序列中未发现功能结构域a.The S-TKc in GSK3B (NovelIso00290.1) predicted；b.The LCR in GSK3B (NovelIso00290.2)；c.The MH1，LCR and DWB in SMAD3 (NovelIso00499.1).The same domains were found in SMAD3 (NovelIso00499.2).However there was no domains found in BCL6 (NovelIso00109.1) sequence图4　基因编码序列功能结构域预测Fig.4　Prediction of the function domains of genes

表5　定量PCR检测结果与转录组比较结果

2.8基因的定量PCR与测序结果比对

通过比较整体的变化趋势和基因表达差异倍数来判断定量PCR与转录组测序结果的一致性。当基因的 FC值在定量PCR和测序结果中都>1或<1时，说明该基因在两种检测中均为上调或下调，表达量增加或降低，这两种情况均表明定量PCR检测结果与转录组测序是一致的。经过比对，所选10个基因的表达倍数及整体变化趋势与测序结果均保持一致，表明转录组测序结果的可靠性。

3讨论

在胚胎培养液中添加GSH能提高牛体外受精胚胎的囊胚率和发育能力[10，19]，然而其分子机制并没有完全揭示。本研究采用单细胞测序技术，揭示GSH添加后牛胚胎整个基因组的基因表达模式变化，获得大量的差异表达基因，其中14个差异基因富集到谷胱甘肽的代谢通路上，此外差异基因还富集到26条其他相关信号传导和代谢通路上，如柠檬酸循环(TCA循环)、溶酶体以及氨基酸代谢等。可见在培养液中添加GSH后，除自身的代谢途径外还引发了其他的对8～16-细胞期胚胎发育至桑椹胚有重要作用的生物学过程，促进胚胎的体外发育，显著提高体外胚胎的囊胚率。

转录组测序因其具有全方位、高通量等优点而广泛应用于动植物研究中，尤其在基因功能及其调控网络的揭示方面有着重要地位，然而单细胞测序在动物胚胎发育调控方面的研究还在起步阶段。2011年果蝇单囊胚的转录组测序揭示基因mRNA的表达丰度与性别决定机制的关系[20]。人类早期胚胎的转录组测序阐明了外胚层细胞和胚胎干细胞分子调控机制的差异，证实胚胎干细胞的多能性，为人类早期胚胎发育和干细胞研究提供借鉴[21]。对于家畜牛单胚胎的转录组测序首次发现等位基因的序列存在差异，同时获得大量的新SNP 位点[22]；而且，附植前的胚胎基因组在转录水平上发生巨大变化[23]。综上表明，转录组测序技术在基因表达、基因变异和早期胚胎发育的转录调控中有着广泛的应用前景。

本研究中所有测序结果的各项质量指标均在正常范围内。测序样品中碱基错误率均控制在0.04%以下，Q20、Q30的比例均在88%以上，GC含量也在40%左右。测序结果与牛基因组比对发现，8个样品中均有91%以上的序列比对到牛基因组上，其中唯一比对的序列占总数的86%，说明测序数据覆盖整个牛基因组。另外，至少有67%以上的序列分布到外显子区，少数的则定位到内含子和基因间区，可能是出现的新转录本或非编码RNA，或者是目前已知基因的内含子中尚未发现的新外显子。然而，转录组测序数据量巨大，每一个生物信息学分析步骤都可能出现假阳性结果，所以有必要采用定量PCR对其进行验证，这也是评定测序结果的重要标志之一。本研究对10个基因的RNA表达水平进行定量检测，得到的结果均与测序结果相吻合，充分说明测序的可靠性和准确性。

BCL6与基因表达调控和B细胞受体信号通路有关。测序中发现BCL6新转录本缺失了第4、5、6、7、9、10外显子序列，在第8个外显子之后发现长度为360 bp的新外显子，说明GSH可能影响BCL6基因表达情况。

GSK3B可调节牛卵母细胞减数分裂进程，参与卵母细胞和卵丘细胞中MAPK信号通路[24]。本研究发现GSK3B(NovelIso00290.1)转录本比已知序列多出一个外显子(长114 bp)，GSK3B(NovelIso00290.2)转录本中也发现两个新外显子(长度分别为2 187和689 bp)。在蛋白水平，两个新转录本都存在S-TKc元件，它是GSK3蛋白家族中存在的一个保守结构域，可作为分子结合位点(如ATP的结合位点)，或氨基酸多肽的结合位点参与蛋白磷酸化过程，与多数细胞活性有着密切的联系。LCR也是在新转录本中发现的由简单氨基酸组成的蛋白质序列[25]，虽然其功能模型尚未清楚，但有研究表明它在生物学功能中发挥关键作用[26-28]。A.Coletta研究表明，含有LCR结构的蛋白在生物总库的交互数据集(The biological general repository for interaction datasets，the BioGrid)中分布更加广泛[25]。

SMAD3是SMAD蛋白家族成员，该家族蛋白根据功能分为受体调节型、通用型和抑制型3类，SMAD3则是受体调节型。测序结果中SMAD3(NovelIso00499.1)新转录本中发现长度为1 695和193 bp的新外显子，而SMAD3(NovelIso00499.2)转录本中发现长度为551和193 bp的新外显子。新转录本的蛋白结构预测后发现MH1和DWB结构域，前者是SMAD3蛋白中存在的一个保守结构域，主要参与TGF-β信号通路的传导[29-30]；后者是dwarfins蛋白家族共有结构域，与哺乳动物中TGF-β磷酸化和细胞生长的调控有关[31]。在核因子I(NF-1)或CCAAT盒结合的转录因子(CTF)中也发现DWB结构域的存在。

综上表明，本研究应用RNA-seq技术对GSH处理的8～16-细胞和桑椹期的体外胚胎进行测序，在验证数据准确性和代表性的基础上，筛选与GSH处理相关的差异表达基因，对其生物学功能和代谢途径进行富集分析；发现大量的新转录本，也对其RNA和蛋白结构进行生物信息学预测；为从分子水平揭示GSH提高胚胎抗氧化能力和促进胚胎发育的机制奠定研究基础。

参考文献(References)：

[1]DEVASAGAYAM T P，TILAK J C，BOLOOR K K，et al.Free radicals and antioxidants in human health：current status and future prospects[J].JAssocPhysiciansIndia，2004，52：794-804.

[2]BROOKER R J.Genetics：analysis and principles[M].New York：McGraw-Hill Science，2011.

[3]CALZI F，PAPALEO E，RABELLOTTI E，et al.Exposure of embryos to oxygen at low concentration in a cleavage stage transfer program：reproductive outcomes in a time-series analysis[J].ClinLab，2012，58(9-10)：997-1003.

[4]FISCHER B，BAVISTER B D.Oxygen tension in the oviduct and uterus of rhesus monkeys，hamsters and rabbits[J].JReprodFertil，1993，99(2)：673-679.

[5]LEESE H J.Metabolic control during preimplantation mammalian development[J].HumReprodUpdate，1995，1(1)：63-72.

[6]LUBERDA Z.The role of glutathione in mammalian gametes[J].ReprodBiol，2005，5(1)：5-17.

[7]TAKAHASHI M，NAGAI T，OKAMURA N，et al.Promoting effect of beta-mercaptoethanol oninvitrodevelopment under oxidative stress and cystine uptake of bovine embryos[J].BiolReprod，2002，66(3)：562-567.

[8]DE MATOS D G，FURNUS C C.The importance of having high glutathione (GSH) level after bovineinvitromaturation on embryo development effect of beta-mercaptoethanol，cysteine and cystine[J].Theriogenology，2000，53(3)：761-771.

[9]GARDINER C S，SALMEN J J，BRANDT C J，et al.Glutathione is present in reproductive tract secretions and improves development of mouse embryos after chemically induced glutathione depletion[J].BiolReprod，1998，59(2)：431-436.

[10]SUN W J，PANG Y W，LIU Y，et al.Exogenous glutathione supplementation in culture medium improves the bovine embryo development afterinvitrofertilization[J].Theriogenology，2015，84(5)：716-723.

[11]PERREAULT S D，BARBEE R R，SLOTT V L.Importance of glutathione in the acquisition and maintenance of sperm nuclear decondensing activity in maturing hamster oocytes[J].DevBiol，1988，125(1)：181-186.

[12]TIMME-LARAGY A R，GOLDSTONE J V，IMHOFF B R，et al.Glutathione redox dynamics and expression of glutathione-related genes in the developing embryo[J].FreeRadicBiolMed，2013，65：89-101.

[13]LERTRATANANGKOON K，WU C J，SAVARAJ N，et al.Alterations of DNA methylation by glutathione depletion[J].CancerLett，1997，120(2)：149-156.

[14]PIRAS V，TOMITA M，SELVARAJOO K.Transcriptome-wide variability in single embryonic development cells[J].SciRep，2014，4：7137.

[15]DE DOMENICO E，OWENS N D，GRANT I M，et al.Molecular asymmetry in the 8-cell stage Xenopus tropicalis embryo described by single blastomere transcript sequencing[J].DevBiol，2015，408(2)：252-268.

[16]BRACKETT B G，OLIPHANT G.Capacitation of rabbit spermatozoainvitro[J].BiolReprod，1975，12(2)：260-274.

[17]YOUNG M D，WAKEFIELD M J，SMYTH G K，et al.Gene ontology analysis for RNA-seq：accounting for selection bias[J].GenomeBiol，2010，11(2)：R14.

[18]KANEHISA M，ARAKI M，GOTO S，et al.KEGG for linking genomes to life and the environment[J].NucleicAcidsRes，2008，36：D480-D484.

[19]ALI K，ABD-ALLAH S，KANDEIL M，et al.Impact of exogenous glutathione duringinvitromaturation，fertilization and development of buffalo oocytes[M].Fifth Scientific Conference，Ras Elbar，2007：1-12.

[20]LOTT S E，VILLALTA J E，SCHROTH G P，et al.Noncanonical compensation of zygotic X transcription in early Drosophila melanogaster development revealed through single-embryo RNA-seq[J].PLoSBiol，2011，9(2)：e1000590.

[21]YAN L，YANG M，GUO H，et al.Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells[J].NatStructMolBiol，2013，20(9)：1131-1139.

[22]CHITWOOD J L，RINCON G，KAISER G G，et al.RNA-seq analysis of single bovine blastocysts[J].BMCGenomics，2013，14：350.

[23]JIANG Z，SUN J，DONG H，et al.Transcriptional profiles of bovineinvivopre-implantation development[J].BMCGenomics，2014，15：756.

[24]UZBEKOVA S，SALHAB M，PERREAU C，et al.Glycogen synthase kinase 3B in bovine oocytes and granulosa cells：possible involvement in meiosis duringinvitromaturation[J].Reproduction，2009，138(2)：235-246.

[25]COLETTA A，PINNEY J W，SOLIS D Y，et al.Low-complexity regions within protein sequences have position-dependent roles[J].BMCSystBiol，2010，4：43.[26]FONDON J W III，GARNER H R.Molecular origins of rapid and continuous morphological evolution[J].ProcNatlAcadSciUSA，2004，101(52)：18058-18063.

[27]VERSTREPEN K J，JANSEN A，LEWITTER F，et al.Intragenic tandem repeats generate functional variability[J].NatGenet，2005，37(9)：986-990.

[28]PHATNANI H P，GREENLEAF A L.Phosphorylation and functions of the RNA polymerase II CTD[J].GenesDev，2006，20(21)：2922-2936.

[29]ATTISANO L，LEE-HOEFLICH S T.The Smads[J].GenomeBiol，2001，2(8)：3010.

[30]BROWN K A，PIETENPOL J A，MOSES H L.A tale of two proteins：differential roles and regulation of Smad2 and Smad3 in TGF-beta signaling[J].JCellBiochem，2007，101(1)：9-33.

[31]YINGLING J M，DAS P，SAVAGE C，et al.Mammalian dwarfins are phosphorylated in response to transforming growth factor beta and are implicated in control of cell growth[J].ProcNatlAcadSciUSA，1996，93(17)：8940-8944.

(编辑程金华)

doi:10.11843/j.issn.0366-6964.2016.07.008

收稿日期：2016-01-22

基金项目：国家科技支撑项目(2012BAD12B01-2)；基本科研业务费重点项目(2013ywf-zd-2)；家畜胚胎工程与繁殖创新团队(ASTIP-IAS06-2016)

作者简介：于学颖(1989-)，女,河北廊坊人，硕士，主要从事家畜胚胎工程研究，E-mail: xueyingyu2015bj@163.com *通信作者：杜卫华,博士,副研究员,E-mail：dwh@iascaas.net.cn；杨凌,博士,副教授,E-mail:yangling@hebeu.edu.cn

中图分类号：S823

文献标志码：A

文章编号:0366-6964(2016)07-1363-10

Transcriptome of Bovine IVF Embryos Treated with Glutathione

YU Xue-ying1，2，GUO Qin-qin2，HAO Hai-sheng2，SUN Wei-jun2，ZHAO Xue-ming2，ZHU Hua-bin2，YANG Ling1*，DU Wei-hua2*

(1.CollegeofAgriculture，HebeiUniversityofEngineering，Handan056000，China；2.InstituteofAnimalScience，ChineseAcademyofAgriculturalSciences，Beijing100193，China)

Abstract:To elucidate the effects of exogenous glutathione (GSH) on embryos transcriptome，bovine IVF embryos at 8-16-cell and morula stage were collected for high-throughput RNA sequencing (RNA-seq).Quantitative determination of genes expression was used to verify the RNA-seq data.Gene ontology and KEGG pathway enrichment analysis of differential expression genes were performed.The accuracy and validity of RNA-seq data were confirmed with combination of error rate of base，mapping of reads to bovine genome，distribution region analysis and consistency between qPCR and RNA-seq results of ten genes.A total of 4 100 genes were differentially expressed between embryos at two stages and 3 952 genes were identified as expressed in embryos treated with GSH.A high enrichment of genes involved in 84 GO terms，such as ATP synthesis，respiratory chain，DNA synthesis and translation.Pathway analysis revealed 27 pathways involved in cell adhesion，TCA cycle，GSH metabolism and amino acid metabolism.Additionally，2 known conserved domains (S-TKc，MH1) and 2 new domains (LCR，DWB) were found in 5 new transcripts by prediction，compared with known transcripts.To sum up，GSH treatment can result into the comprehensive change at transcriptional level of bovine IVF embryos and improve their development consequently.

Key words:bovine IVF embryo；RNA-Seq；qRT-PCR；novel isoform；GO analysis；KEGG pathway analysis