基于RNA-Seq技术分析不同产蛋量庄河大骨鸡垂体和下丘脑的可变剪接事件

马智勇，李美成，马巍，王春强

（锦州医科大学畜牧兽医学院，辽宁锦州 121001）

0 引言

RNA-seq技术的产生和发展为生物信息学领域的研究带来了新的契机，其可分别检测有、无参考序列的转录组，进一步发现和寻找新的信息[1]，其产生的大量数据为可变剪接的研究带来了新的思路和方法，被广泛应用于功能基因的挖掘[2]。Gilbert[3]于1978 年首次发现基因的AS现象，是指真核生物转录过程中，基因的一个mRNA前体通过不同的剪接方式切除内含子，重新拼接外显子产生成熟的mRNA 的过程。近年来的试验研究结果陆续表明，可变剪接在许多物种的多个生物学过程中起着重要作用，如遗传信息的稳定性及翻译效率[4]，动植物生长、发育、信号转导以及生物和非生物胁迫的反应等方面[5-8]。

下丘脑-垂体-性腺(HPG)轴控制家禽的生殖活动，是调控家禽生殖发育的主要枢纽[9]。雏鸡性成熟后进入产蛋期，而其性成熟的主要特征是生殖系统的发育成熟，即卵泡发育成熟，正常排卵，形成鸡蛋排出体外。生殖激素在这个过程中发挥着重要作用，而调控生殖激素分泌的中心是HPG轴[10-11]。来自机体内外的影响因素会通过刺激HPO轴或者卵巢、卵泡等靶器官来调控鸡的内分泌激素，从而调控鸡的产蛋性能。

庄河大骨鸡以体大、蛋优、肉美而闻名，是中国辽宁省的地方蛋鸡品种，是中国优良的肉蛋兼用型地方品种鸡[12]。然而，近年来，由于系统选育的缺乏导致其产蛋性能遗传进展缓慢[13-14]，因此，提高产蛋率是繁育工作的重点之一。目前调节母鸡产蛋机制的研究集中在卵巢中不同卵泡发育阶段基因调控的变化[15]，对蛋鸡的垂体和下丘脑的研究相对较少。因此，本研究基于RNA-seq 技术对高产和低产蛋鸡下丘脑和垂体中的AS事件和新转录本的预测进行分析。旨在寻找参与调控母鸡产蛋过程的新基因，为加快庄河大骨鸡产蛋性能的遗传进展提供新思路。

1 材料与方法

1.1 试验动物及分组

试验用6只母鸡均采自于锦州医科大学大骨鸡试验场，屠宰后取垂体和下丘脑组织迅速置于液氮中，并在-80℃下储存备用。根据母鸡在产蛋高峰期39～42周内产蛋的数量，将母鸡分为高产组和低产组（L组3只和N组3只），L组和N组在此期间分别产蛋25个和14个，组内各只母鸡产蛋数量相同。

1.2 RNA分离，文库构建及测序

利用Trizol试剂从垂体和下丘脑组织样品中提取总RNA。Nanodrop 2000 分光光度计（北京科誉兴业科技发展有限公司）检测RNA纯度。通过Qubit 2.0精确定量RNA浓度，并通过Agilent2100生物分析仪（安捷伦科技中国有限公司生产2100 生物分析仪）评估RNA 完整性。 样品经过测试后，利用NEBNext®UltraTMRNA 文库制备试剂盒，试剂盒由安捷伦科技有限公司提供制备文库。于Agilent Bioanalyzer 2100 系统上评估文库质量。文库制备物在Illumina HiSeq 2500 平台上测序，使用PE 150 配对末端测序策略。

1.3 数据分析

