基于全长转录组测序的牛脂肪组织可变剪接事件分析

夏晗，李凡，彭玲伟，杜玉琴，郭爱珍，杨利国，周扬

华中农业大学动物科学技术学院、动物医学院/农业动物遗传育种与繁殖教育部重点实验室，武汉 430070

肌内脂肪的含量以及沉积分布影响牛肉的嫩度、系水力、剪切力值、风味和多汁性，适当提高牛肉肌内脂肪含量可改变牛肉品质[1-2］。但是，当皮下脂肪的积累较多时会降低饲料报酬，造成资源损耗[3］。脂肪沉积受复杂的调控网络影响，Utrera 等[4］针对包括脂肪沉积能力的牛胴体性状遗传力进行了72 项研究，也未能完全了解其遗传机制。

可变剪接是剪接体对前体mRNA 选择性的修饰，使同一基因可以表达不同成熟mRNA 剪接异构体的现象[5］。可变剪接被认为在基因中广泛存在，有研究显示，在高通量测序结果中95%以上的人类基因存在可变剪接，由此产生表达多种不同的剪接异构体[6］。可变剪接使基因表达具有更多可能性，展示了基因功能的多样性[7］。研究显示可变剪接的发生会导致蛋白结构的变化，导致同一基因的不同蛋白亚型在功能上出现差异[8］；且有研究表明，可变剪接会对牛的脂肪相关基因表达产生重要影响[9］。例如，NOVA 依赖性的可变剪接会调节白色脂肪组织褐变[10］、质量减轻导致的皮下脂肪中TCF7L2 可变剪接调节可导致脂肪组织中的高血糖和胰岛素作用受损[11］等。所以了解可变剪接事件的发生对于明晰脂肪沉积调控网络有着重要的作用。

基于PacBio 平台的全长转录组测序是第三代测序技术的代表，能直接测序得到RNA 转录本的全长片段，可以规避短片段带来的缺失信息的影响，提供更完整的转录组信息，更好地进行差异表达基因分析及功能注释[12］。同时，长读技术有利于更精确、更准确、更多地发现可变剪接事件，有助于提高对脂肪沉积调控网络的认识[13］。目前，第三代测序技术已经被广泛应用于植物的遗传信息挖掘，并取得了较好的结果，例如：赪桐的种质资源挖掘[14］、紫娟茶树叶片呈色机理及物质合成途径[15］、掌叶鱼黄草的遗传信息挖掘[16］等。但是，目前基于全长转录组的研究较少，尤其是牛类动物繁殖周期长、生长缓慢更是限制了其研究的进展。本研究利用基于PacBio 测序平台的第三代测序技术对夷陵黄牛脂肪组织遗传信息进行挖掘，以期为正确认识可变剪接对脂肪沉积调控网络的机制提供理论参考。

1 材料与方法

1.1 试验材料

供试牛为湖北省宜昌地区、年龄在3 岁、发育正常、体态良好的3 头健康夷陵黄牛。屠宰后，使用消毒手术刀，分别在每头牛相同部位采集5 cm3大小的皮下脂肪、腹腔脂肪、肌间脂肪样品。使用75%乙醇彻底清洗组织，PBS 清洗残留在组织上的乙醇，吸干水分后于−80 ℃液氮保存，用于总RNA 提取。

1.2 总RNA提取

使用TRIZOL 法提取夷陵黄牛脂肪组织总RNA。提取后的RNA 需使用1%凝胶电泳检测样品纯净度，并经过Nano Drop 仪器检测，确定总RNA的OD260/OD280值在1.9~2.1、OD260/OD230值在1.4~1.8。并使用Qubit 对RNA 进行定量，Agilent 2100、Bioanalyzer 等进一步进行质量检测，提取的总RNA质量满足要求后，低温保存，待后续进行文库构建。

1.3 文库构建及测序

构建文库所使用的RNA 需要使用Qiagen 试剂盒进行回收纯化预处理，以获得更高质量的RNA。随后使用SMARTer®PCR cDNA Synthesis Kit 试剂盒反转得到cDNA。依据cDNA 质量优化扩增程式，使用KAPA HiFi PCR Kits 试剂盒进行扩增，对获得的全长cDNA 片段使用SMRT bell template prep kit 1.0试剂盒构建文库。同时，对cDNA 进行末端修复、损伤修复处理，并在全长cDNA 的双端添加stem loop sequencing adapter 进行筛选，文库构建成功后，利用PacBio Sequel 测序平台开展全长转录组测序。

1.4 全长转录组数据组装与筛选

对建库测序得到的原始数据，首先使用SMRTLINK 5.1 软件对其进行质控，然后使用SNR（Signal Noise Ratio）筛选剔除低质量数据及接头序列，得到环形一致序列（circular consensus sequence，CCS）。去除含5' primer、3' primer、poly A 片段，获得的全长转录组序列经过校正和去冗余后，可用于后续分析。

1.5 可变剪接检测

试验使用ICE（iterative isoform-clustering）中的algorithm 模块进行聚类，使用DAGCon 得到一致性序列，并通过arrow 进行校正。将二代数据校正后的全长序列与参考基因组比对，去除融合基因及非注释数据后，将比对到参考基因组同一基因上的不同isoform 进行提取，即为可变剪接。根据可变剪接不同的形成方式进行分类。

1.6 全长转录组功能注释与分类

获得高质量的全长转录组序列后，使用KOG 蛋白质真核同源数据库、GO 基因本体数据库、KEGG代谢通路数据库、NR 非冗余蛋白数据库、Swiss Prot精准蛋白数据库对预测的Isoforms进行注释，用以评估脂肪组织中基因的功能注释信息及其参与的代谢通路及调控网络，并判断基因的物种相似性。

2 结果与分析

2.1 测序质量分析

对采集的夷陵黄牛的皮下脂肪、腹腔脂肪、肌间脂肪样品使用PacBio 平台第三代测序技术进行全长转录组测序，共获得19.81 G的原始数据，筛除接头和长度小于50 bp的原始离线数据后，得到13 014 954条Subreads，平均长度为1 474 bp（Full Passes≥1、最小预测准确度0.8），核验下机数据中的CCS（circular consensus sequence）的质量，共获得779 345 个CCS reads，CCS 总碱基数为1 812 898 181 bp，CCS Read平均长度2 320.43 bp，平均Full Pass 数量为11。通过验证5′引物、3′引物和poly-A尾的保留，获得高质量的FLNC（full-length nonchimeric reads）640 205条，占比82.15%，平均长度为2 091 bp，以上结果说明全长转录组测序数据质量较好，可以用于后续试验。

2.2 全长转录组Reads 聚类与校正

利用来源于相同样品的二代测序数据对PacBio平台的第三代测序reads 进行校正[17］。校正数据后，我们获得了779 345 个consensus reads。使用ICE 去除缺失序列和冗余序列后，进一步使用Lordec 软件，以二代数据作为参照对356 772 条一致性序列结果进行校正，获得可用于进一步分析的一致性序列总长为790 043 023 bp。

2.3 可变剪接分析

校正后的数据与参考基因组进行比较，得到与2 095 个融合基因相关的2 768 条reads。将剩余数据进一步过滤去冗余。使用Match Annot 软件将比对后结果与注释信息进行比较，并将三代注释结果和原基因组结果进行合并，得到基因注释的PB iso⁃form 63 400 条，非基因区序列有36 480 条。将比对到参考基因组cDNA 区域的序列进行检测，共识别到可变剪接事件50 520 个。其中2 406 个基因发生外显子跳跃、1 184 个基因发生外显子替换、1 262 个基因发生外显子接受、4 085个基因发生内含子保留、2 365个基因发生外显子位置替换、4 143个基因发生其他可变剪接事件（表1）。

表1 可变剪接类型统计表Table 1 Statistical table of alternative splicing types

将检测到可变剪接事件的基因进行GO 富集分析，显示与83 个条目显著相关（P<0.05），其中富集到分子功能类的GO 条目有22 个，细胞组分类的GO条目14 个，生物过程类的GO 条目47 个。在细胞组分类别中最显著是“薄膜涂层”，其次是“涂层膜”“黏着斑”条目。在分子功能类别中最显著的是“分子功能”，其次是“连接”“蛋白质结合”条目。在生物学过程类别中最显著的是“生物调节”，其次是“细胞过程的调节”“生物过程的调节”，另外，我们还检测到与脂质合成及代谢相关过程的显著富集，包括“糖异生”“葡萄糖代谢过程”“己糖生物合成过程”“单糖生物合成过程”“糖酵解过程”“丙酮酸代谢过程”“己糖代谢过程”等。 这说明可变剪接对于脂肪沉积调控网络可能有着重要的影响（图1）。

图1 GO及KEGG功能富集图Fig.1 Functional enrichment of GO and KEGG

对检测到可变剪接事件的基因进行KEGG 富集分析，以期研究可变剪接对基因之间互作关系的影响，结果显示有27条通路显著富集，其中直接与脂质合成及代谢相关的有15条通路，例如“脂肪细胞中脂解的调节”“脂肪消化吸收”“脂肪细胞因子信号通路”“PPAR 信 号 通 路”“AMPK 信 号 通 路”“PI3KAkt 信号通路”等。这说明可变剪接的参与导致脂肪沉积调控网络的复杂化。

2.4 新基因功能注释分析

为了获得更直观的注释结果，使用NR、KOG、GO、KEGG 和Swissport 5 个数据库进行注释，在全部80 756 条序列中有69 259 条被注释，如表2 所示。如图2 显示至少有16 572 个基因被5 个数据库所注释，有26 653个基因至少在1个数据库中被注释。其中占全部基因数量81.61%的65 902个基因被NR 数据库所注释，根据序列同源性比对结果，其中28 530条与普通牛同源、10 028条与牦牛同源、3 265条与平原野牛同源、2 970 条与绵羊同源、2 431 条与亚洲水牛同源以及与其他物种同源的18 678 条，如图3 所示。这说明了夷陵黄牛在品种形成过程中主要受到普通牛和瘤牛的影响。

图2 功能注释韦恩图Fig.2 Venn diagram of functional annotation

图3 NR分类注释统计图Fig.3 NR classification annotation statistics

表2 功能注释统计表Table 2 Function note statistics table

KOG 数据库注释到41 004 条序列在基因区域，占全部序列的50.78%。根据序列信息可以分为26组，其中一般功能预测条目注释到8 406 条序列，信号转导机制条目7 221 条，注释最少的为未知蛋白组（仅为52 条）。从分组可以发现，注释到细胞代谢发育等通用途径的序列占绝大多数。例如：胞内运输、分泌和囊泡运输、脂质运输和代谢、能量生产和转化分别注释到2 822、1 620、920 条序列，这些途径证明了在夷陵黄牛脂肪组织中进行着高频率的能量代谢以及物质跨膜运输，这可能与脂肪细胞中脂质积累相关（图4）。

图4 KOG 分类注释统计图Fig.4 KOG classification annotation statistics

GO 数据库注释到的基因有26 653 个，占全部基因数量的33%，富集到GO条目共1 280个，在分子功能条目上富集到559 个条目，30 672 条序列，其中以蛋白质结合、核酸结合条目聚集到最多的序列，分别为5 123 条和2 706 条；生物过程富集到507 个条目15 039条序列，其中转录调控及蛋白磷酸化聚集序列最多，分别为1 495 条和1 423 条；细胞组分富集到214 个条目7 907 条序列，其中膜和膜的组成部分条目聚集的序列最多，分别为1 142 条和1 721 条。其中本研究检测到与细胞能量储存与代谢等相关条目的显著富集，可能与脂肪组织复杂的生物学过程及活跃的细胞代谢有关（图5）。

图5 GO分类注释统计图Fig.5 GO classification annotation statistics

KEGG 的功能注释结果发现，夷陵黄牛脂肪组织全长转录组数据中共注释到基因45 984 个，占全部基因的56.94%，KEGG 将注释得到的基因信息分为5 个大类：细胞过程、环境信息处理、遗传信息处理、代谢、有机系统，这些类别分别富集到10 745、9 516、5 789、7 005、1 655 个 基 因，如 图6 所 示。KEGG 除富集到“细胞生长和死亡”“细胞运动”“免疫系统”“神经系统”等通路外，还富集到如“能量代谢”“聚糖生物合成和代谢”“脂质代谢”等与脂肪合成、能量代谢高度相关的路径，这些显著富集的通路可能参与脂肪沉积网络的调控。

图6 KEGG 分类注释统计图Fig.6 KEGG classification annotation statistics

将80 756 条序列与Swiss-Prot 数据库蛋白库进行比对，共注释得到69 259 个基因，占全部基因数量的85.76%。采用angel 软件，对序列做CDS 和Pro⁃tein 预测，共发现94 458 条CDS 序列及其对应的pep序列，并且比对到UTR 非编码区域的155 214 条序列。将数据与基因注释结果对比优化后，得到3'or 5'延长区域共29 448 处，分别对应21 145 条reads，7 176 个基因。

3 讨 论

可变剪接已经被证实在95%以上的人类基因中存在，可变剪接的存在使得基因的功能存在不确定性，脂肪调控网络变得更为复杂[7］。三代测序技术可以弥补二代测序序列短、完整性差的缺陷，直接获得全长序列，从而准确地表征基因结构特点，同时使用RNA-seq 对全长转录组数据进行校正，在进一步提高PacBio 循环共有序列（CCS）质量的同时，降低长序列的高碱基错误率，提升可变剪接检出率和准确率[18］。本研究利用三代测序技术对夷陵黄牛脂肪组织（皮下脂肪、腹腔脂肪、肌间脂肪）的总RNA 进行全长转录组测序。通过数据质控和分析，一共确定了779 345 个CCS reads，平均长度为2 320.43 bp，共识别可变剪接事件50 520 个，超过在牛上发现全部可变剪接事件的33.4%，说明可变剪接可能与脂肪沉积复杂的调控网络紧密联系[19］。前5%的发生可变剪接基因GO 和KEGG 富集分析结果显示，83 个条目显著富集，其中与脂质代谢相关过程的的有“糖异生”“葡萄糖代谢过程”“糖酵解过程”“丙酮酸代谢过程”“己糖代谢过程”等。KEGG 结果显示，这些可变剪接相关基因涉及到的通路共有27 条，与脂质合成及代谢相关的15 条，其中“PPAR 信号通路[20］”“AMPK 信 号 通 路[21］”“PI3K-Akt 信 号 通 路[22］”“MAPK 信号通路[23］”“Ras 信号通路[24］”“Rap1 信号通路[25］”“ECM-受体相互作用通路[26］”“cGMP-PKG信号通路[27］”“Apelin信号通路[28］”等均有文献报道。这些可变剪接事件的发现，是对现有脂肪可变剪接数据库的补充，也为充分正确认识脂肪沉积调控网络提供了新的理论依据补充。

本研究共发现80 756 个新基因序列，KOG、KEGG、NR、Swiss Prot、GO 数据库分别注释到其中41 004、45 984、65 902、57 880、26 653条序列，在物种同源性比对中发现夷陵黄牛主要与普通牛和瘤牛同源，属于南方黄牛的血统组成，与夏小婷等[29］的研究结果一致。但是小众动植物普遍存在基因组信息研究匮乏等现象，在进行同源比对时可参考的基因组序列信息有限，容易导致比对结果的偏离[30］。尽管已经证明了全长转录组在缺少基因组信息时有获取数据的优势，但长读技术依旧受到原始样品质量的影响[31］。本研究中，KOG 数据库注释到的序列占全部序列的50.78%，其中聚集到一般性功能预测分组的序列数量最多，说明维持细胞一般性功能的序列占据主要部分。GO 数据库注释的序列占比33%，主要与维持细胞正常形态及功能相关，其中分子功能组以蛋白质结合、核酸结合条目聚集到最多的序列，生物过程组中转录调控及蛋白磷酸化聚集序列最多，细胞组分中膜的组成部分聚集的序列最多。KEGG 将注释得到的基因信息分为5个大类：细胞过程10 745 条、环境信息处理9 516 条、遗传信息处理5 789 条、代谢7 005 条、有机系统1 655 条，除富集到“细胞生长和死亡”“细胞运动”“免疫系统”、“神经系统”等主要通路外，还富集到如“能量代谢”“聚糖生物合成和代谢”“脂质代谢”等与脂质合成及代谢高度相关的路径。除此之外，我们还使用Swiss-Prot数据库蛋白库注释到69 259 条序列，发现94 458 条可能的CDS 序列及其对应的pep 序列。这些数据可以为现有脂肪沉积相关研究进行数据补充，并为后续研究提供参考信息，为后续的研究提供宝贵资源。