王进娥，麻俊渊，陈舟，李书元，李洋，霍生东，杨艳梅，扎西英派

（西北民族大学生命科学与工程学院，甘肃 兰州 730030）

牦牛生活在青藏高原及其毗邻的高山、亚高山地区，具有极强的高原适应性和耐粗饲等特点，有“高原之舟”的美誉，是宝贵的遗传基因库［1］。发情季节主要集中在气候温暖、牧草丰盛的暖季（7 月～10月），发情周期为16～25 d，妊娠期为250～260 d，产犊期在每年的4～7 月份［2-3］。由于受光照、群体、气温、海拔、营养、基因等影响，牦牛性成熟较慢，生殖能力低下。一头成年牦牛的平均繁殖率仅为48.61%，其中超过一半为两年一胎或三年两胎［4］。一般情况下，雌性牦牛发情起止时间平均为36 h，变动范围为24～48 h［5］。妊娠期牦牛的卵巢通常不规则，卵巢表面出现妊娠黄体维持妊娠［6］。卵巢黄体从妊娠中期的36 个月开始下降，但组织结构仍相对紧密［7］。

牛津纳米孔（nanopore）技术是基于纳米孔单分子实时电信号的第三代测序技术，其原理是纳米孔内有共价结合的分子接头，当单个碱基或DNA分子通过纳米孔通道时，会使电荷发生变化，从而短暂地影响流过纳米孔的电流强度。由于化学结构的差异，A、C、G、T这4种不同碱基通过纳米孔时会产生不同强度的电流，通过灵敏的电子设备可以检测到电流变化，进而可以识别DNA 链上的碱基完成测序［8］。自2012 年牦牛全基因组序列测定［9］，转录组学技术在牦牛繁殖调控领域已经取得了较多突破性的研究成果。兰道亮等［10］基于RNA-Seq 技术对牦牛卵巢进行转录组测序分析发现：6321个新转录本，预测2267 个新转录本具有编码蛋白能力，新转录本GO分类注释显示与繁殖性状相关的新转录本最多。为深入了解牦牛繁殖机能提供了基础数据；字向东等［11］为探明不同发育阶段牦牛胚胎的转录组差异，对IVF 生产的牦牛2-细胞、4-细胞、8-细胞、桑椹胚和囊胚5 个发育阶段的胚胎进行RNA-Seq 转录组测序，比较分析发现：胚胎转录本中转录起始区域可变剪接和转录结束区域可变剪接所占比例最大；蒲思颖等［12］利用RNA-seq 技术分析了囊胚玻璃化冷冻前后转录组变化，结果发现：在冷冻前后牦牛囊胚中分别检测到9827 和13567 个转录本，这些研究揭示了牦牛胚胎细胞不同发育阶段差异表达基因状况，为牦牛发育分子调控机制及完善胚胎体外培养技术奠定了基础。然而，目前利用第三代测序技术对牦牛妊娠期方面尚未见研究成果，本研究填补了该技术在牦牛妊娠期的研究空白。

在牦牛的发情周期中，机体在发情期和妊娠期的调控是不同的，但都处于一个完整的发情周期内，在时间进程上是有联系的，因此利用第三代测序技术，对牦牛发情期与妊娠期卵巢组织的转录组进行测序和分析，揭示牦牛发情期与妊娠期卵巢组织差异表达转录本及发生的AS 和APA 的数量，获得差异表达转录本的功能、分类和相关信号通路，对现有的牦牛基因组和转录组注释信息进行完善，为开展牦牛发情期与妊娠期相关基因的研究及分子调控机制研究奠定基础。

1 材料与方法

1.1 试验动物

2019年9月，在青海省海盐县牧民牦牛群中（海拔约3 200 m）选取16 头6 岁左右母牦牛，在自然放牧条件下每天早、晚2次跟群观察母牦牛发情表现，每次3 h以上。每次跟群观察期间，母牦牛会出现2次及以上试情公牛爬跨时呆立不动的“静立反射”发情现象。同时母牦牛小便频繁，阴道流出的粘液较多，甚至有少量血液流出，有以上现象的牦牛确定处于发情期。以发情牦牛为对照组，直到2019 年10月，连续两次直肠检查后，卵巢表面不规则，一侧有妊娠黄体出现，确定为妊娠期。选择发情期和妊娠期牦牛各3头，屠宰前精神、营养状况良好，经临床常规检查，确认没有可见疾病，屠宰后，迅速取出2侧卵巢，用生理盐水清洗，且发情期与妊娠期3头牦牛卵巢分开标记，并储存在液氮中，送北京百迈克生物科技有限公司测序，公司将发情期与妊娠期每头牦牛的2 侧卵巢混合提取RNA，然后将每头牦牛提取的RNA进行单独测序。

