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曲线均基本重合,表明碱基组成平衡,测序正常。
样品名称 | 原始读数 | 干净读数 | 总定位基因数目 | 外显子读数 | 映射正链读数 | 映射负链度数 |
LP_1 | 68317482 | 66328810 | 57680205(86.96%) | 66.60% | 27816090(41.94%) | 27768078(41.86%) |
LP_2 | 72359294 | 69806998 | 60386063(86.5%) | 62.10% | 29087549(41.67%) | 29092824(41.68%) |
LP_3 | 63889132 | 62116040 | 5366 920(86.4%) | 64.40% | 25918157(41.73%) | 25814561(41.56%) |
NP_1 | 64341034 | 61029484 | 52964336(86.78%) | 64.10% | 25448135(41.7%) | 25489645(41.77%) |
NP_2 | 69054884 | 67531462 | 59190144(87.65%) | 65.20% | 28511031(42.22%) | 28463500(42.15%) |
NP_3 | 62110264 | 60924682 | 52608969(86.35%) | 64.60% | 25342586(41.6%) | 25363479(41.63%) |
LH_1 | 79744218 | 77900 688 | 68487901(87.92%) | 69.50% | 32899952(42.23%) | 32968715(42.32%) |
LH_2 | 65534460 | 62447384 | 54556009(87.36%) | 66.60% | 26194279(41.95%) | 26265470(42.06%) |
LH_3 | 79835258 | 78021862 | 67701041(86.77%) | 68.20% | 32517332(41.68%) | 32596024(41.78%) |
NH_1 | 76681332 | 74134976 | 64570880(87.1%) | 70.70% | 30962806(41.77%) | 31058491(41.89%) |
NH_2 | 49332078 | 46982632 | 41205550(87.7%) | 67.10% | 19804054(42.15%) | 19834071(42.22%) |
NH_3 | 73617140 | 71760796 | 61408353(85.57%) | 66.90% | 29489160(41.09%) | 29529056(41.15%) |
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事件。
基因文库 | 基因 | AS事件类型 | AS事件总数 | 差异事件总数 | 差异剪接基因 |
垂体 | 8720 | SE | 13804 | 549(264:285) | 548 |
MXE | 1464 | 159(91:68) | 159 |
RI | 512 | 32(15:17) | 32 |
A3SS | 397 | 18(7:11) | 18 |
A5SS | 281 | 11(7:4) | 11 |
下丘脑 | 9448 | SE | 15933 | 610(312:298) | 613 |
MXE | 1913 | 142(71:71) | 142 |
RI | 504 | 31(19:12) | 31 |
A3SS | 394 | 24(12:12) | 24 |
A5SS | 277 | 6(4:2) | 6 |
以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)。
基因文库 | 基因 | SE | MXE | RI | A3SS | A5SS |
垂体 | SLMAP | 25 | 14 | 0 | 0 | 0 |
CIB1 | 19 | 9 | 0 | 0 | 0 |
NRCAM | 20 | 4 | 0 | 0 | 0 |
CADM1 | 20 | 12 | 0 | 0 | 0 |
CARMIL3 | 18 | 3 | 1 | 0 | 0 |
下丘脑 | SLMAP | 23 | 10 | 0 | 0 | 0 |
CIB1 | 19 | 9 | 0 | 0 | 1 |
NRCAM | 25 | 13 | 0 | 0 | 0 |
EPB4IL3 | 22 | 6 | 0 | 1 | 0 |
ABI3BP | 19 | 4 | 0 | 0 | 0 |
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个功能分类。
垂体组织文库中,24个DSGs富集于3个KEGG途径,其中MAPK信号传导途径显著富集(P<0.05),12个DSGs参与MAPK途径。下丘脑组织文库中,22个DSGs富集于3个KEGG途径,粘附连接途径显著富集(P<0.05),7个DSGs参与该途径(表4)。
基因文库 | 途径 | 描述 | P值 | DSGs |
垂体 | gga04010 | MAPK signaling pathway | 0.029472063 | 12 |
gga04260 | Cardiac muscle contraction | 0.081935761 | 5 |
gga04261 | Adrenergic signaling in cardiomyocytes | 0.086900369 | 7 |
下丘脑 | gga04520 | Adherens junction | 0.011715804 | 7 |
gga04530 | Tight junction | 0.062176625 | 6 |
gga04020 | Calcium signaling pathway | 0.083233215 | 9 |
2.4 新转录本预测和基因结构优化
构建的12个文库中共获得6758个新转录本,长度分布于170~35112 bp之间。根据外显子的预测结果,新转录本包含1~68个外显子。部分新转录本的结构和功能预测见表5。构建的文库中共12486个基因结构得到优化,6230个基因的5’端发生延伸,6256个基因的3’端发生延伸,部分已注释基因的优化结果见表6。
染色体编号 | 来源 | 长度 | 正负链 | 外显子数 | 功能预测 |
4 | Novel Gene | 21214 | + | 1 | nucleic acid binding |
22 | Novel Gene | 11154 | - | 4 | histone acetyltransferase KAT2B isoform X2 |
2 | Novel Gene | 8241 | + | 6 | heat shock transcription factor, X-linked-like |
12 | Novel Gene | 1896 | - | 15 | hypothetical protein ENH_00063870 |
6 | Novel Gene | 34081 | + | 28 | purine ribonucleoside binding |
15 | Novel Gene | 7214 | - | 41 | nucleoside binding |
17 | Novel Gene | 5900 | + | 68 | extracellular matrix structural constituent |
基因ID | 染色体编号 | 正负链 | 原始跨度 | 延伸跨度 |
ENSGALG00000000003 | 1 | + | 20174067-20177667 | 20157937-20179163 |
ENSGALG00000000081 | 16 | - | 205312-210738 | 204525-211382 |
ENSGALG00000000129 | 11 | + | 19064511-19086448 | 19064481-19107184 |
ENSGALG00000000233 | 27 | - | 4203293-4208774 | 4201773-4212041 |
ENSGALG00000000309 | 26 | + | 789839-794797 | 787968-795597 |
ENSGALG00000000296 | 5 | - | 27201732-27283175 | 27195404-27283845 |
ENSGALG00000000853 | 17 | + | 10878731-10887869 | 10872754-10887869 |
ENSGALG00000000862 | 23 | - | 1766445-1795825 | 1765749-1795881 |
ENSGALG00000000978 | 21 | + | 905170-916809 | 898302-932294 |
ENSGALG00000001042 | 19 | - | 615601-657520 | 613509-657908 |
2.5 ESR1基因的RT-PCR检测
分析中发现雌激素受体基因ESR1具有较多的AS事件,产生4种类型AS(图2)。RT-PCR结果发现存在3个不同长度的片段(分别为845、624、509 bp)(图3),表明垂体中的ESR1确实存在于不同的AS事件中。
图3 ESR1基因的RT-PCR检测注:M:DL2000 DNA Marker;1:垂体;2:对照 |
Full size|PPT slide
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参与了Ca
2+结合或释放通道活性的生物学过程。目前,普遍认为钙信号传导是促性腺激素释放激素(GnRH)调节促性腺激素合成和释放的主要机制。研究表明,下丘脑钙离子信号传导途径异常会影响中枢神经系统的发育,调节下丘脑进食神经元的兴奋性,从而影响动物的进食以减少能量的摄入,进而影响营养物质的吸收和运输
[15]。同时,大量研究证实,钙离子在调节GnRH快速分泌促性腺激素中起重要作用
[19-20]。因此,我们推测下丘脑Ca
2+相关通路异常可能降低底物合成效率,影响促性腺激素的合成和分泌。
本研究中,低产蛋鸡垂体中的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事件对母鸡的产蛋数量起着至关重要的作用,为蛋鸡产蛋性能的研究奠定一定的基础。
{{custom_sec.title}}
{{custom_sec.title}}
{{custom_sec.content}}
[1] Cutter S L,Vulnerability to Environmental Hazards[J].Progress in Human Geography,1996, 20: 529-539.
[2] Mejia-Navarro M. Geological Hazards, Vulnerability, and Risk Assessment Using GIS: Model for Glenwood Springs, Colorado[J]. Geomorphology, 1994, 10 (1-4): 331-354.
[3] Blaikei P, Cannon T, Davis I, et al. At Risk: Natural hazard, people vulnerability and Disasters[M]. London: Routledge, 1994: 147-167.
[4] 黄崇福,自然灾害风险评价:理论与实践[M],北京:科学出版社,2005: 96-98.
[5] 史培军,论灾害研究的理论与实践[J].南京大学学报(自然科学版),1991,11:37-42.
[6] 史培军,再论灾害研究的理论与实践[J].自然灾害学报,1996, 5(4): 6-17.
[7] 陈寿军,沿海地区雷电活动特征及风险区划研究[J].南京师范大学硕士论文, 2013.
[8] 蒋勇军,况明生,匡鸿海,等.区域易损性分析、评估及易损度区划-以重庆市为例[J].灾害学, 2001, 16(3): 59-64.
[9] 刘岩,李征,程向阳等.安徽省雷电灾害风险区划[J].南京信息工程大学学报,2014, 6(2): 163-168.
[10] 尹娜,肖稳安,雷灾易损性分析、评估及易损度区划[J].热带气象学报,2005, 21(4): 441-449.
[11] 严春银,江西省雷电灾害易损性分析及其区划[J].江西科学,2006, 24(2): 131-135.
[12] 袁湘玲,纪华,程琳.基于层次分析模型的黑龙江省雷电灾害风险区划[J].暴雨灾害,2010, 29(3): 279-283.
[13] 程萌.菏泽市雷暴的形成及其天气形势分析[J].现代农业科技,2016, 4: 232-234.
[14] 中国气象局雷电防护管理办公室,全国雷电灾害汇编(内部资料)[M].2004~2013.
[15] 菏泽市统计局.菏泽统计年鉴2013[M].北京:中国统计出版社,2013.
[16] 王以彭, 李结松, 刘立元. 层次分析法在确定评价指标权重系数中的应用[J]. 第一军医大学学报, 1999, 19(4): 377–379.
[17] Saaty T. The Analytical Hierarchy Process[M]. New York: Mcgraw-Hill, 1980.
[18] 唐巧玲,山东省雷电活动特征及雷电灾害风险区划研究[J].兰州大学硕士论文,2013.
[19] 张会,张继权,韩俊山.基于GIS技术的洪涝灾害风险评估与区划研究[J].自然灾害学报,2005, 14(6): 141-146.
[20] 张京红,田光辉,蔡大鑫,等.基于GIS技术的海南岛暴雨洪涝灾害风险区划[J].热带作物学报,2010, 31(6): 1014-1019.
[21] 李楠,任颖,顾伟宗,等.基于GIS的山东省暴雨洪涝灾害风险区划[J].中国农学通报,2010, 26(20): 313-317.