肉鸡肌肉与脂肪组织基因组差异剪接基因分析

2018-10-30孙艳发钟礼鹏许卫华李虹仪黄其春

畜牧兽医学报 2018年10期

张敏，王杰，孙艳发*，钟礼鹏，李焰，许卫华，李虹仪，张茂，黄其春

(1. 龙岩学院生命科学学院，龙岩 364000； 2. 中国农业科学院北京畜牧兽医研究所，北京 100193； 3. 福建省家畜传染病防治与生物技术重点实验室，龙岩 364000)

可变剪接(alternative splicing，AS)是mRNA通过不同的方式产生多种不同剪接RNA转录产物的过程[1]，是高等真核生物蛋白质多样性的重要基因调解机制之一[2]。人基因组中有约95%的多外显子基因通过100 000多种不同的AS事件产生转录产物[3]。新一代RNA测序(RNA-seq)技术已经证实RNA转录产物的可变剪接是人基因组中编码蛋白质多样性的重要机制[4]。通过对基因组中AS事件的鉴定，已经发现一些基因的异常AS可导致人类多种疾病[5-7]。

畜禽基因和基因组中一些AS事件逐渐被发现和证实。研究人员通过设计跨越Lmbr1基因外显子4的引物，发现与丝羽乌骨鸡趾型呈现显著相关的Lmbr1基因中Lmbr1-α为其主要转录产物，而Lmbr1-β为次要转录产物[8]。经重叠延伸RT-PCR结合巢式PCR法扩增，在东北民猪中发现TLR4基因的3个可变剪接体[9]，鸭TLR4基因中也发现了两种不同的剪接体[10]。研究人员通过PCR扩增和5′RACE 方法获得白色和黑色山羊皮肤组织Agouti基因部分基因组序列及mRNA序列，经过拼接比对，在白色山羊中发现5种Agouti基因可变剪接体，而黑色山羊中发现2种可变剪接体[11]。运用RNA-seq技术，畜禽基因组的AS事件得到了全面的鉴定和分析。通过对绒山羊肌内脂肪细胞在成熟前后RNA-seq数据分析发现，绒山羊肌内脂肪细胞中存在A3SS、A5SS、SE和RI 4种主要可变剪接类型，其中A3SS最为普遍[12]。牛基因组的21 755个基因中4 567(21%)个基因发生了可变剪接[13]。猪睾丸基因组中存在92 738个AS事件，以TSS、TTS、SE和AE 4种类型为主；通过分析发现，与睾丸素代谢相关的基因中63个基因发生了可变剪接，以TSS和TTS为主[14]。这些AS事件的鉴定与分析有利于全面研究基因组中AS事件和蛋白质的多样性。由于AS事件大多数表现出组织特异性变异，其中10%～30%的AS事件易受细胞、组织和条件特异性调节的影响[3]。上述的研究并未对组织间差异剪接，特别是差异剪接基因(differential splicing gene，DSG)进行深入的分析。差异剪接分析描述了两个样本间AS位点使用情况和剪接异构体的差异，这对基因组中AS及其调控机制的研究至关重要，同时揭示了差异基因表达分析错过的功能多样性[15]。DSG指不同组织样品间因差异剪接而导致基因亚型的相对丰度有差异的基因。因此，DSG分析对于研究组织间AS及其调控机制具有重要的意义。

本研究以AA肉鸡RNA-seq测序数据为基础，通过分析鉴定肉鸡肌肉和脂肪组织基因组中DSG，分析这些DSG功能和其信号通路，为进一步解析肉鸡基因组中组织间AS事件及组织间基因表达的特异性奠定理论基础。

1 材料与方法

1.1 试验样品

AA肉鸡母鸡苗由福建正大食品股份有限公司提供，于龙岩学院动物饲养房饲养至42日龄。肉鸡饲养管理和屠宰参考肉鸡生产技术规范(DB11/T 328-2005)与肉鸡屠宰质量管理规范(NY/T 1174-2006)。3只肉鸡禁食12 h后进行屠宰，取胸肌和腹脂组织样品转入2 mL RNase-free离心管中，标记后迅速置于液氮中保存备用。

