牛卵泡发育相关下调基因的筛选
2020-10-14赵园园赵成刚安清明武渊勋孟金柱
赵园园，赵成刚，安清明，武渊勋，孟金柱
(1.铜仁学院，贵州 铜仁 554300；2.山西普源泰奶牛养殖有限公司，山西 晋中 030800)
在牛和人等单胎动物中，有腔卵泡的生长呈典型“波”的模式[1-3]。在每个卵泡波开始之前，血清中促卵泡素(FSH)会短暂升高，并刺激一群小的有腔卵泡出现。其中，优势卵泡(DF)在FSH浓度下降的情况下会继续生长到排卵卵泡大小，同时产生大量的雌二醇(E2)，其余较小的一些从属卵泡(SF)则会失去产生雌二醇的能力，从而导致卵泡闭锁[4-6]。单胎动物卵泡的选择是一种进化上的保守机制，对控制每次妊娠的后代数量至关重要[7-8]。尽管垂体促性腺激素在调节卵泡发育波中的关键作用已被充分证实[9-12]，但在每个卵泡波中，最终只选择一个卵泡成为DF的分子调控机制尚不清楚。
转录组测序技术已被应用于探索多种物种的转录调控机制，如奶牛[13]、绵羊[14]和猪[15]。LI等[16]通过对水牛不同大小卵泡(直径<5 mm、5～8 mm、8～12 mm及>12 mm)颗粒细胞(GCs)进行转录组测序，认为卵泡的成熟和排卵可能受到免疫控制，包括同种异体排斥反应、趋化因子信号通路、自然杀伤细胞介导的细胞毒性、吞噬体、抗原处理和呈递。对牛出现偏差前的第一大卵泡(PDF1)和出现偏差时的最大卵泡(ODF1)GCs进行转录组测序得到41个表达下调基因，42个表达上调基因，MAPK13、CPXM1和ARID4B在卵泡发育过程中可能起抑制作用[17]，但对牛卵泡闭锁调控机制的研究较少。因此，本研究经过对牛DF和SF中的GCs进行深度测序，筛选出卵泡发育下调相关基因，对分析牛卵泡GCs凋亡的关键调控因子具有重要意义，进一步为阐明单胎动物卵泡闭锁的调控机理奠定基础。
1 材料与方法
1.1 实验动物样品采集在山西普源泰奶牛养殖有限公司(山西晋中)，选取6头待淘汰的健康荷斯坦奶牛送往屠宰场屠宰后，分别采集每头牛卵巢上正常发育的最大卵泡(8～10 mm)与第二大卵泡(5～8 mm)，作为候选DF和SF置于4℃杜氏磷酸缓冲液中迅速带回实验室。
1.2 卵泡液激素测定采集到的每头牛候选DF和SF分别放到6个盛有生理盐水的培养皿中，用注射器抽取卵泡液后，分别用E2ELISA试剂盒和PROG ELISA试剂盒(上海蓝基)测定雌激素(E2)和孕激素(P)。
1.3 总RNA提取、文库构建及测序将获得的候选DF卵泡液中E2/P>1，候选SF卵泡液中E2/P<1所对应的3对卵泡用眼科剪剪开并用细胞刮刮取GCs。获得的GCs经RNA提取、纯化、完整性检测、浓度测量后，交北京诺和致源生物信息科技有限公司进行文库构建、上机测序。
1.4 数据处理及分析用FastQC软件将获得的原始reads(raw reads)进行质量评估，并经过一系列处理(去除质量低的、含N的和带接头的序列)从而获得高质量的干净reads(clean reads)。将获得的干净序列通过从头组装得到(singleton)和重叠群(contigs)，此时得到unigenes序列。通过SOAP V2.0软件将unigenes序列与牛的参考基因组数据库进行比对。对获得的mRNA，使用DESeq2差异表达分析，同时筛选出牛卵泡发育相关表达下调基因。DAVID软件对获得的表达下调基因进行GO及KEGG信号通路分析。
1.5 反转录及引物合成使用北京全式金生物技术有限公司的反转录试剂盒将1.3得到的总RNA反转录成cDNA，具体反应条件参照说明书：42℃，15 min；85℃，5 s 。通过Primer 5.0将筛选出具有代表性的表达下调4个基因(PPP1R14A、OLA1、QRFPR和RMRP)进行引物设计(表1)，内参采用β-actin，并交北京华大基因有限公司合成。
表1 供试引物序列
2 结果
2.1 牛优势卵泡与从属卵泡筛选结果经过对候选DF和SF卵泡液中E2和P 的测定，发现3对卵泡卵泡液中E2/P>1，E2/P<1，确定其为DF和SF，可以采用这3头牛的卵泡GCs进行测序(表2)。
2.2 牛卵泡GCs中总RNA原始数据过滤情况分析通过深度测序共获得93 304 489条原始reads，有质量低的、含N的和带接头的reads，为了保证信息分析质量，经过FastQC软件处理得到92 019 573条高品质干净的reads，详细结果见表3。
表2 DF和SF卵泡液E2和P的测定结果
表3 牛卵泡颗粒细胞中总RNA原始reads数据过滤情况
2.3 牛卵GCs中差异表达基因筛选将测序结果与牛的参考基因组数据库进行比对，共获得32 346个基因， 将DF和SF的FPKM进行标准化后，发现FPKM≥0.5 的基因13 243个，表4列出了表达量最高的10个基因及其功能。通过DESeq2软件对13 243个基因进行差异表达分析，当P<0.05时共得到696个差异表达基因，其中194个在DF中表达显著上调，502个表达下调(图1)。
2.4 牛卵GCs中表达下调基因GO及KEGG分析通过对牛卵泡GCs中表达下调基因进行GO功能分析，对可能抑制卵泡发育的基因进行分类。应用DAVID软件对获得的502个表达下调基因进行GO分析,共分为3大类104组，其中参与生物学过程(Biological process)占61.5%，参与细胞内信号级联的基因数居首位(15个)；与细胞组分有关基因占25.0%,其中参与细胞内无膜界细胞器和无膜界细胞器的基因最多(分别为30个)；与分子功能相关的基因占13.5%，其中参与钙离子结合的基因最多(17个)(图2)。
表4 牛DF 和 SF中表达量最高的10个基因
图1 DF和SF中差异表达基因火山图 横坐标代表mRNA在DF/SF中表达倍数变化；纵坐标代表mRNA表达量变化的统计学显著程度；散点代表各个mRNA；蓝色圆点表示无显著性差异的mRNA；红色圆点表示显著下调的差异mRNA；绿色圆点表示显著上调的差异mRNA
为了进一步获取与抑制牛卵泡发育相关的信号通路，再次使用DAVID软件对表达下调基因进行KEGG信号通路分析，共发现5条通路(表5)，其中基因富集最为显著的是癌症通路。
2.5 QRT-PCR分析通过深度测序集合功能分析，最后从中筛选出4个具有代表性的在牛卵泡发育过程中可能会起抑制作用的基因(表6)。QRT-PCR结果显示，EGR1、PPP1R14A和QRFPR在DF和SF中的表达趋势与深度测序结果一致，且EGR1在SF中的表达量极显著地高于DF (P<0.01)；PPP1R14A和QRFPR在SF中的表达量显著地高于DF (P<0.05)，OLA1在SF和DF中的表达趋势与深度测序结果相反，但差异不显著(图3)。
3 讨论
哺乳动物在发育过程中，卵巢上最大的卵泡，即未来的DF持续生长、成熟并排出卵子，而未来的SF由于生长速率和雌二醇分泌降低，从而走向闭锁[18-21]。研究发现，许多细胞死亡配体(包括肿瘤)诱导的caspases激活坏死因子(TNF) a、Fas和TNF相关的凋亡诱导配体包括GCs的凋亡导致了卵泡的闭锁[22]。但在每个卵泡波中，最终只选择选择一个卵泡成为DF的分子调控机制尚不清楚。
表5 表达下调基因KEGG pathway分析结果
表6 抑制牛卵泡发育的候选基因
图2 牛DF和SF中下调表达基因GO分析
图3 候选表达下调基因在牛DF和SF中的相对表达量 *表示差异显著(P<0.05);**表示差异极显著(P<0.01)
IRELAND等[23-25]把直径≥6 mm的牛卵泡分为E2活跃型卵泡和E2不活跃型泡。E2活跃型卵泡卵泡液中E2>P，且具有健康卵泡的生物学特性；雌激素不活跃卵泡卵泡液中E2
LI等[17]通过对牛PDF1和ODF1 GCs进行转录组测序，发现CPXM1、ARID4B和MAPK13可能在卵泡发育过程中起抑制作用；对猪的闭锁卵泡和健康卵泡GCs测序，发现GADD45A等11个基因可能会抑制卵泡GCs的增殖，进而导致卵泡闭锁[26]。本研究通过对牛DF和SF GCs进行高通量测序，筛选出502个下调表达基因。分别对下调表达基因进行GO、KEGG分析，筛选出4个可能在牛卵泡发育过程中起抑制作用的基因，QRT-PCR验证分析结果显示，EGR1、PPP1R14A和QRFPR在SF和DF中表达趋势与高通量测序结果完全一致。
在小鼠卵巢中的研究发现，老年小鼠卵巢中EGR1的表达显著高于年轻小鼠卵巢，免疫组化染色显示闭锁卵泡内EGR1染色呈强阳性，尤其是凋亡的GCs，同时在小鼠原代GCs中，EGR1的上调抑制了细胞增殖，促进了细胞凋亡[27-28]。在测定真菌毒素DON对牛GCs功能的影响试验中发现， DON通过激活MAPK通路，从而导致转录因子EGR1和EGR3的mRNA丰度随时间和剂量的增加而增加，极大程度地抑制了雌二醇和孕酮的分泌，促进了GCs的凋亡[29]。关于PPP1R14A、QRFPR和OLA1在卵泡中的研究目前尚未见报道，本研究发现PPP1R14A、QRFPR在SF中的表达量显著高于DF，推测其在牛卵泡发育过程中起抑制作用，可能会导致卵泡闭锁。
