绵羊不同发育阶段背最长肌组织中可变剪接的鉴定与分析

鲍晶晶 浦亚斌 马月辉 赵倩君

（中国农业科学院北京畜牧兽医研究所，北京 100193）

在真核生物转录过程中，pre-mRNA通过多种不同的方法切除内含子，重新拼接外显子产生成熟的mRNA，这个过程称为可变剪接（Alternative splicing，AS）［1］，可变剪接是高等真核生物蛋白质多样性和遗传多样性的重要调控机制［2］。人基因组中95%以上的多外显子基因能发生可变剪接［3］。单个基因通过可变剪接产生的mRNA异构体的数量从2到数千个不等［4］。可变剪接的基本类型包括外显子跳跃（Skipped exons，SE），5′端可变剪接位点（Alternative 5′ splice sites，A5SS），3′端 可 变 剪接位点（Alternative 3′splice sites，A3SS），互斥外显子（Mutually exclusive exons，MXE）和内含子保留（Retained introns，RI）［5］。通过可变剪接产生的mRNA和蛋白质异构体在结构、稳定性、亚细胞定位和功能上是不相同的，且一些基因的异常可变剪接会导致发育无法正常进行［6］。Blue等［7］在刚出生和成年小鼠心脏中鉴定了由可变剪接产生的3个特异性基因包涵体重链、包涵体轻链-a和转运蛋白结合蛋白-1，并表明包涵体重链基因的剪接调控被抑制时，小鼠的体重和肌肉重量就会增加。PKM基因的一种剪接体SRSF3能调节结肠癌细胞的生长和肿瘤特异性能量代谢的维持［8］。调节可变剪接事件是心肌细胞表达功能细胞骨架和肌节蛋白的必要步骤。大约3%的特发性扩张型心肌病患者存在RNA结合蛋白RBM20的突变，RBM20是可变剪接的组织特异性调节剂［9］。锌指蛋白ZNF148的两种剪接异构体在大肠癌细胞凋亡和侵袭中发挥不同作用，ZNF148能促进CRC细胞增殖、侵袭和迁移，ZNF148表现出相反的作用，由此得出ZNF148的两种剪接异构体可能对恶性生物学活性产生相互拮抗作用［10］。有研究表明约1/3的mRNA转录本会被无义衰变降解，以调节正常基因表达或靶向降解异常表达的基因［11］。这些研究表明AS形成的功能各异的转录本在肌肉组织发育及疾病发生过程中扮演着重要作用。

可变剪接事件在畜禽动物中的研究也陆续被报道。在木川乌骨鸡中鉴定了TYR基因的5个可变剪接异构体，并证明其在黑色素形成机制中扮演不同角色［12］。对小尾寒羊和陶塞特羊的卵巢组织中差异AS事件基因功能注释到与繁殖力、发育相关的通路［13］。Agouti是一个与毛色相关的基因，在白色山羊中发现了5种可变剪接异构体，而在黑色山羊中只发现了2种异构体［14］。许多涉及性腺发育、激素代谢、昼夜节律的基因通过选择性剪接被差异调控［13］。但对多浪绵羊肌肉不同发育过程中可变剪接的时序变化及其作用机制未见报道。本研究基于多浪绵羊不同发育阶段背最长肌的RNA-seq测序数据，通过rMATS软件［15］分析首次鉴定了不同发育阶段背最长肌组织中剪接事件和差异剪接基因（Differential spilicing gene，DSG），通过GO和KEGG分析这些DSG的功能，旨为研究肌肉不同发育阶段AS事件及其调控机制提供理论依据。

1 材料与方法

1.1 材料

本研究所用多浪绵羊来自中国农业科学院北京畜牧兽医研究所昌平试验基地。随机选择健康的妊娠90日龄的绵羊胎儿（F90）3只，出生后30日龄羔羊（L30）3只和3岁成年羊（A3Y）3只，屠宰后分别采集9只多浪绵羊的背最长肌组织立即保存于液氮中。

1.2 方法

1.2.1 组织总RNA的提取和质量检测 根据RNeasy® Plus Universal Mini Kit（Qiagen）试剂盒的操作说明进行绵羊背最长肌组织总RNA的提取。用1%的琼脂糖凝胶电泳检测总RNA是否降解以及是否有污染，用Nanodrop ND-1000检测RNA的纯度，Agilent 2100 bioanalyzer精确检测RNA的浓度和完整度（RIN），将RNA浓度≥100 ng/μL，总量≥3 μg，RIN值≥7.0的RNA样品用于转录组建库测序。

1.2.2 cDNA文库的构建和测序 9个绵羊背最长肌组织的cDNA文库构建和测序均由北京贝瑞和康生物技术有限公司完成。测序平台是Illumina Hiseq2500 V4，测序模式为125PE。用NGSQC Toolkit v2.3.3软件［16］过滤掉含有接头的reads，过滤掉N含量大于10%的reads，过滤掉低质量的reads，得到 clean reads。