1.2 差异表达转录本（DET）筛选

在差异表达转录本检测过程中，将发情期和妊娠期样品转录本进行两两比较，通过对差异显著性P值进行校正获得错误发现率（false discovery rate，FDR）。两组样品转录本表达差异倍数大于等于2（FC≥2），且FDR<0.01视为差异转录本筛选标准。筛选得到的差异转录本采用GO（gene ontology）和KEGG（kyoto encyclo-pedia of genes and genomes）进行深入分析。

1.3 基因的AS分析

可变剪切（alternative splicing，AS）是指一个基因转录形成的前体mRNA（pre-mRNA）通过不同的剪切方式产生不同成熟mRNA 的过程［13］。AS的类型主要包括可变转录终止位点（alternative3’splicesite，A3），可变转录起始位点（alternative 5’splicesite，A5），外显子跳跃（skipping exon，SE），内含子保留（retained intron，RI），可变外显子（mutually exclu⁃sive exons，MEE）。利用Astalavista 软件鉴定牦牛发情期与妊娠期基因的AS 类型，采用软件默认参数。

1.4 基因的可变多聚腺苷酸化分析

可变多聚腺苷酸化（alternative polyadenylation，APA）是指通过poly-A尾巴的AS，源于同一基因的不同转录本含有不同长度的3’非翻译区（untranslated region，UTR）。通过TAPIS pipeline 对牦牛发情期与妊娠期基因的APA位点进行鉴定，采用默认参数。

2 结果与分析

2.1 差异基因的鉴定

与发情期相比，妊娠期牦牛卵巢组织中共鉴定得到差异表达基因有1 636个，妊娠期较发情期上调基因有746个，下调基因有890个（图1）。

图1 牦牛妊娠期卵巢差异表达基因数目统计Figure 1 Statistical result of the number of differentially expressed genes in the ovary of yak during pregnancy

2.2 差异转录本分析

与发情期相比，妊娠期牦牛卵巢组织中共鉴定得到差异表达转录本3 562个，妊娠期较发情期上调转录本有1 693个，下调转录本有1 896个（图2）。

图2 牦牛妊娠期卵巢差异表达转录本数目统计Figure 2 Statistical result of the number of differentially expressed transcripts in the ovaries of yaks during pregnancy

2.2.1 差异表达转录本GO 注释及富集 根据比较分析牦牛发情期和妊娠期差异表达转录本的表达量，对差异表达转录本进行了GO注释分析。注释到GO 数据库共有56 个功能条目，发现差异表达转录本在生物学过程、细胞组分及分子功能三大类中的较多GO类别上均有注释。在生物过程的分类中共有21 个功能条目在GO类别出现富集。注释数量最多的是细胞过程（cellular process）、单一组织过程（single-organism process）和 生 物 调 节（biological regulation）。在牦牛繁殖方面，187 个DETs与繁殖、繁殖过程和节律过程有关。繁殖与繁殖过程均包括88个DETs，其中60个上调，28个下调；节律过程中有11个DETs，其中7个上调，4个下调（图3-A）。在细胞组分（cellular component）分类中共有18个功能条目在GO 类别出现富集。主要集中于细胞（cell）、细胞组件（cell part）和细胞器（organelle）（图3-B）；在分子功能（molecular function）分类中共有17 个功能条目在GO类别出现富集。有1419 个DETs与结合（binding）分子功能类别相关，其中858 个上调，561 个下调；有774 个DETs 与催化活性（ catalytic activity）相关，其中450 个上调，324 个下调（图3-C）。

图3 牦牛发情期与妊娠期差异表达转录本GO注释Figure 3 Go annotation of the transcript differentially expressed during estrus and pregnancy in yak

2.2.2 差异表达转录本KEGG 注释及富集 由图4 可知，在KEGG 数据库中鉴定到1 558 个DETs。在牦牛发情期和妊娠期卵巢中203 条通路注释到KEGG 数据库，主要与细胞过程、环境信息加工、遗传信息处理、代谢和生物系统有关。其中代谢类最多（11 个亚类），其次是生物系统类（10 个亚类）。相比之下，与遗传信息处理相关的较少（6个亚类）。细胞过程类包括4 个亚类。此外，环境信息过程的3 个亚类途径相关。其中PI3K-Akt信号通路、细胞内吞作用、剪切体和氧化磷酸化是关键通路（表1）。牦牛繁殖相关通路还有MAPK信号通路、雌激素信号通路（estrogen signaling pathway）、孕酮介导卵母细 胞 成 熟 （progesterone-mediated oocyte maturation）、细胞周期（cell cycle）、催产素信号通路（oxytocin signaling pathway）和卵母细胞减数分裂（oocyte meiosis）等。

