基于RNA-Seq 高通量测序技术对肉牛下丘脑基因表达分析及新基因挖掘研究

雷 虹，孙 鏖，胡雄贵，李付强，张 爽，易康乐*

(1.湖南省畜牧兽医研究所，湖南 长沙 410130； 2.湖南天华实业有限公司，湖南 涟源 417100； 3.吉林省永吉县口前镇畜牧兽医站，吉林 永吉 132100)

在家畜动物发情周期和生殖活动中，下丘脑-垂体-卵巢(the hypothalamic-pituitary-ovary，HPO)轴起着十分重要的作用[1-3]。经典的繁殖免疫学理论中，雌性动物的繁殖活动过程要受到下丘脑、垂体前叶和性腺之间复杂的相互作用和调控。一般认为，下丘脑部分泌促性腺激素释放激素(GnRH)并形成脉冲，驱动了整个性腺轴由大脑到卵巢的信号调节通路。 垂体和卵巢组织中的配子发生和类固醇生成都要受到GnRH的调节。 垂体促性腺激素上的受体与GnRH 具有极高结合亲和力，从而刺激促性腺激素脉冲式的释放。在GnRH 的影响下，卵泡刺激素(FSH)和黄体生成激素(LH)对卵泡的发育和卵母细胞的成熟进行调控。与此同时，脉动式的GnRH 释放方式对于动物的正常生殖功能起到了至关重要的作用[4-5]。因此，鉴于GnRH 在雌性动物繁殖活动中至关重要的作用，体现出下丘脑部在动物的发情调控及生殖活动中起到了关键主导作用。没有下丘脑分泌的GnRH，促性腺和性腺就不起作用[6-7]。下丘脑的激素分泌活动受多种因素影响，例如品种、营养和季节，而这些因素是如何影响下丘脑生理活动的分子机制仍不清楚。

在动物研究领域中，基于高通量测序 RNA-Seq 技术的转录组学的运用越来越深入[8-9]。在肉牛遗传标记基因功能注解数据库中，大部分的基因功能注释是通过人类、小鼠和其他哺乳动物基因库中比对后预测得到的,而肉牛特有特异基因的挖掘和注释可能还不全面或存在遗漏。因此，本研究拟对不同年龄阶段的安格斯母牛下丘脑组织进行采集后，进行mRNA水平转录组分析。通过应用RNA-Seq 技术对其进行高通量转录组测序分析，筛查相关表达基因，筛选出与性腺轴调控繁殖过程相关的基因与信号通路，为进一步了解肉牛的环境适应性、分子育种及探析繁殖性能等诸多方面提供依据。

1 材料与方法

1.1 样品采集在湖南涟源国家安格斯牛核心育种场选取9头半同胞(同父异母)、体型发育基本一致、健康的9头安格斯母牛，6，18和30月龄各3头。于2018年4月9日进行屠宰。肉牛屠宰后立即采集卵巢组织样本,置于液氮中速冻,并于-70℃保存备用,用于提取组织总RNA。本研究不涉及动物伦理及福利问题,所有样本均采集自政府审查批准的相关屠宰场。

1.2 RNA定量和鉴定测试RNA样品的纯度、浓度和完整性，以确保将合格的样品用于转录组测序。每个样品总共1μg RNA用作RNA样品制备的输入材料。使用Illumina的NEBNext UltraTMRNA库制备试剂盒(美国NEB)生成测序库，并将索引代码添加到每个样品的属性序列中。对文库质量进行评估。

1.3 聚类和排序使用BMKCloud(www.biocloud.net)进行RNA-Seq分析。根据制造商的说明，使用TruSeq PE Cluster Kit v4-cBot-HS(Illumia)对样品数据在cBot簇生成系统上进行聚类分析，并将文库的制备物在Illumina系统平台上测序，从而生成配对末端的读取数据。

1.4 质量控制Fastq格式的原始数据(原始读取)首先通过内部perl脚本进行处理。在此步骤中，通过从原始数据中删除包含适配器的读取，包含ploy-N的读取和低质量的读取来获得clean data。同时，计算了clean data的Q20、Q30、GC含量和序列重复水平。所有下游分析均基于高质量的clean data。

1.5 比较分析低质量序列读数将从数据集中删除。原始序列在数据处理后被转换为clean reads。然后将这些clean reads进行参考基因组定位，进一步分析和注释仅具有完全匹配或一个错配的读数。使用Hisat2工具软件绘制参考基因组图[10-11]。