1.2 组织总RNA提取

采用TRIzol试剂(Invitrogen，Carlsbad，美国)提取组织样品总RNA，使用0.8%的琼脂糖凝胶电泳检测总RNA是否降解和受污染，并通过Agilent 2100 生物分析仪使用RNA 6000 Nano试剂盒(Agilent Technologies，美国)检测总RNA的浓度、完整度(RIN)、28S/18S比值和片段大小。用于构建RNA-seq文库和测序的总RNA均满足以下条件：总RNA量≥6 μg，RIN≥8.0，28S/18S≥1.0。

1.3 肉鸡胸肌和腹部脂肪组织RNA-seq文库的构建和测序

RNA-seq文库构建和测序均由深圳华大基因科技有限公司完成。分别取3只肉鸡的胸肌和腹脂组织RNA样品各3份，构建6个文库，使用Illumina Hiseq 2500测序平台进行测序，测序获得长度为150 bp的序列(reads)。去除带接头的reads、未知碱基N含量大于5%的reads、去除低质量的reads(质量值低于 15 的碱基占该reads总碱基数的比例大于 20% 的reads为低质量的reads)，得到clean reads。

1.4 肉鸡胸肌和腹部脂肪组织参考基因组比对及AS事件和DSG检测

得到clean reads之后，使用HISAT软件[16]将clean reads比对到鸡参考基因组序列。使用StringTie软件[17]对每个样品进行转录本重构，然后用Cuffmerge软件[18]将所有样品的重构信息整合在一起，再使用Cuffcompare软件将整合后的转录本与参考注释信息进行比较，计算转录本数量。

使用rMATS软件[19]检测样品中5种AS事件，分别为外显子跳跃(skipped exon，SE)、可变供体位点(alternative 5′ splicing site，A5SS)、可变受体位点(alternative 3′ splicing site，A3SS)、互斥外显子(mutually exclusive exons，MXE)和单内含子保留(retained intron，RI)，并检测肌肉和脂肪组织样品间的DSG。rMATS软件是一款用于RNA-seq数据检测DSG的软件，它可以计算inclusion亚型和skipping亚型的相对丰度(图1)，同时计算相应的P-value和FDR值来衡量样品间剪接差异的显著性。本研究中使用FDR≤0.05为判断标准，小于该阈值的DSG定义为显著差异的DSG。

I. 内含子异构体读长; S. 跳跃异构体读长；j. 连接长度；e1，e2. 外显子长度；r. 读长；lI. 内含子有效长度；ls. 外显子异构体有效长度I. Reads of the inclusion isoform; S. Reads of the skipping isoform; j. Junction length; e1，e2. Exon length; r. Read length; lI. Effective length of the inclusion isoform; ls. Effective length of the skipping exon isoform图1 AS事件中各种剪接类型的相对丰度计算Fig.1 Calculation of relative abundance of different splicing types in AS

1.5 肉鸡胸肌和腹部脂肪组织DSG的GO功能富集和KEGG分析

使用KOBAS 3.0软件[20-21]进行DSG的注释、GO功能富集分析、信号通路(pathway)分析，并对相关信号通路和基因进行进一步分析。进行GO功能富集分析时，以基因组为背景，统计每个GO分析中基因与背景基因之间的差异，对P值进行校正，定义校正后的P<0.05为这些基因的功能在该GO类别中显著富集。进行KEGG分析时，以基因组为背景，根据注释结果以及功能分类，将DSG进行生物通路分类，校正后的P<0.05时，表明DSG富集在一条通路中。同时对通路中DSG的表达量进行检测。具体方法为：与参考基因组比对之后，使用RSEM软件[22]计算基因和转录本的表达水平，采用校正的P≤0.05来筛选肌肉和脂肪组织基因的差异表达基因。

1.6 数据统计分析