表1 部分差异表达转录本富集KEGG通路Table 1 Section of enriched KEGG pathways of DETs

图4 牦牛发情期与妊娠期差异表达转录本KEGG注释图Figure 4 KEGG annotation of the transcript differentially expressed during estrus and pregnancy in yak

2.3 牦牛发情期与妊娠期基因的AS数量及类型

通过对牦牛发情期与妊娠期卵巢样品进行AS事件比较，在发情期中共鉴定到2 479次AS事件，包括1 214（48.97%）次SE，433（17.47%） 次A3，401（16.18%） 次RI，358（14.44%）次A5 和73（2.94%）次MEE（图5-A）。在妊娠期中共鉴定到2 053 次AS事件，包括853（41.55%） 次SE，528（25.72%） 次RI，341（16.61%） 次A3，289（14.08%） 次A5 和42（2.04%） 次MEE。SE 为主要剪切类型，MEE 为稀有剪切类型（图5-B）。

图5 牦牛发情期与妊娠期AS事件类型Figure 5 Pie chart of variable splice event types

2.4 牦牛发情期与妊娠期基因的APA位点数量及特征

牦牛发情期共鉴定到9 189 个基因含有1 个及以上的APA位点，其中1 个APA位点的基因数量最多，达到3 701（39.9%）个；2和3 个APA位点的基因数量次之，分别为1 873（20.33%）个和1 166（12.67%）个；另有1 218（13.25%） 个基因有5 个以上的APA位点（图6-A）；牦牛妊娠期共鉴定到8 924个基因1个及以上的APA位点，其中1个APA位点的基因数量最多，达到3 908（43.73%）个；2 和3 个APA 位点的基因数量次之，分别为2 006（22.43%）个和1 077（12.03%）个；另有829（9.29%）个基因有5个以上的APA位点（图6-B）。

图6 牦牛发情期与妊娠期基因APA位点个数分布Figure 6 Number distribution of APA loci

3 讨论

目前，基于第3 代测序技术AS 和APA 研究少且主要集中在人类和少数模式物种［14］。对于牦牛发情期与妊娠期基因AS 和APA 的相关信息几乎空白。本研究结果表明，牦牛发情期与妊娠期AS类型均以SE 为主，其次为A3 和RI，MEE 最少。这与John Wiley等［15］研究结果AS事件中，植物中多以RI为主导，动物中多以SE 占主导是一致的。本研究中，在牦牛发情期基因中鉴定到的AS 事件数量（2 479 次）多于妊娠期基因中的AS事件数量（2 053次），此外牦牛发情期基因中发生APA 的基因数量（9 189 个）稍多于妊娠期基因中发生APA的基因数量（8 924 个），说明牦牛发情期中基因的AS和APA更为活跃。鉴定到的AS事件和APA位点信息可进一步完善牦牛发情期和妊娠期基因组和转录组注释信息，也为其他动物的相关研究提供可参考信息。

3.1 周期节律与繁殖过程

在本研究中，GO 分析结果显示，节律过程涉及由NRIP1 上调的DETs，由SOD 下调的DETs。核受体相互作用蛋白1（NRIP1），也称为受体相互作用蛋白140（RIP140），是通过与核受体、转录因子和其他辅助调节因子相互作用调节基因表达的关键因子［16］。RIP140 通过与RORS 的相互作用参与生物钟基因的调节［17］。赵梦静等［18］发现肝脏RIP140 的昼夜节律性，RIP140上调Cyp2b10在小鼠和AML-12 细胞中的表达，并减弱其节律。超氧化物歧化酶（SOD）在体内主要是清除超氧阴离子自由基。有研究发现，成年人24 h血浆SOD参数变化呈现明显的近日节律特征，其峰值分别出现在晚上和下午［19］。这都进一步为NRIP1、SOD分别参与牦牛发情期与妊娠期节律上调与下调提供了依据，但其调控机制尚未清楚，需进一步探究。周期节律的变化对牦牛的繁殖过程有重要影响［20］。光照季节性变化决定动物是否发情，由此将季节性繁殖动物分为长日照和短日照动物［21］。光照周期的模式、光照周期变化的趋势以及光照周期信号表达的节律的周期性都影响着动物的繁殖活动［22-23］。

3.2 激素与繁殖过程

