基于RNA-Seq 技术的犏牛囊胚冷冻前后单核苷酸多态性和可变剪切分析

郑 杰, 张国忠, 胡艳梅, 伍学一, 杨宗富, 字向东

(1. 攀枝花市妇幼保健院生殖医学中心, 攀枝花 617000;2. 西南民族大学动物科学国家民委重点实验室, 成都 610041)

随着人类辅助生殖技术的飞速发展，胚胎超低温保存已成为体外受精(in vitrofertilization，IVF)临床治疗中广泛应用的方法[1]。同时，胚胎超低温保存技术在女性生育力保存、供卵、降低多胎妊娠和卵巢刺激的风险等方面具有重要意义，是目前唯一被美国生殖医学会认可的临床开展生育力保存的方法[2-3]。玻璃化冷冻法(vitrification)因具有冷冻效率高、操作简便、能有效避免冰晶的形成等优点，成为目前应用最广泛、冷冻效率最高的胚胎超低温保存方法[4]。然而即便选用玻璃化冷冻法对胚胎进行超低温冷冻，胚胎也会不可避免地遭受一定程度的冷冻损伤[1]。Zheng 等[5]研究了571个冷冻囊胚移植周期后发现，经两次冷冻后的囊胚相比新鲜囊胚，移植后的流产率较高，二者存在显著差异。另外还有研究表明，玻璃化冷冻会在一定程度上影响生物细胞增殖、分化和凋亡等相关基因表达，从而间接影响胚胎的发育能力[4-6]。目前，临床上一般都用胚胎发育情况、临床妊娠率和活产率等指标评估冻融胚胎质量，但不能从分子水平上对胚胎进行评估，不能有效地阐述胚胎的冷冻损伤机制。

单核苷酸多态性(single nucleotide polymorphism,SNP)指在DNA 或RNA 水平上，等位基因发生突变而引起DNA序列多态性的现象，包括颠换、置换、缺失和插入。SNP 是转录组中十分常见的变异类型，其特点是数量多、分布广、可稳定遗传[7]。可变剪切(alternative splicing, AS)则是指采用不同剪切方式使得mR N A 前体产生多种mRNA剪切异构体的过程。AS普遍存在于真核生物中，具有调节基因表达和蛋白质多样性的重要作用，人类AS 的基因已达到总基因的95%[8]。对SNP和AS的研究能使人们更便捷、准确地从分子水平来探究冻融囊胚的损伤机制。近年来，RNA转录组测序(RNA transcriptome sequencing，RNASeq)技术快速发展，除了可以准确分析基因的表达水平外，还可用于发掘相关的功能基因，预测新转录本，以及分析基因的SNP 和AS 等[9]。

哺乳动物如牛、羊、猪等的胚胎形态、大小、结构和功能与人类胚胎相似，哺乳动物胚胎的研究结果能对人类胚胎研究提供一定的理论基础和技术支持。目前，研究者已成功利用RNA-Seq技术在牦牛和犏牛胚胎发育调控[10-11]，以及牦牛和犏牛胚胎冷冻损伤[12-13]方面做了初步研究，结果提示新鲜胚胎经玻璃化冷冻后SNP和AS发生了相应变化，但这些研究并未对冷冻前后的SNP和AS 进行比较分析。因此，本研究选用牦牛卵母细胞和娟姗牛精子进行胚胎体外培养，再对杂交二代犏牛囊胚进行玻璃化冷冻实验，利用RNA-Seq技术对玻璃化冷冻前后犏牛囊胚SNP的数目和分布特征，以及AS 的数目、事件类型和重要差异表达基因(differentially expressed genes，DEGs)进行比较分析，以期为后续人类胚胎冷冻前后相关基因功能分析和挖掘提供有效数据和理论基础，为完善人类胚胎的玻璃化冷冻技术提供新思路。

1 材料与方法

1.1 材料

1.1.1 卵巢与冻精来源 牦牛卵巢来源于四川省广汉市向阳屠宰场，娟姗牛冷冻精子购自四川省家畜改良站。

1.1.2 试剂 促卵泡生成素(follicle-stimulating hormone, FSH)(Folltropin®-V)和促黄体素(luteinizing hormone, LH)(Lutropin®-V)均购自加拿大Bioniche公司；青霉素-链霉素双抗和胎牛血清均购自美国Gibco 公司；雌二醇(β-estradiol)、动物细胞培养液Medium 199(10×)和丙酮酸钠均购自美国Sigma 公司；透明质酸酶(HYASE，10×)、洗精液(SpermRinseTM)、洗精受精液(G-IVFTMPLUS)、卵裂胚培养液(G-1TMPLUS)和囊胚培养液(G-2TMPLUS)均购自瑞士Vitrolife 公司；玻璃化冷冻试剂盒(V T 101)和玻璃化解冻试剂盒(VT102)为日本Kitazato 公司产品。