使用IBM SPSS Statistics 19.0软件(SPSS Inc., Chicago, IL, USA)对肌肉和脂肪组织基因组中的可变剪接数量进行独立样本T检验和多重比较(SLD法)，数据结果以“平均数±标准差”表示。

2 结果

2.1 肉鸡肌肉和脂肪组织AS事件鉴定

测序原始数据经过质量控制后，肉鸡胸肌3个组织样分别获得6.69、6.71和6.75 G，腹脂3个组织样分别获得6.68、6.69、6.73 G可用于注释分析的序列片段(clean reads)。得到clean reads之后，使用HISAT软件[16]将clean reads比对到鸡参考基因组序列。肌肉组织样品比对率分别为69.18%、69.49%和68.47%，脂肪组织样品对比率分别为74.60%、73.10%和73.75%。相同组织样品间均匀的比对率表明样品之间的数据具有可比性。

肉鸡肌肉和脂肪组织基因组中AS事件鉴定结果见表1。由表1可知，肉鸡肌肉和脂肪组织基因组中分别检测到(5 966.00±111.66)和(6 757.00±156.51)个AS，脂肪组织中AS数量比肌肉组织中高13.26%(P<0.05)。共检测到的5种AS类型中，除MXE外，脂肪组织中4种类型剪接事件(SE、A5SS、A3SS和RI)的数量均显著高于肌肉组织(P<0.05)；且以SE为主，分别占肌肉和脂肪组织中总AS事件数量的54.92%和52.67%。

表1肉鸡肌肉和脂肪组织基因组中AS事件鉴定

Table1Theresultofalternativesplicing(AS)eventsidentifiedinmuscleandfattissuesgenomeinbroilers

项目Item肌肉组织Muscle tissue脂肪组织Fat tissueP值P-valueAS事件总数Total number of AS 5 966.00±111.666 757.00±156.510.002AS类型AS typeSE3 276.33±65.623 558.67±72.110.008A5SS579.67±20.23700.67±15.950.001A3SS989.00±13.861 205.33±39.460.001MXE276.33±11.06265.67±18.230.435RI844.67±22.501 026.67±35.300.002

2.2 肉鸡肌肉和脂肪组织DSG鉴定

经rMATS软件[19]检测分析发现，AA肉鸡肌肉和脂肪组织中检测到513个显著的DSGs，其中SE、A5SS、A3SS、MXE和RI 5种类型的DSG数量分别为314、34、17、54和94个(图2)。DSG中发生SE可变剪接的基因数量最多；其次为RI和MEX；发生这3种类型的DSG占总AS事件数量的90.06%。

2.3 肉鸡肌肉和脂肪组织DSG的GO注释

为了进一步了解肉鸡肌肉和脂肪组织DSG的功能，使用KOBAS 3.0在线软件[20-21]对所有DSG基因进行了基因注释、基因本体(gene ontology，GO)注释和功能分析。选取校正P≤ 0.05的DSG作为GO富集分析结果(图3)。通过GO富集分析发现，肉鸡肌肉和脂肪组织中513个DSGs中的93个DSGs显著富集在细胞组分(cellular component)和分子功能(molecular function)类别中的11个GO项内，同一基因多富集在多个GO项中。影响肉鸡肌肉和脂肪组织DSG涉及的细胞组分包括细胞(GO：0005623)、细胞部分(GO：0044464)、细胞内(GO：0005622)、细胞外(GO：0044424)、细胞器(GO：0043226)、胞内细胞器(GO：0043229)、细胞质(GO：0005737)、胞浆部分(GO：0044444)8个GO项；还有少量DSG分类到分子功能中，包括蛋白质结合(GO：0005515)、磷蛋白结合(GO：0051219)和蛋白质磷酸化氨基酸结合(GO：0045309)3个GO项。

图2 肉鸡肌肉和脂肪组织差异可变剪接基因数量Fig.2 The number of DSG between muscle and fat tissues in broilers