1.2.3 RNA-seq数据的比对与组装 从NCBI数据库中下载绵羊的参考基因组（Oar_v4.0）以及基因组注释GTF文件（http：//www.ncbi.nlm.nih.gov/genome/）。将clean reads 通过HISAT［17］软件比对到绵羊的参考基因组上，用StringTie［18］软件将比对上的reads进行转录本的组装并计算每个基因和isoform的表达水平，然后通过merge［19］将所有样品的转录本组装在一起再次呈递给StringTie软件重新估算转录本丰度，然后将数据载入ballgown包［20-21］进行差异基因分析。

1.2.4 可变剪接的鉴定及差异分析 使用rMATS［15］软件对每个样品存在的可变剪接类型及相应表达量进行检测。rMATS是一款可利用RNA-Seq数据自动检测和分析差异可变剪接的软件。rMATS软件的统计模型是根据计算P值和错误发生率（False discovery rate，FDR）来鉴定差异的可变剪接。本实验采用FDR＜0.05作为标准来筛选差异剪接基因。

1.2.5 差异剪接基因的功能富集分析 对差异剪接基因采用 g：profiler（https：//biit.cs.ut.ee/gprofiler/）在线分析软件［22-23］进行GO富集分析，采用DAVID（https：//david.ncifcrf.gov/）在线分析软件［24-25］进行KEGG富集分析，采用P＜0.05作为显著富集标准。

2 结果

2.1 可变剪接事件的鉴定

RNA-seq原始数据经过质控后，每个背最长肌组织样品均产生不低于10 GB的clean data，GC含量平均为48.6%且Q30＞87%，测序数据满足后续试验要求。通过HISAT软件将clean data比对到绵羊参考基因组上，9个样品83.8%-88.5%的clean reads能比对到绵羊的参考基因组上。将经过Stringtie软件组装的转录本与绵羊参考基因组进行比较，绵羊妊娠90日龄胎儿时期（F90）、出生后羔羊时期（L30）和成年3岁时期分别鉴定出13 923、11 959和12 164个可变剪接事件（图1）。其中SE占各个发育阶段鉴定出的可变剪接类型的比例最高，达到69.19%-72.19%，A5SS、A3SS、MXE、RI分别占各个阶段可变剪接事件总量的6.80%-7.75%，1.04%-1.19%，5.80%-6.12%，4.49%-5.18%。从图中可以看出不同发育阶段的绵羊背最长肌组织中的可变剪接事件的数量存在差异，且在胚胎时期的数量最多，出生后30 d和成年3岁时期的变化不大，这也与绵羊出生前后肌肉的发育相一致。

图1 绵羊背最长肌组织不同发育阶段可变剪接事件的数量

2.2 背最长肌组织的差异剪接基因鉴定

对rMATS软件分析结果以FDR＜0.05为标准进行差异筛选，共筛选到6 935个显著的差异剪接基因（图2），F90与L30比较组在SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为1 856，155，310，199，25；F90与A3Y比较组在SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为1 958，165，302，236，28；L30与 A3Y比 较 组 在 SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为1 216，120，184，152，29。差异剪接基因中SE事件的占比最高（约72.41%），RI的占比最少。

图2 绵羊背最长肌组织不同发育阶段差异可变剪接基因数量

2.3 差异剪接基因的功能富集分析

对获得的差异剪接基因进行GO富集分析，多浪绵羊胎儿组与羔羊组背最长肌组织差异剪接基因富集到的前20个GO terms（图3），多浪绵羊胎儿组与成年羊组背最长肌组织差异剪接基因富集到的前20个GO terms（图4），多浪绵羊羔羊组与成年羊组背最长肌组织差异剪接基因富集到的前20个GO terms（图5）。3个比较组都富集到了与肌肉显著相关的GO terms，如横纹肌发育（GO：0014706）、肌肉结构发育（GO：0061061）、肌细胞分化（GO：0042692）、肌膜（GO：0042383）、发育进程（GO：0032502）。

图3 绵羊胎儿组与羔羊组肌肉组织差异剪接基因的GO富集结果

对3个比较组的差异剪接基因的KEGG通路分析也显著富集到与肌肉显著相关的通路，见表1-表3，如胰岛素信号通路、Wnt信号通路、MAPK信号通路、肌动蛋白细胞骨架调控。其中胞吞作用、胰岛素信号通路、焦点粘连、泛素介导的蛋白水解是这3个比较组中都显著富集到的共同信号通路。细胞的胞吞作用参与分子物质运输转运；胰岛素信号通路参与肌肉的发育；焦点粘连参与细胞膜受体和肌动蛋白骨架之间的结构连接；泛素介导的蛋白水解调控细胞周期；肌动蛋白细胞骨架调控维持细胞形态和调控细胞粘连。

图4 绵羊胎儿组与成年羊组肌肉组织差异剪接基因的GO富集结果

3 讨论