1.2 方法

1.2.1 犏牛囊胚的体外生产 将屠宰场采集的牦牛卵巢放于35 ℃盛有DPBS缓冲液的恒温保温壶内，确保壶内温度恒定，2 h 内送回实验室，继续后续实验。挑选直径为3～8 mm 的卵泡，抽取其卵泡液，挑选外周至少被3 层颗粒细胞包裹且细胞质均匀的卵丘卵母细胞，置于38.5 ℃、5% CO2及饱和湿度的培养箱中培养24 h，使牦牛卵母细胞体外成熟(in vitromaturation, IVM)[14]。

选取质量好的成熟卵母细胞，以20 枚/滴的量转入G-IVFTM微滴中，然后加入已获能的娟姗牛精子，置于38.5 ℃、含5% CO2和90% N2且饱和湿度的三气培养箱内培养22～24 h 以完成IVF。用0.2%透明质酸酶去除受精卵上的颗粒细胞，转移至G-1TM微滴并置于三气培养箱中培养72 h，最后转移至G-2TM微滴并置于三气培养箱中培养48 h，得到犏牛囊胚。

1.2.2 样品的收集 选取经IVF 发育良好且未孵化的3 枚犏牛新鲜囊胚作为新鲜组，使用玻璃化冷冻试剂盒对其余犏牛新鲜囊胚进行玻璃化冷冻操作。将囊胚在液氮中储存3 d 后，再使用玻璃化解冻试剂盒解冻囊胚，选取解冻后完全复苏、未孵化的3 枚犏牛冻融囊胚作为冻融组。用PBS缓冲液清洗两组囊胚，再分别移至内含细胞裂解液的两个单细胞采集管中，超低温保存。

1.2.3 文库的构建及测序 在相同实验条件下，采用TRIzol 法[15]分别提取新鲜组和冻融组犏牛囊胚的总RNA，然后由安诺优达基因科技(北京)有限公司采用Smart-Seq2 方法[16]进行反转录扩增。检测扩增产物的分布情况，判定其质量，以20 ng cDNA 为模板进行后续文库的构建。用Bioruptor超声系统(比利时Diagenode公司产品)将cDNA打断成200 bp 左右的小片段，然后进行cDNA 末端修复、加碱基A、测序接头拼接等后续步骤。选取接头产物，进行PCR 扩增。扩增产物经2%琼脂糖凝胶电泳检测后，使用CWBIO 凝胶回收试剂盒回收DNA，完成最终文库的构建。最后在HiSeqTM2500 测序平台上运行双端测序程序(PE100)，对犏牛新鲜囊胚和冻融囊胚两个文库进行高通量测序。