1.6 基因表达水平的定量基因表达水平的定量通过每百万个作图片段转录本的每千碱基片段来估算基因表达水平。

1.7 基因功能注释用于基因功能的注释源于以下几类数据库：蛋白质家族数据库(Pfam)和蛋白质直系同源簇数据库(KOG/COG)KEGG；Ortholog数据库(KO)和 Gene Ontology 数据库(GO)；人工注释和审查的蛋白质序列数据库(Swiss-Prot)； NCBI非冗余蛋白质序列数据库(Nr)和非冗余核苷酸序列数据库(Nt)[12-15]。

2 结果

2.1 总RNA提取及文库构建对9头牛的下丘脑组织进行总RNA提取。提取后的mRNA用核酸蛋白测定仪进行测定，在总RNA 纯度(D260/D280)大于2.0以上、质量浓度大于10 mg/L的情况下，进行下一步试验。应用qPCR方法对文库质量进行定量验证。

2.2 转录组测序质量分析对牛下丘脑的9个样品的转录组测定数据分析，共获得59.74 Gb clean data，各组的样品clean data均达到6.64 Gb，Q30碱基百分比平均在94.60%及以上(表1)。分别将各样品的clean reads与指定参考基因组数据库进行比对，比对率从94.48%～94.68% 不等(表2)。统计了各样品mapped reads在指定的参考基因组外显子、内含子和基因间区的数目，对各样品基因组不同区域上mapped reads的分布情况进行绘图(图1)。

表1 测序数据统计表

表2 样品测序数据与所选参考基因组的序列比对结果统计表

图1 不同样本测序序列与参考基因组序列各位置比对百分比

2.3 表达基因分析

2.3.1表达基因筛选 对于9个样品数据，采用并集的方法，对表达阈值≥0.1 的基因进行筛选，并将相同的样品进行混池后进行比较后，维恩图绘制(图2)。结果显示，6月龄组共有21 089个基因表达，18月龄组共有20 627个基因表达，30月龄共有20 782个基因表达。有18 979个基因在所有月龄段中都有表达，有1 706个基因只在2个月龄段中共同表达，有2 149个基因只在其中的1个月龄段中表达。

图2 不同组织表达基因维恩图

2.3.2表达基因细胞注释和富集分析 对表达基因进行数据库的功能注释结果显示，下丘脑组织共有2 621个基因在GO数据库中得到注释(图3)，13 453个基因在KEGG数据库中得到注释，12 930个基因在KOG数据库中得到注释(图4)。

图3 表达基因GO富集柱状图

图4 表达基因KOG富集柱状图

2.3.3表达基因的细胞功能与代谢通路分析 对所有筛选出表达基因进行GO数据库的标注分析，富集较多的为细胞过程、细胞组成和附着黏连。通过GO数据库二级分类注解分析表明，共发现12类与器官发育、性激素反应、受精活动和胚胎发育相关。其中新基因主要集中在动物器官发育、单个组织繁殖过程和授精等细胞和生物过程中(表3)。

表3 卵巢发育、受精活动和胚胎发育相关表达 基因GO分类表

将所有表达基因标注到KOG数据库中，富集较多的为普通功能预测(general function prediction only)、信号传导机制(signal transduction mechanisms)和翻译后修饰、蛋白转换、分子伴侣(posttranslational modification,protein turnover,chaperones)(表4、图4)。

表4 相关表达基因KOG分类表

将所有表达基因标注到KEGG数据库中。下丘脑组织共富集到300个信号通路。对富集最显著的前20条信号通路进行绘图(图5)。在所有的信号通路中，筛选与繁殖相关的信号通路进行归类整理。结果显示，在cAMP信号通路、Jak-STAT信号通路、卵巢类固醇激素信号通路等16个信号通路中有参与表达(表5)。

图5 所有表达基因KEGG富集散点图(前20)

表5 与繁殖相关KEGG信号通路统计表

3 讨论

在鸡、番鸭、绵羊等物种中均有利用转录组测序技术对下丘脑组织进行测序的相关研究报道[16-19]。但以不同年龄阶段母牛为研究对象，进行系统得下丘脑转录组测序的相关报道不多。本试验以安格斯牛为研究对象，通过选取初情期左右的安格斯牛母牛下丘脑组织进行RNA-Seq和生物信息学分析，希望深度挖掘安格斯牛与生长发育、生殖激素调控和繁殖性能相关的分子标记信息。