图3 肉鸡肌肉和脂肪组织基因组DSG的GO富集结果Fig.3 The result of GO analysis for DSG in genome between muscle and fat tissues in broilers

2.4 肉鸡肌肉和脂肪组织DSG的通路和相关基因分析

以鸡基因组为参照，使用KOBAS 3.0在线软件[20-21]对肉鸡肌肉和脂肪组织DSG进行了通路分析。分析发现，31个DSGs富集在肌动蛋白细胞骨架调节、焦点粘连、丙酮酸代谢和细胞内吞作用4条信号通路中。肌动蛋白细胞骨架调节与焦点粘连两条信号通路关系密切，其中10个基因均在两条信号通路之中(表2)。肌动蛋白细胞骨架在维持细胞形状、调控细胞粘连等方面具有重要功能；焦点粘连则参与细胞膜受体和肌动蛋白骨架之间的结构连接；丙酮酸代谢作为三羧酸循环的重要组成部分，在调解物质代谢中具有重要作用；细胞内吞作用在分子运输转运方面具有重要作用。

为了验证DSG在肌肉和脂肪组织中的表达情况，以脂肪组织为参照，经RSEM软件[22]分析，信号通路中的31个DSGs有13个在肌肉组织中上调表达，4个下调表达；其中14个DSGs(占信号通路基因中48.39%)表达量在肌肉和脂肪组织中未发生显著的表达量变化(表3)。

表2肉鸡肌肉和脂肪组织基因组中DSGs的通路分析结果

Table2TheresultofpathwayanalysisforDSGsingenomebetweenmuscleandfattissuesinbroilers

代谢通路PathwayDSG数量Number of DSGsKEGG ID校正的P值Corrected P-valueDSG名称DSG name肌动蛋白细胞骨架调节Regulation of actin cytoskeleton13048100.00CRK,ENAH,FGFR1,ITGA3,ITGA7,ITGB1,MRAS,MYL9,P38C-CRK,PPP1R12A,PPP1R12B,PXN,ROCK2焦点粘连Focal adhesion11045100.02CCND3,CRK,ITGA3,ITGA7,ITGB1,MYL9,P38C-CRK,PPP1R12A,PPP1R12B,PXN,ROCK2丙酮酸代谢Pyruvate metabolism5006200.03ACSS2,HAGH,LDHD,PDHA1,PKM细胞内吞作用Endocytosis12041440.04AP2B1,CAPZB,CCDC53,CLTCL1,CYTH1,EPN3,FAM21A,GIT2,H3KBP1,SPG21,USP8,VPS29

表3肉鸡肌肉和脂肪组织基因组中DSGs差异表达结果

Table3TheresultofdifferentialexpressionlevelsforDSGsingenomebetweenmuscleandfattissuesinbroilers

序号N基因名称Gene name基因IDGene ID差异表达倍数Fold change校正的P值Corrected P-value上/下调Up/Down-regulation1ACSS2419158-0.411.59×10-38Down2AP2B1417525-2.420.01*3CAPZB3964180.535.78×10-28Up4CCDC53418094-5.210.42*5CCND34199280.862.26×10-13Up6CLTCL14167650.364.86×10-53Up7CRK107054794-2.370.08*8CYTH1422085-1.132.78×10-9*9ENAH374180-246.060.99*10EPN31008592900.304.38×10-46Up11FAM21A423772-0.739.97×10-14Down12FGFR1396516-61.850.96*13GIT23740350.992.78×10-5Up14HAGH4165370.456.14×10-37Up15ITGA3373946-0.501.69×10-21Down16ITGA71017490980.414.32×10-21Up17ITGB1374058-0.994.90×10-8Down18LDHD4156890.473.38×10-10Up19MRAS3951490.720.00Up20MYL9396215-1.103.46×10-6*21P38C-CRK417553-1.580.60*

(转下页 Carried forward)

使用RSEM软件计算基因和转录本的表达水平，采用校正的P≤0.05来筛选肌肉和脂肪组织的差异表达基因；Up表示基因在肌肉组织中上调，脂肪组织中下调，Down表示基因在肌肉组织中下调，在脂肪组织中上调；*表示基因表达量差异不显著