1.2.4 转录组数据分析 经过转录组文库质检、成簇与测序，运用CASAVA软件[17]对Illumina测序所得的原始数据图像进行分析，得到原始序列数据，即Raw Reads。为保证信息分析数据的质量，需要对原始数据进行一系列数据过滤步骤，最终得到Clean Reads，再进行后续分析。采用测序饱和度分析和测序随机性分析，对测序数据进行评估。运用TopHat 软件[18]将测序过滤后得到的Clean Reads 与牦牛参考基因组(BosGru-v2.0)[19]进行比对分析。采用DEGseq v1.18.0 软件[20]进行基因差异表达分析，界定犏牛新鲜囊胚和冻融囊胚二者文库间的DEGs。将DEGs 向GO(Gene Ontology)数据库(http://www.geneontology.org/)各个条目进行映射，以校正后P＜0.05 为阈值，找出DEGs 中显著富集的GO 条目。然后通过与KEGG(Kyoto Encyclopedia of Genes and Genomes)数据库进行比对，分析DEG s 涉及的信号通路(Pathways)。采用SAMtools-0.1.19[21]和ASprofile v1.0.4 软件[22]分别对两样本数据进行SNP 和AS 分析。

1.2.5 测序结果的实时荧光定量PCR验证 为验证测序结果的准确性，采用实时荧光定量PCR方法验证基因的表达情况，其中随机选取4个DEGs，并以H2A 作为内参基因。用Primer Premier 5软件进行引物设计，引物序列和产物长度如表1所示。PCR 总反应体系为10 μL，包括上下游引物各0.8 μL、cDNA 模板0.5 μL、ddH2O 2.9 μL和Sso AdvancedTMSYBR®Green Super mix 5 μL。PCR 扩增程序：首先95 ℃预变性3 min；然后95 ℃变性10 s，60 ℃退火20 s，72 ℃延伸60 s，30个循环；最后72 ℃延伸5 min，并4 ℃保存。

表1 实时荧光定量PCR 引物信息Table 1 Primers for real-time fluorescence quantitative PCR

2 结果

2.1 测序产量统计

新鲜囊胚组和冻融囊胚组分别获得52641524条和56521104 条Raw Reads，过滤后分别获得51099116 条和54192358 条Clean Reads，分别占Raw Reads 的97.07%和95.88%。Q30 可用于反映测序的碱基质量水平，质量值(Q)越高代表碱基被测错的概率越小。本次新鲜囊胚组和冻融囊胚组测序的Q30 百分比分别为92.02%和91.91%，表明此次测序质量水平高，文库构建完整，可进行后续数据分析。

2.2 测序饱和度分析和随机性分析

如图1A 和图1B 所示，将数据按照10%、20%......100%比例进行测序，并统计在每个测序数据量下获得的基因的表达量。根据基因的表达量不同，将基因分为4个等级，统计不同数据量下基因定量的准确性。横坐标为不同深度的数据量，纵坐标为在该深度下准确估计表达量的基因百分比。犏牛新鲜囊胚和冻融囊胚两样本检测到的基因数随着测序数据量的增多而升高，当测序量比例达到30%时，可明显看出测序数据量已达平台期，增长速度已基本趋于平缓，表明在此时检测到的基因数已基本达到饱和值，两样本的测序量已经基本覆盖到细胞中表达的全部基因，测序结果准确可靠。

测序结果随机性(均一性)则是通过检测统计每个基因在不同位置的深度来衡量确定的。因为基因长度各不相同，所以将基因外显子区域平均划分成100个窗口，统计所有基因在该窗口的深度平均值。如图1C 和图1D 所示，两样本的打断随机性(即均一性)良好，Reads在基因各等分的分布较为均匀，不偏向于基因的特定区域，说明数据可靠，可用于后续数据分析。

图1 测序饱和度(A 和B)和随机性(C 和D)分析Figure 1 Sequencing saturation (A and B) and randomness (C and D) analysis

2.3 测序结果比对分析

运用TopHat软件将犏牛新鲜囊胚和冻融囊胚两样本的Clean Reads 与牦牛参考基因组的Reads 进行比对，发现其匹配总数比例分别为80.00%和80.19%，说明该样品无外源物种污染。

2.4 DEGs 的GO 功能分析

GO 注释分类后显示，犏牛新鲜囊胚和冻融囊胚中共有10532个DEGs，被归类于23个生物过程(biological process，BP)、22个细胞组成(cellular component，CC)和22个分子功能(molecular function，MF)为主的3个GO 分类条目中(图2)。在BP 分类中，细胞过程类别的富集最显著；在CC 分类中，细胞组成类别的富集最显著；在MF 分类中，整合类别的富集最显著。

2.5 DEGs的KEGG通路分析

对犏牛新鲜囊胚和冻融囊胚的DEGs 进行KEGG 信号通路分析，结果表明，共有12315个DEGs 注释到318条信号通路。显著富集的前5条通路(Top 5 pathways)如表2 所示，其中剪接体通路的富集程度最高，且犏牛冻融囊胚相比新鲜囊胚的大部分DEGs 表达水平上调。

2.6 基因变异分析

采用SAMtools-0.1.19软件进行基因变异检测结果分析，结果显示，犏牛新鲜囊胚和冻融囊胚两样本分别检测到116681 和224750个SNP 位点，6668 和10729个插入/缺失(InDel)位点，总变异个数为123349 和235479个。犏牛冻融囊胚样本中SNP位点数和InDel位点数均明显多于新鲜囊胚样本。

图2 差异表达基因(DEGs)的GO 功能分类Figure 2 GO functional classifications of differentially expressed genes (DEGs)

表2 差异表达基因(DEGs)富集前5 位的KEGG 通路Table 2 Top 5 of enriched KEGG pathways of differentially expressed genes (DEGs)

根据碱基替换方式不同，SNP位点可分为转换型(transition)和颠换型(transversion)。根据犏牛新鲜囊胚和冻融囊胚两样本检测到的SNP 位点，统计每个突变类型的频率(图3)。在新鲜囊胚样本中，发生频率最高的是转换型SNP(A/G、C/T)位点，约占总位点数的72%，而颠换型SNP(A/C、A/T、C/G、G/T)位点约占28%，转换和颠换的比例(Ts/Tv)为2.57。在冻融囊胚中，发生频率最高的是转换型SNP 位点，约占总位点数的71%，颠换型SNP 位点约占29%，转换和颠换的比例为2.45。在新鲜囊胚样本的转换型SNP 中，C/T 类型数量略多于A/G 类型数量；但在冻融囊胚样本中，A/G 类型数量略多于C/T 类型数量。两样本中C/T 类型与A/G 类型所占比例都无较大差异，且颠换型SNP 中A/T 所占比例最少。

图3 单核苷酸多态性(SNP)突变频率分布Figure 3 Mutation frequency of single nucleotide polymorphism (SNP)

2.7 基因的AS 分析

AS 广泛存在于真核生物中，是调节基因表达和蛋白质多样性的重要机制。AS 事件主要分为外显子跳跃和盒式外显子跳跃、内含子滞留和多重内含子滞留、可变5’端或3’端剪切、转录起始区域可变剪切(alternative splicing of transcription start sites，TSS)、转录结束区域可变剪切(alternative splicing of transcription termination sites，TTS)共5 种类型[22](图4A)。

图4 可变剪切事件类型(A)和分析统计图(B)Figure 4 The types (A) and analysis statistics (B) of the alternative splicing (AS)

用ASprofile 软件分析各样本基因的AS 事件(图4B)。从犏牛新鲜囊胚和冻融囊胚中分别检测出49388 和65241个AS 事件。在两个样本检测到的AS 事件类型中，TSS 是最普遍的类型，新鲜囊胚中TSS 数量为22962个，占总AS 事件的46.5%；冻融囊胚中TSS 数量为28652个，占总AS 事件的43.9%。其次是TTS 类型，新鲜囊胚中TTS数量为22048个，占总AS事件的44.6%；冻融囊胚中TTS 数量为27196个，占总AS 事件的41.7%。无论从总数还是各类型来看，犏牛冻融囊胚的AS 事件均多于犏牛新鲜囊胚。

2.8 测序结果的实时荧光定量PCR 验证

随机选取HSP110、RPL31、REX1 和BAD共4个DEGs，采用实时荧光定量PCR 方法验证这些基因的表达情况，并采用Pfaffl 法[23]分析4个基因在两样本中的相对表达量。结果表明，实时荧光定量PCR和RNA-Seq两种方法计算得到的差异倍数基本一致(图5)。PCR 验证说明此次RNASeq 结果准确可靠。

图5 实时荧光定量PCR 验证RNA-Seq 检测出的差异表达基因(DEGs)Figure 5 qRT-PCR to verify the differentially expressed genes (DEGs) detected by RNA-Seq

3 讨论

本研究采用SAMtools-0.1.19对犏牛新鲜囊胚和冻融囊胚两样本测序数据进行SNP 检测分析，结果显示犏牛冻融囊胚样本中SNP位点数和InDel位点数都明显多于犏牛新鲜囊胚样本，说明胚胎经玻璃化冷冻后，其基因突变的概率明显增大。犏牛新鲜囊胚和冻融囊胚两样本中转换型SNP比例都明显高于颠换型，Ts/Tv 比值分别为2.57 和2.45。Ts/Tv 理论比值是0.5，而本研究中Ts/Tv比值要远大于理论值，这与Kraus 等[24]和Aslam等[25]在其他哺乳动物中的研究结果相似；这一现象称之为“转换偏差”。转换偏差这种现象的产生，可能是因为物种长期自然进化的选择，有利于保持编码蛋白原有结构，以减少有害突变形成的概率[26-27]。也有人认为，嘌呤和嘧啶的结构和代谢等内在特征决定，才导致转换偏差的产生[28]。本研究中，犏牛新鲜囊胚样本中C/T 类型数量略多于A/G 类型。这是因为DNA 中的5-甲基胞嘧啶(5-methylcytosine，5mC)有较高的频率突变成胸腺嘧啶(thymine，T)，所以一般在SNP 中C/T类型所占比例最高[29-30]。然而犏牛冻融囊胚样本中A/G 类型数量略多于C/T 类型，这与张得芳等[31]利用高通量测序技术对不同海拔下花叶海棠叶片进行转录组测序研究的结果相似，推测可能是样本自身GC 含量较高所致。犏牛囊胚经冷冻操作后，不管是SNP 和InDel 位点数明显增多，还是发生转换偏差，都提示哺乳动物胚胎经超低温冷冻保存后突变概率明显增大，为有效降低冷冻程序带来的损伤，SNP在胚胎玻璃化冷冻过程中起重要的调控作用。

本研究从犏牛新鲜囊胚和冻融囊胚两样本中分别检测出49388 和和65241个AS 事件，说明二者含有大量的AS 基因。此外，犏牛新鲜囊胚和冻融囊胚发生AS 的数量存在差异，表明新鲜囊胚经历胚胎冷冻和复苏等外界环境的变化后，AS 基因的表达存在特异性。研究表明，除了不同环境下AS 存在基因特异性外，生物不同发育时期[32]及不同组织[33]中均可发现AS 基因特异性。Montelli 等[34]研究发现，当生物处于不同发育时期或在不同组织结构中，AS 能有效调控基因表达产物的功能及结构。这些基因的AS 类型非常丰富，目前共鉴别出12 种类型，其中以TSS、TTS、外显子跳跃和盒式外显子跳跃、内含子滞留和多重内含子滞留、可变5’端或3’端剪切类型为主。本研究检测到的AS 类型中，TSS 比例最高，在犏牛新鲜囊胚中发生率为46.5%，犏牛冻融囊胚中发生率为44.6%；其次是TTS，犏牛新鲜囊胚和冻融囊胚中的发生率分别为43.9%和41.7%。这与字向东等[10]利用RNA-Seq 技术对牦牛IVF 胚胎发育转录组的检测结果一致。但是，另有学者通过高通量测序技术分析犏牛睾丸转录组时发现，内含子滞留事件发生率最高(40%)，其次是外显子跳跃事件(23%)[35]。也有文献表明，猪[36]和京海鸡卵巢[37]中内含子滞留事件发生率最低。有研究者认为，内含子滞留事件在动物中发生率低，而在植物中最常见[38]。这些研究表明，在不同生物、不同发育时期[32]、不同组织[33]间所发生的AS类型均存在差异性，这可能是由在不同生物、不同组织和不同生长发育时期真核生物的基因和蛋白质数量差异较大所导致。

对犏牛新鲜囊胚和冻融囊胚两样本的DEGs进行KEGG 通路分析发现，剪接体通路富集程度最高，也进一步证明了胚胎经超低温冷冻保存后发生了大量AS 事件。在剪接体通路中犏牛冻融囊胚上调表达最为显著的DEGs 是以mRNA 前体剪接因子(pre-mRNA processing factor，Prp)4 和Prp19 为代表的大部分mRNA 前体剪接因子。李昂等[39]研究发现，Prp19 会对HeLa 等细胞产生不良反应，导致细胞凋亡，并可能会激活前体RNA剪接体，从而引发DNA 损伤效应[40]。Prp4 可通过影响一些基因的表达，间接使其内含子剪切效率明显降低，从而可能发生DNA 突变[41]。而内含子剪切效率的降低，还会引起mRNA 降解，降低蛋白质的表达水平，从而影响蛋白质的功能。所以，从犏牛新鲜囊胚和冻融囊胚两样本检测出大量AS 事件，以及通过KEGG 通路分析发现剪接体通路的显著富集，这些结果都表明，胚胎经超低温冷冻保存后发生了大量AS 事件，而一些Prp的上调可能会对胚胎产生不利影响。

有关人类胚胎的临床研究表明，冻融胚胎的发育情况、临床妊娠率及活产率与新鲜胚胎相比无明显差异，甚至更优于新鲜胚胎[42-43]。临床上大多数是以形态学为基础判断冻融胚胎发育情况是否良好，并在临床妊娠率、活产率等临床指标上判断冻融胚胎和新鲜胚胎的差异是否显著。从临床指标上一般分析认为冻融胚胎和新鲜胚胎并未存在显著差异，但这并不能证明新鲜胚胎经超低温冷冻保存后不会造成一定程度的冷冻损伤。有研究指出，冻融胚胎和新鲜胚胎的临床妊娠率和活产率无明显差异，甚至冻融胚胎更优于新鲜胚胎，这可能是由于在新鲜周期和过度刺激周期中子宫内膜的容受性受损，而并非冻融胚胎质量优于新鲜胚胎[44]，这也从另一面与本研究相契合，证明胚胎经超低温冷冻保存后，并不会改善新鲜胚胎的质量，反而会不可避免地造成一定程度的冷冻损伤。