本研究对9头半同胞母牛的下丘脑组织，按不同年龄段分组后进行混池建库测序。测序结果显示，6月龄组共有21 089个基因表达，18月龄组共有20 627个基因表达，30月龄组共有20 782个基因表达。有18 979个基因在所有组织中都有表达，分别占到6，18和30月龄基因表达总数的89.99%，92.01%和91.32%。结果显示，6月龄下丘脑的基因表达数量要多于其他月龄段，说明肉牛在生长发育的初期，下丘脑的基因调控活动较其他时期更为活跃。

为了解各月龄段下丘脑组织表达基因相关功能，对测序筛选出的表达基因进行了KOG、GO和KEGG数据库比对注释分析。

通过KOG数据库注释分析发现，有14 459个基因被注释。其中肉牛各月龄段的下丘脑组织中富集最多的是普通功能预测，为3 029个，占到全部注解数的20.95%。

通过GO数据库注释分析发现，在6月龄组中富集到的GO条目为51个，18月龄组中富集到的GO条目为52个，30月龄组中富集到的GO条目为53个。在GO的分子功能分析中，3个年龄段富集最多的都是细胞过程、细胞组成和分子结合绑定。在对表达基因进行GO数据库二级分类注解分析后，发现在所有涉及器官发育、生殖调控相关的调控分类中，很大比例的基因(拷贝)涉及到动物器官发育过程，而且这些基因大部分都是新拷贝基因。说明下丘脑组织与肉牛生长发育过程密切相关，且有大量新基因有待深入研究与探讨。

对所有筛查出的表达基因在进行了GO数据库细胞功能注释分析后又进行了KEGG信号通路分析。其中富集最多的是癌症通路。出现这一结果的原因可能有在所有基因组功能注释的功能基因中，与癌症和疾病相关的基因数量和信号通路相对较多。其次参与到癌症通路的大部分基因可能同时涉及和参与了细胞的分化、增殖和器官组织的发育。但癌症通路在性腺轴生理活动中确切功能及作用尚不明确，有待进一步的研究和验证[20]。

与此同时发现，在3个月龄段的下丘脑组织都有表达的并与动物繁殖性状相关信号通路有17条。其中GnRH 信号通路是性腺轴主导的调控动物生殖活动的经典传导通路[21-22]。同时本研究注意到，Notch信号通路、神经营养因子信号通路和血管内皮生长因子信号通路的参与基因在性腺轴中检出频率最高，说明这些信号通路可能也通过性腺轴对雌性动物生殖活动产生影响。

研究表明，血管内皮生长因子(VEGF)家族及其受体在血管生成过程中起到重要的调控作用。VEGF家族共有7个因子，分别是VEGF-A、VE-GF-B、VEGF-C、VEGF-D、VEGF- E和胎盘生长因子(PIGF)-l 、PIGF-2，同时还有2个受体分别为KDR和FLT1[23]。研究已证实，VEGF家族成员在雌性哺乳动物的卵泡发育、卵母细胞的成熟、黄体的形成与功能的维持、子宫的周期性变化等生理过程中都起着重要的作用[24-26]。Notch 信号通路参与了生物的很多生理活动，并通过与细胞间的互作方式来调控机体各组织器官的分化和发育[27]。现已发现Notch 信号通路的主要参与因子由4种受体因子和5 种配体共同组成[28]。在成体干细胞的不对称分裂过程，与细胞凋亡和衰老过程中，Notch 信号通路起到了关键作用[29]。其中在卵巢生殖干细胞增殖分化及卵巢衰老过程中Notch 信号通路所起到的调控作用已被证实[30]。近年来的研究结果证实，Hippo 信号通路在机体组织和器官发育过程中发挥重要的调控作用[31-32]。现有的研究表明Hippo 信号通路参与到了动物生殖生理活动过程中，在卵泡的调控生长、原始卵泡的启动、卵巢生殖干细胞功能调节和卵巢组织延缓衰老等过程中起到了重要的调控作用[33-35]。

本研究通过对不同月龄的肉牛下丘脑组织进行转录组测序，对数据进行比对、聚类等分析后得出以下结论：血管内皮生长因子信号通路、Notch信号通路和Hippo 信号通路可能通过下丘脑和HPO轴参与了肉牛繁殖相关调控活动。