Fastq 软件对原始数据进行质控和评估。该研究中鸡参考基因组序列（90 版）由基因组网站下载(ftp://ftp.ensembl.org/pub/current_fasta/gallus_gallus/dna/Gallus_gallus.Gallus_gallus-5.0.dna.toplevel.fa.gz)。Hisat2v2.0.5构建参考基因组索引，将配对末端过滤后读数与参考基因组比对。基于所有数据的基因组定位结果，利用Cufflinks软件组装转录物。

1.4 可变剪接事件分析

使用 rMATS(http://rnaseq- mats.sourceforge.net/index.html)软件比较结果中每个样本的AS 事件进行分类、数量统计以及对差异可变剪接事件进行分析。根据P值和错误发生率(False discovery rate，FDR)鉴定差异AS 事件。以FDR＜0.05 作为差异AS 事件的筛选标准。

1.5 差异剪接基因(DSG)的功能富集分析

应用DAVID 在线功能注释软件(http://david.ncifcrf.gov/summary.jsp)进行差异可变剪接基因的GO注释和KEGG途径的富集分析。

1.6 新转录本的预测和基因结构的优化

使用rMATs 软件分析样品基因剪接与原有剪接模型之间的差异，对映射数据中的读数进行新可变剪接事件预测。使用cufflinks 拼接得到的基因注释文件，与原有基因注释文件进行比较，找出原有注释中未包含的基因并对基因的位置进行优化，补充并修改原有注释文件。

1.7 RT-PCR验证雌激素受体1(ESR1)的AS事件

使用逆转录酶HIScript III RT 试剂盒(Vazyme Biotech Co.，Ltd)、oligo-dT 引物对总RNA 进行逆转录分析。使用DNAman 根据不同转录物的保守区域在单个反应中扩增出4 个剪接变体。用于RT-PCR 的ESR1的引物序列如下：正向 5'-TTGATGGTGCTTTGAG- 3' 和反向5'-GAATGCCAGGTTCTGT-3'，由北京华大基因生物公司完成引物合成。

2 结果

2.1 构建文库的测序

为研究高产和低产蛋鸡垂体和下丘脑中的AS 事件，本试验对庄河大骨鸡的垂体和下丘脑组织进行了RNA-seq 分析。构建了12 个cDNA 文库，共产生约8.25 亿原始读数，经质量控制后共有约7.99 亿干净读数，用于后续分析检测。85.57%～87.92%的序列被定位到参考基因组上。构建的12个文库中，读数映射到“+”链(41.09%～42.23%)和“-”链上(41.15%～42.32%)比例几乎相等（表1），AT和GC曲线均基本重合，表明碱基组成平衡，测序正常。

表1 测序分析

2.2 AS及DSGs的检测

6 个垂体组织文库中，从8720 个基因中检测到16458 AS 事件（平均每个基因可检测到1.89 个AS 事件），差异AS事件总数为769，其中，上调差异AS事件为384，下调差异AS事件为385，可识别5种不同类型的AS事件，按发生频率由高到低分别为：SE(83.87%)、MXE(8.90% )、RI(3.11% )、A3SS(2.41% )、A5SS(1.71%)。基因SLMAP（39AS 事件）、CIB1（28AS 事件）、NRCAM（24AS事件）产生AS事件频率较高，且主要为SE和MXE类型。

6 个下丘脑组织文库中，从9448 个基因中检测到19021 AS 事件（平均每个基因可检测到2.01 个AS 事件），差异AS事件总数为813，其中，上调差异AS事件为418，下调差异AS 事件为395，可识别5 种不同类型的AS事件，按发生频率由高到低分别为：SE(83.77%)、MXE(10.06% )、RI(2.64% )、A3SS(2.07% )、A5SS(1.46%)。基因NRCAM（38AS 事件）、SLMAP（33AS 事件）、CIB1（28AS事件）产生AS事件频率较高，主要类型同垂体文库检测结果（表2）。综上，12 个组织文库中，SE类型发生的比例均为最高，A5SS发生的比例均为最低；基因NRCAM、SLMAP、CIB1在两类组织文库中均发生较高频率的AS事件。

表2 AS分类、数量及DSGs

以FDR＜0.05为标准对rMATS软件分析结果进行差异筛选，垂体组织文库中共筛选到644 个显著的差异剪接基因，在SE、MXE、RI、A3SS、A5SS事件中鉴定到的差异剪接基因数分别为548、159、32、18、11；下丘脑组织文库中共筛选到680 个显著的差异剪接基因，在SE、MXE、RI、A3SS、A5SS 事件中鉴定到的差异剪接基因数分别为613、142、31、24、9；两类文库差异剪接基因中，SE 占比例均为最高，A5SS 占比例均为最低，趋势同AS事件相似（表2、表3）。

表3 参与AS事件的基因（5例）

2.3 基因的功能注释和富集分析

为进一步阐明DSGs的功能，本试验进行了GO和KEGG 途径富集分析，以注释DSGs 并研究其分布。垂体组织文库中，644 个DSGs 中的329 个基因得到功能注释，基因富集在生物过程、细胞组分、分子功能的占比分别为27%、33%、40%（图1A）。68%的DSG 同时富集于3 个功能分类。下丘脑组织文库中，680 个DSGs中的337个基因得到功能注释，基因富集在生物过程、细胞组分、分子功能的占比分别为26%、52%、22%（图1B）。64%的DSG同时富集于3个功能分类。

图1 GO分类

垂体组织文库中，24 个DSGs 富集于3 个KEGG途径，其中MAPK信号传导途径显著富集(P＜0.05)，12个DSGs 参与MAPK 途径。下丘脑组织文库中，22 个DSGs富集于3个KEGG途径，粘附连接途径显著富集(P＜0.05)，7个DSGs参与该途径（表4）。

表4 KEGG途径富集分析以注释DSGs

2.4 新转录本预测和基因结构优化

构建的12个文库中共获得6758个新转录本，长度分布于170～35112 bp之间。根据外显子的预测结果，新转录本包含1～68个外显子。部分新转录本的结构和功能预测见表5。构建的文库中共12486个基因结构得到优化，6230个基因的5’端发生延伸，6256个基因的3’端发生延伸，部分已注释基因的优化结果见表6。

表5 新转录本结构注释结果（7例）

表6 已知基因结构优化（10例）

2.5 ESR1基因的RT-PCR检测

分析中发现雌激素受体基因ESR1具有较多的AS事件，产生4种类型AS（图2）。RT-PCR结果发现存在3个不同长度的片段（分别为845、624、509 bp）（图3），表明垂体中的ESR1确实存在于不同的AS事件中。

图2 4种AS类型

图3 ESR1基因的RT-PCR检测

3 讨论

本研究基于RNA-seq 技术，以鸡基因组序列（版本90）作为参考基因组，构建了鸡垂体和下丘脑文库。检测到5 种类型的AS 事件(SE、RI、A5SS、A3SS、MXE)，均为普遍认同的AS的5种基本形式。AS事件的类型在不同生物中发生情况各不相同：据报道，在哺乳动物特别是小鼠中，SE 是AS 事件的最常见类型[8,16]；但在植物研究中发现，RI 是水稻和竹子的主要AS 事件类型，分别占总数的47%和27.43%[13,17]；在对低等生物的研究中，已发现RI 是根瘤菌的主要AS 事件类型，占总数的84.3%[14]。本实验对蛋鸡的垂体和下丘脑的研究表明，SE占AS事件的比例最高，分别为83.87%和83.77%，A5SS 的比例最低，分别占AS 事件的1.71%和1.46%，与大多数哺乳动物相同。

家禽可以通过HPG 轴合成和释放激素以调节繁殖活动。家禽中的HPG 轴与哺乳动物中的HPG 轴差异较大，表现出可塑性[18]。本研究中，低产蛋鸡下丘脑中的22 个DSGs 被定位到粘连、紧密连接和钙信号传导途径，26 个DSGs 参与了Ca2+结合或释放通道活性的生物学过程。目前，普遍认为钙信号传导是促性腺激素释放激素(GnRH)调节促性腺激素合成和释放的主要机制。研究表明，下丘脑钙离子信号传导途径异常会影响中枢神经系统的发育，调节下丘脑进食神经元的兴奋性，从而影响动物的进食以减少能量的摄入，进而影响营养物质的吸收和运输[15]。同时，大量研究证实，钙离子在调节GnRH 快速分泌促性腺激素中起重要作用[19-20]。因此，我们推测下丘脑Ca2+相关通路异常可能降低底物合成效率，影响促性腺激素的合成和分泌。

本研究中，低产蛋鸡垂体中的24 个DSGs 定位于心肌细胞中的MAPK信号通路、心肌收缩或肾上腺素能信号通路。研究表明，GnRH与GnRHR的结合可激活MAPK途径，但机制尚不完全清楚，GnRHRs可能参与了许多蛋白质信号复合物的组装[14]。在小鼠垂体细胞中，GnRHRs与ERK2共定位于脂筏[21]，ERK 的GnRH 活化也参与钙离子募集过程，这可能部分解释了GnRH在垂体中的作用机理。

下丘脑中合成的GnRH通过作用于垂体前叶以将其天然的高亲和力跨膜受体GnRHR结合在细胞表面上，从而刺激促性腺激素的信号级联[22]。迄今为止，已经使用常规的生化和生物信息学工具鉴定了许多GnRH 和GnRHR基因[23-24]。已在哺乳动物中鉴定出GnRHR家族的3个成员[25]，但尚未对家禽中的GnRHR基因进行系统的研究。本研究中确定了垂体中5 个GnRHR的AS事件，其中4个为SE、1个为MXE。

在GO分析中，垂体的4个DSGs，即神经生长因子(NGF)，视黄酸8(STRA8)，EPH 受体A5(EPHA5)和mutS 同系物4(MSH4)注释为卵泡发育的生物学过程。神经生长因子(NGF)AS 事件在低产和高产母鸡垂体中有所不同。虽然关于NGF对鸡垂体生殖功能的影响少有报道。但近年来，NGF在哺乳动物的非神经系统中的作用已成为研究热点。如，已经在NGF敲除小鼠模型中系统地研究了NGF在卵泡形成和早期卵泡发育中的作用[26]。缺乏NGF基因的小鼠卵巢表现出严重的异型增生，卵母细胞不能被颗粒细胞包围而形成卵泡结构。初级卵泡和次级卵泡的数量也显著减少，表明NGF对于次级卵泡的形成至关重要[27]。NGF可以影响卵巢中某些基因的表达。NGF的丢失会导致FSHR基因在卵巢中的表达显著降低。相反，当使用外源性NGF治疗新生野生型小鼠卵巢时，FSHR的表达显著增加[28]。需要进一步的研究来确定母鸡垂体中NGF的不同AS事件是否以内分泌方式影响母鸡卵泡的发育，并进一步影响产蛋过程。STRA8编码视黄酸反应蛋白。小鼠中STRA8的同源蛋白已被证明参与精子发生和卵子发生中减数分裂起始的调控[29]。但小鼠和母鸡STRA8蛋白之间的特异性差异表明，这些同源物在功能上并不相同。编码STRA8的基因被认为在母鸡产蛋过程中起作用。但有研究表明，EPHA5和MSH4参与卵巢卵泡颗粒细胞和卵泡内膜细胞的细胞增殖和分化以及激素的合成和分泌[27-30]。EPHA5和MSH4在垂体中的作用和调节机制尚未见报道。不同AS形式的具体功能需要进一步研究。

在垂体中发现了许多AS 事件，如据报道雌激素受体1(ESR1)在垂体中发挥作用。雌激素对于动物卵泡的发育至关重要，并且2 种雌激素受体具有独立的功能[31]。雌激素受体(ERs)是配体激活的转录因子，当结合在细胞质中时，会转移到细胞核，引起以ER为中心的信号传导过程[32]。雌激素受体的2 种亚型ESR1和ESR2存在于垂体中，但受体敲除模型表明ESR1是雌激素负反馈调节的主要受体[33]。此外，雌激素可负面调节女性的垂体功能并调节其生殖活动[34]。此外，有研究发现ESR1基因多态性与中国大骨鸡的产蛋特性密切相关[35]。解剖母鸡后，发现蛋鸡卵巢上的大卵泡数量明显高于高产蛋鸡组。由于ESR 表达是在雌激素刺激下发生的，因此我们推测排卵紊乱会导致大卵泡的持续存在，ESR1的不同亚型可能通过介导的负反馈调节垂体功能，进而影响母鸡的产蛋性能。

4 结论

AS 事件类型在庄河大骨鸡下丘脑和垂体中分布与大多数哺乳动物相同；低产蛋鸡下丘脑和垂体中的DSGs 主要定位于钙信号途径和MAPK 途径，参与促性腺激素的合成和分泌，但机制尚不清楚；垂体中的ESR1存在于不同的AS事件类型中。综上，AS事件对母鸡的产蛋数量起着至关重要的作用，为蛋鸡产蛋性能的研究奠定一定的基础。