1978年，Gilbert第一次发现可变剪接以来［26］，已证明可变剪接参与高等真核生物的基本生命调节过程并已成为真核基因调控领域的研究热点。可变剪接是转录后基因表达调控的关键步骤，它通过使单个基因产生不同的RNA异构体而有助于蛋白质和功能的多样性。可变剪接通过细胞和组织特异性调节细胞分化和机体发育，可变剪接的失调会导致人类各种遗传疾病，包括早衰症、强直性营养不良、额颞叶痴呆、囊性纤维化和癌症［27］，并且可变剪接最近已成为预后和治疗癌症的标志［28］。可变剪接通过控制特定细胞类型和特定时间来表达特定的RNA异构体［29］。LPIN1基因有2种RNA剪接体，LPIN1α主要作用于前脂肪细胞的分化，而LPIN1β则会促进前脂肪细胞增殖［30］。在陕北绒山羊睾丸组织中首次鉴定了CTNNB1的4个转录本异构体，命名为Ctnnb1-A/-B/-C/-D，Ctnnb1-C在睾丸中表达最高，表明其在雄性生育中起重要作用［31］。NFIX是NFI家族的成员，在肌肉和大脑发育中发挥重要作用，有研究在奶山羊中鉴定并验证了NFIX基因的2个新转录本（NFIXa和NFIXb），正常NFIX转录本在肺中高表达，NFIXa异构体在胰腺中高表达，而NFIXb异构体在肺和胰腺中高表达。此外，NFIXa异构体在肝脏、脾脏、脂肪、肠道和睾丸中的表达水平显著高于正常NFIX和NFIXb异构体，NFIXa异构体在大脑中的表达显著高于NFIXb异构体［32］。说明可变剪接事件具有时空特异性和组织特异性。高通量转录组测序技术已成为研究RNA转录产物的可变剪接的可靠工具。冉茂良等［33］对猪60胚龄、90胚龄、30日龄和180日龄4个不同发育时期的睾丸进行转录组测序鉴定出92 738个可变剪接事件，筛选到63个发生可变剪接的基因与睾丸素代谢密切相关。张敏等［34］利用RNA-seq技术在肉鸡肌肉与脂肪组织中分别鉴定到5 966和6 757个AS事件，并检测到513个参与脂肪和肌肉发育的显著差异的剪接基因。郭家中等［35］在简州大耳羊早期胎儿（妊娠45 d、60 d、105 d）到新出生（出生后3 d）阶段背最长肌中鉴定了137 308个可变剪接事件，4个阶段均发生可变剪接事件的基因经GO富集分析表明显著富集在细胞表达调控和物质合成等过程，并且发现肌肉早期发育阶段发生的剪接事件多与晚期时期。黄艳群等［36］研究发现与丝羽乌骨鸡趾型显著相关的是Lmbr1基因转录本Lmbr1-α。

图5 绵羊羔羊组与成年羊组肌肉组织差异剪接基因的GO富集结果

表1 绵羊胎儿组与羔羊组肌肉组织差异剪接基因的KEGG通路分析

表2 绵羊胎儿组与成年羊组肌肉组织差异剪接基因的KEGG通路分析

表3 绵羊羔羊组与成年羊组肌肉组织差异剪接基因的KEGG通路分析

基因组中差异可变剪接会导致细胞组成和功能的差异。本研究采用rMATS软件鉴定可变剪接事件及差异剪接基因，在妊娠90日龄与出生后30日龄、妊娠90日龄与成年3岁绵羊以及出生后30日龄与成年3岁绵羊背最长肌组织比较组中分别鉴定到了2 545，2 689和1 701个差异剪接基因。从研究结果可以看到出生前后背最长肌组织的可变剪接事件存在差异，由此可见可变剪接具有发育阶段特异性，且可变剪接在出生前的行动更活跃，这是与哺乳动物出生前主要是肌纤维数量的增加，出生后则是肌纤维的增长和变粗有密切关系［37］。本研究可变剪接事件中SE事件的数量是最多的，这与前人的研究一致［38］，说明SE是生物体中最普遍发生的可变剪接事件。对差异剪接基因进行功能富集分析，GO显著富集到了与肌肉发育相关通路上且共同富集到细胞内部分（GO：0044424），参与细胞内物质合成，KEGG共同富集到了胞吞作用、胰岛素信号通路、焦点粘连、泛素介导的蛋白水解4条信号通路上，这些信号通路与肌肉发育密切相关，表明可变剪接在绵羊背最长肌的发育过程中发挥重要作用。差异可变剪接基因与差异表达基因不一样，本研究中约26%的差异可变剪接基因的表达水平是差异的，说明差异可变剪接基因与差异表达基因各自发挥功能共同维持细胞生命活动［39］。

4 结论

多浪绵羊背最长肌组织在妊娠90日龄胎儿时期（F90）、出生后羔羊时期（L30）和成年3岁（A3Y）时期分别鉴定出13 923、11 959和12 164个可变剪接事件，其中SE占各个发育阶段鉴定出的可变剪接类型的比例最高，约占70.69%，A5SS、A3SS、MXE和RI分别约占7.28%、11.18%、5.96%和4.84%。绵羊背最长肌组织在F90_vs_L30、F90_vs_A3Y和L30_vs_A3Y比较组中鉴定到差异剪接基因数分别为2 545、2 689和1 701个。本研究中3个比较组的差异剪接基因的GO和KEGG富集分析都显著注释到与肌肉发育显著相关的通路上。