The expression level of genes and transcripts were calculated using RSEM software. Genes with correctedP-value ≤ 0.05 were considered DEGs. Up indicted gene was up-regulated expression in muscle tissue and down-regulated expression in fat tissue; Down indicted gene was down-regulated expression in muscle tissue and up-regulated expression in fat tissue; * indicted the genes expression level were not different between muscle ant fat tissues

3 讨论

自1978年Gilbert[23]发现AS现象以来，Ensembl数据库[24]已经收录了人、小鼠、斑马鱼、鸡、鸭、猪、牛、羊等80多种生物的可变剪接事件。基于高通量RNA-seq技术挖掘基因组AS的应用[25]，使得人、小鼠、猪、牛、果蝇等多个物种基因组的可变剪接事件得到了全面的鉴定和分析[14]。冉茂良等[14]研究发现，AS事件具有发育时期特异性。徐铁山等[26]在北京鸭肌肉和皮脂中分别发现了5 958.00和7 528.67个AS。本研究中，在肉鸡肌肉和脂肪组织基因组中分别检测到5 966.00和6 757.00个AS，其中脂肪组织中AS数量显著高于肌肉组织。由此可见，可变剪接事情具有组织特异性，脂肪组织发生AS事件的频率远远高于肌肉组织。本研究发现，肌肉和脂肪组织中DSG发生SE可变剪接的基因数量最多。SE是选择性剪接最简单和最普遍存在的一种方式，这与Li等[27]对猪肌肉组织中差异剪接类型的研究相一致，说明SE可变剪接类型在肉鸡肌肉和脂肪组织中发生频率较高。

差异剪接分析对于理解细胞和器官的生理过程十分重要[28]。为了探索肌肉和脂肪组织生理过程的差异，本研究进行了肉鸡肌肉和脂肪组织间差异剪接分析。肉鸡肌肉和腹脂组织基本组成差异较大，前者由肌纤维构成，后者则属于白色脂肪，由脂肪细胞构成[29]。本研究肌肉和脂肪组织间的DSG的功能大部分与细胞基本结构相关。由此可见，差异剪接在调控肌肉和脂肪组织细胞的基本组成上具有重要作用。信号通路分析发现，DSG富集在肌动蛋白细胞骨架调节、焦点粘连、丙酮酸代谢和细胞内吞作用4条信号通路中，这些通路与肌肉组成、脂肪代谢等密切相关。由此可见，肌肉和脂肪组织细胞组成和行使功能的差异，部分源于基因组中差异可变剪接。可变剪接调控的基因通常不同于差异表达所调控的基因。因此，在基因组水平描述基因调控时，差异剪接与差异基因表达通常互相补充[15]。本研究发现，约50%的DSG在肌肉和脂肪组织中的表达水平未发生显著变化，说明肌肉和脂肪组织中DSG和差异表达基因发挥各自的功能，共同调控肌肉和脂肪组织的基本组成、生理功能等。

对于DSG可以通过剪接因子对整个基因组的可变剪接进行研究，也可以针对某一基因进行具体的结构和功能的研究。剪接因子在mRNA前体剪接过程中发挥重要作用。研究发现，可以通过剪接因子SRSF10对鸡基因组的剪接作用来研究整个基因组中的可变剪接[30]。在与本试验相关的研究中发现，对脂肪代谢具有重要调控功能的Apelin信号通路对肉鸡腹脂的沉积具有重要的调控作用(未发表)。其中Apelin前体由 77 个氨基酸组成，编码基因通过可变剪接可指导合成多种成熟的Apelin活性多肽片段：Apelin-36、Apelin-17、Pyr-apelin-13、Apelin-13和Apelin-12等[31]。因此，对Apelin基因的可变剪接、时空表达规律以及功能验证等研究，能够加深对肉鸡中该基因DSG、基因功能的深入理解。