牦牛是季节性繁殖动物。发情季节主要集中在7月～10月［2］。本研究中，GO分析结果显示，在生物过程组分中，繁殖与繁殖过程均有88 个DETs，其中60 个上调，28 个下调。热休克蛋白（heat shock pro⁃tein，HSP）是ATP 依赖的高度保守的分子伴侣，能利用ATP水解产生的能量参与蛋白质的正确折叠，从而在细胞周期运行，细胞增殖、迁移、凋亡等过程发挥重要作用。HSP90能与醛固酮、雄激素、雌激素和孕激素等类固醇激素的受体结合，在促性腺激素活性的调节中发挥重要作用。有研究发现，HSP90参与胚胎的发育过程，核移植时胚胎的低生存率与HSP90的低表达有关系，将HSP90基因导入到核移植胚胎中可提高胚胎的生存率［24］。在本研究生殖过程中，HSP90基因上调，且牦牛妊娠期表达含量高于发情期，提示HSP90 上调可提高胚胎生存率。在KEGG 分析中，HSP90 在PI3K/Akt信号通路、孕酮介导卵母细胞成熟和雌激素信号通路中上调，这与GO分析结果一致。HSP70首先是在胚胎发生中被鉴定且起负调节维持激素-受体的平衡。有研究发现，过量的HSP70对胚胎发生产生不利影响。据报道妊娠早期外周血淋巴细胞中HSP70 的细胞内浓度与自然流产或胎儿死亡等不良妊娠结局直接相关［25］。本研究中HSP70在雌激素信号通路中下调，表明有利于维持妊娠。P38在MAPK信号通路中上调，p38MAPK的磷酸化现象存在于大多数哺乳动物的整个妊娠阶段，且随着妊娠的进展这种现象逐渐加强［26］。PI3K/Akt 信号通路成分和一些下游效应分子存在于卵泡中，在形成、生长、排卵和黄体化过程中发生变化［27］。cAMP可以抑制卵母细胞囊泡破裂（oocyte vesicle rupture，GVBD）的发生，维持卵母细胞减数分裂的停止，而维持减数分裂停止所需要的cAMP 是由卵母细胞自身产生的［28］。Cdc2 在孕激素介导的卵母细胞成熟过程中上调。细胞分裂周期基因（cell division cycle，cdc基因）控制着细胞周期的启动及各时相的转换，其中以Cdc2 基因最重要，该基因是唯一在两个控制点Gl→S、G2→M 均起作用的cdc 基因［29］。P34Cdc2激酶的活化是细胞分裂的增殖信号，它具有启动DNA复制和诱发有丝分裂的双重作用，在细胞周期调节中起重要作用［29］。

3.3 氧化磷酸化与牦牛繁殖过程

氧化磷酸化是线粒体中的物质在氧化时释放能量供给ADP 和无机磷合成ATP 的偶联反应，是机体重要的能量代谢途径。本研究氧化磷酸化通路中的多个基因表达下调，主要集中在NADH脱氢酶亚基 家 族 成 员（NADH，ND1，ND2，ND5，ND6，Ndufa7，Ndufb2，Ndufb5）和细胞色素c 氧化酶家族成员（cytochrome c oxidase，COX1，COX2，COX3，COX7A）。目前这些基因家族在牦牛发情期与妊娠期卵母细胞发育成熟过程中的相关作用未见报道，这些基因家族均为线粒体氧化磷酸化通路中的关键基因，本身发挥能量代谢作用，这些基因在氧化磷酸化通路中的下调，推测其功能可能与降低相关能量代谢，阻止进一步发育，维持发情期的分裂阻滞状态，同时减小氧化磷酸化供能反应给细胞带来的氧化损伤起作用。

3.4 AS事件与牦牛卵巢发育

KEGG 通路分析还发现，剪接体通路富集程度较高，进一步证明了牦牛发情期与妊娠期转录本发生了大量AS 事件。在剪接体通路中表达上调最为显著的DETs 是mRNA 前体剪接因子（pre-mRNA processing factor，Prp）19 和43 为 代 表 的 大 部 分mRNA 前体剪接因子。李昂等［30］研究发现，Prp19对HeLa 等细胞产生不良反应，导致细胞凋亡，并会激活前体RNA 剪接体，因而出现DNA 损伤效应。Prp43 及其辅助因子突变抑制剪接体激活缺陷的剪接因子突变，通过影响一些基因的表达，间接使其内含子剪切效率降低，引起DNA 突变［31］，内含子剪切效率的降低，使得mRNA 降解，蛋白质表达水平降低，从而影响蛋白质的功能。所以，从牦牛发情期与妊娠期检测出大量AS 事件，以及通过KEGG 通路分析发现剪接体通路的显著富集，这些结果都表明，牦牛在发情与妊娠期间发生了大量AS事件，而一些Prp的上调可能会对卵巢发育产生不利影响。

4 结论

综上所述，基于第3代测序技术对牦牛发情期与妊娠期卵巢组织的转录组进行测序和分析，揭示出牦牛不同生理状态下卵巢组织差异表达转录本的数量，获得了差异表达转录本的功能、分类和相关信号通路。为丰富牦牛卵巢组织转录组信息，展牦牛发情期与妊娠期相关基因的研究及分子调控机制奠定了基础。