豹猫转录组从头组装及组织特异性表达分析

蒋兰， 张雪艳， 王俊茵， 李静

(1. 四川大学生命科学学院，生物资源与生态环境教育部重点实验室，成都610065；2. 峨眉山景区管委会，峨眉山生物多样性保护研究所，四川峨眉山614200)

豹猫Prionailurusbengalensis属食肉目Carnivora猫科Felidae豹猫属，地理分布极为广泛，从亚洲南部和西南部的热带雨林到中国北方和俄罗斯远东地区的森林(Hughes，1998；Woolfetal.，2002)，在中国，根据地理分布豹猫被分为5个亚种(Rossetal.，2015)。尽管豹猫资源丰富，但人类活动使其栖息地不断缩小，同时对其皮毛的需求也导致豹猫的生存受到严重威胁(Hughes，1998)。目前中国的野生豹猫种群数量正在下降，分布范围也逐渐缩小，一些原有的分布区域，如宁夏、青海、山东等地区的豹猫数量已非常稀少。《中国濒危动物红皮书》(汪松，1998)将豹猫列为易危种。

目前国内外有关豹猫的文献报道集中在遗传多样性和系统发育等研究。研究人员分别利用12s RNA、cytb序列和cDNA序列探讨西表岛猫Prionailurusiriomotensis和豹猫系统发育关系，结果表明二者分歧时间较近，西表岛猫应为豹猫的一个亚种(Masudaetal.，1994；Suzukietal.，1994)。Masuda和Yoshida(1995)基于cytb序列认为马岛猫Prionailurusbengalensiseuptilurus也是豹猫的一个亚种。Saka等(2018)发现西表岛和马岛的豹猫种群的MHC多样性远低于家猫Feliscatus种群，暗示岛屿上的豹猫由于近亲繁殖或地理隔离，其对病原体的抵抗力降低。Ko等(2018)基于微卫星分析了韩国豹猫的遗传多样性，发现其平均等位基因数和杂合度分别为3.8和0.41，多样性低于其他12种猫科动物57个种群的平均值。白素英等(2004)利用随机扩增多态性DNA(RAPD)研究中国6个豹猫种群结构发现，各种群的关系与种群间的地理距离及亚种划分一致。目前国内关于豹猫的研究集中在分类、分布、食性以及解剖等方面(张淑云等，1988；白素英等，2004；Xiongetal.，2016)，而关于其基因表达的研究还未见报道。

转录组测序技术(RNA-seq)是一种广泛应用于研究基因表达调控的分析方法，通过有效识别和量化RNA转录本，从而获得其特定生理状况的转录组信息(Fengetal.，2012；Spradlingetal.，2013)。RNA-seq使得从整体水平分析动物基因结构及基因功能成为可能，更有利于发现生物学过程、免疫反应基因、疾病产生过程中的基因表达谱变化以及免疫受体多样化机制(Wang & Kirkness，2005)。由于RNA-seq技术对样本新鲜程度要求较高，一定程度上阻碍了其在野生动物相关研究中的应用。本研究基于1只意外死亡的野生豹猫个体，采用RNA-seq对其大脑、心脏、肾脏、肝脏、肺和骨骼肌 6个组织进行转录组测序，通过从头组装、注释，报道了一个高质量的豹猫转录组。这些数据为理解豹猫的遗传背景及基因表达谱提供了宝贵的资源，也为这一小型猫科动物的基因组注释和深入开展保护遗传学研究奠定了基础。

1 材料与方法

1.1 样品采集与测序

豹猫各器官的组织样品均来自于1只在峨眉山被汽车意外撞死的雌性成年个体，被管理人员发现后提供给本实验室。共采集了6个不同器官的组织样品，分别为大脑、心脏、肾脏、肝脏、肺和骨骼肌，新鲜的组织样品迅速置于液氮中保存，其余样品-80 ℃储存。样品的总RNA提取和转录组建库测序工作由北京诺禾致源生物信息科技有限公司完成，采用Illumina HiSeq 4000测序平台对6个样品的cDNA文库分别进行150 bp的双端测序。

1.2 质控与过滤

转录组组装前需对原始序列进行质量控制，过滤标准为去除含N比例大于10%、低质量(质量值<5)碱基数比例大于50%、含有接头污染的序列，最终得到高质量的clean reads用于后续分析。

1.3 转录组从头组装

使用Trinity(Haasetal.，2013)对过滤后的clean reads进行denovo组装，仅输出片段长度大于300 bp的contig序列。由于Trinity组装获得的转录组中通常含有大量相似的冗余序列，故最后使用Cd-hit-est (Li & Godzik，2006)去除冗余序列，得到非冗余转录本用于后续分析。

1.4 功能注释

为了获得较为全面的基因功能信息，利用blastx(Camachoetal.，2009)将Unigene与NCBI非冗余蛋白质(NR)数据库、euKaryotic Ortholog Groups(KOG)数据库、Swiss-Prot数据库进行比对，E值设为1E-5。再根据NR数据库的比对结果进行Gene Ontology(GO)和Kyoto Encyclopedia of Genes and Genomes(KEGG)注释的分类统计，分别通过在线工具Web Gene Ontology Annotation Plot(WEGO)(Yeetal.，2006)和KEGG Automatic Annotation Server(KAAS)完成(Moriyaetal.，2007)。在使用WEGO进行分类统计之前，需先将转录本与NR数据库比对的结果导入Blast2GO(Conesaetal.，2005)搜索相关的GO功能条目。

1.5 转录本的表达量和组织特异性分析

2 结果

2.1 测序和质控

共获得171 370 884条150 bp的raw reads，共计51.4 Gb。质控过滤后得到154 489 126条(90.15%)clean reads，共计45.42 Gb，每个组织的数据量为5.54～13.46 Gb(表1)，表明测序数据质量符合要求。

表1 测序结果和质量控制Table 1 Summary of sequencing and quality control

2.2 转录组从头组装

使用Trinity将clean reads进行denovo组装，最终得到369 246条转录本，转录本长度301～31 696 bp，平均长度1 465 bp，Contig N50为2 660 bp(表2)。拼接的序列中，155 985条(42.24%)序列长度在1 000 bp以上，组装质量较高。大多数转录本长度为300～400 bp和1 000～2 000 bp(图1)。使用Cd-hit-est对Trinity组装的转录本去冗余，最终得到269 412条非冗余转录本用于注释分析。

表2 豹猫转录组组装结果Table 2 Summary of transcriptome assembly of Prionailurus bengalensis

2.3 转录本的功能注释

使用blastx将去冗余后的转录本与NR、Swiss-Prot、KOG和KEGG 4个公共数据进行同源性比对注释。NR蛋白数据库注释到转录本110 997条(41.13%)；Swiss-Prot蛋白数据库注释到转录本 98 192条(36.39%)；KOG数据库注释到转录本97 616条(36.17%)；KEGG数据库获得注释信息转录本65 895条(24.42%)。4个数据库注释结果的韦恩图显示，共注释到转录本114 517条(42.44%)，63 860 条转录本在所有数据库中都有注释(图2：A)。

NR数据库比对表明，73%的序列有较强的同源性(E值≤1.0E-30)，其中，25%的序列比对E值为0，39%的为0～1.0E-45，9%的为1.0E-45～1.0E-30(图2：B)。相似度分布显示，83%的序列比对相似度≥80%，仅有17%的为30%～80%(图2：C)。物种间比对结果显示，豹猫的转录本与猫科动物的相似性高，约49%的序列与猫科动物基因组(家猫：23%，金钱豹Pantherapardus：13%，猎豹Acinoyxjubatus：8%，虎Pantheratigris：5%)己经注释的基因同源(图2：D)。

转录本的GO富集显示，33 262条转录本注释到247 376个GO功能条目(图3：A)，包括生物过程(113 399条，45.84%)、分子功能(37 773条，15.27%)和细胞组分(96 204条，38.89%)。生物过程中，细胞过程(20 431条)、代谢过程(16 659条)和生物调节(12 227条)最多；细胞组分中，细胞(18 428条)、细胞部分(18 261条)和细胞器(13 859条)最丰富；分子功能中，最丰富的为连接(18 061条)和催化活性(11 917条)。

根据同源性比对KOG数据库，97 616条(36.17%)转录本分别在25种类别中被注释到(E值≤1E-5)(图3：B)，信号转导机制和一般功能预测的转录本最多，其次是翻译后修饰、蛋白转换、分子伴侣和转录。

KEGG注释结果显示有65 895条(24.42%)转录本被分配到386个KEGG通路上，其中与人类疾病相关的转录本(26 761条)最多，其次是生物系统(22 997条)(图3：C)。在二级分类中，参与信号转导的转录本(10 039条)和传染病：病毒(10 003条)的转录本数量最多。

2.4 转录本的组织特异性表达

转录本TSI值的频数分布图显示(图4：A)，仅17条转录本TSI值<0.15，它们在所有6个组织中表达差异很小，为低组织特异性基因，其中可能包括了一些重要的管家基因；而39.65%(66 763条)的转录本TSI值在0.15～0.85之间，为中等组织特异性基因，它们的表达量在各组织中存在差异，TSI值越大差异越大；60.34%的转录本TSI值>0.85，为高组织特异性基因，其在部分组织中的表达与其他组织存在明显差异，其中24.5%(41 354条)的转录本TSI值达到1.0，它们严格地在特定的某一种组织中表达。由表达量热图(图4：B)可见，这3类转录本在各组织中的表达模式不同。

2.5 各组织高表达转录本

为了解这些组织特异表达基因是否为各组织高表达基因，统计了6个组织中表达量最高的10条转录本(表3)，共包括39条转录本，其TSI值为0.48～1.0(平均0.83)。其中骨骼肌中10条转录本FPKM值均超过10 000，而脑组织中9条转录本FPKM值为1 213.09～3 074.01，仅1条>10 000，显示不同组织中基因表达量差异较大。这些转录本中，6条转录本(T_90962_c0_g2_i7、T_92170_c4_g8_i1、T_92170_c4_g7_i1、T_90443_c2_g7_i2、T_92619_c3_g1_i2和T_87893_c2_g3_i5)在骨骼肌和心脏中高表达，4条转录本(T_90962_c0_g2_i7、T_84299_c4_g5_i2、T_92588_c6_g2_i6和T_92288_c9_g4_i1)在肾脏和肝脏中高表达。39条转录本中26条(61.5%)为高组织特异性表达，包括8条在骨骼肌中高表达，7条在肺中高表达，6条在肝脏中高表达，肾脏、心脏、大脑中高表达的分别有5条、4条、3条；其他13个高表达转录本的组织特异性较低，如COX1在除肺以外的5个组织中都高表达。

表3 豹猫6个组织中表达量最高的10条转录本的NR注释结果及其组织特异性Table 3 NR annotation of the top 10 transcripts with the highest FPKM value in the 6 tissues of Prionailurus bengalensis and TSI analysis

续表3

3 讨论

作为一种分布广泛的小型猫科动物，过去关于豹猫的研究大多关注其种群分布、遗传多样性水平和系统发育关系等(Mukherjeeetal.，2010；Pateletal.，2017)。Ito等(2020)基于基因组SNP分析了豹猫遗传多样性，Bredemeyer等(2021)报道了豹猫的全基因组，但迄今为止关于豹猫基因表达调控相关研究仍非常缺乏。本研究采用RNA-seq对豹猫6个组织的转录组进行了测序、从头组装、功能注释及不同组织的基因表达谱分析，从而为丰富豹猫的遗传资源、辅助基因组注释、开发分子标记等提供了重要数据，也对于理解豹猫的基因表达调控，开展保护遗传学研究奠定了基础。

本研究基于豹猫6个组织的大量转录本从头组装了豹猫转录组，能更全面地体现豹猫基因表达情况，组装质量良好：N50长度为2 660 bp，42.24%的转录本长度超过1 000 bp。研究显示，越长的转录本序列越完整，组装质量越高，更容易获得较多关于基因的信息(Mengetal.，2015)。与发表的东北虎Pantheratigrisaltaica肺转录组(319 bp)和6个组织合并转录本(332 bp)(Luetal.，2016)、小尾寒羊和杜泊羊骨骼肌转录组(735 bp和706 bp)(张春兰，2016)、梭子蟹Portunustrituberculatus雌性(954 bp)和雄性(823 bp)(Wangetal.，2018)转录组序列平均长度相比，豹猫转录组质量更高。GO、KOG和KEGG注释结果显示出豹猫各组织表达不同生理功能相关的多种基因。仅有42.44%(114 517条)的豹猫转录本具有注释信息，超过50%的转录本无法比对到现有数据库，这可能是由于豹猫的基因组和蛋白质相关信息仍不完善，大量转录本可能是豹猫特有的转录本，深入研究这些新转录本可揭示豹猫独特的生理、遗传特征。

豹猫6个组织的TSI分布显示，TSI<0.15的管家基因最少，仅17条(0.01%)，60.34%的转录本TSI>0.85，其表达具有高度组织特异性。这显示尽管共有同一个基因组，但豹猫在不同组织中的基因转录和表达却存在很大差异，以保障不同组织执行不同的生理功能。Su等(2002)统计老鼠45个正常组织、Yanni等(2005)统计人类Homosapiens12种组织的TSI也显示，TSI值为0.9～1.0的高组织特异性表达基因最多，这与本研究结果一致。但不同的是，人类中57%、老鼠中56%的基因都属于中等组织特异性表达，即0.15≤TSI≤0.85。这可能是由于转录本注释质量的差异。人和老鼠具有高质量基因组，转录本注释完善，而豹猫缺乏相关的基因组和蛋白质信息，本研究从头组装的转录本高达26万余条，许多转录本并非全长基因，存在大量难以注释的短序列。

为进一步了解这些组织特异性表达的基因是否也是高表达的基因，统计了豹猫每个组织中表达量最高的前10条转录本。共计39个高表达的基因中，26个都属于在特定组织高表达的基因，这些基因都是与该组织特定功能密切相关的重要基因。如肝脏中特异高表达的RBP4与肝脏的脂肪降解功能相关，人的多组织基因表达图谱也显示该基因在肝脏和脂肪中特异表达(Fagerbergetal.，2014)。肺中特异高表达的SFTPC编码肺表面活性物质蛋白C，它通过降低覆盖肺部的液体的表面张力来维持肺组织的稳定性，SFTPC突变能引发婴幼儿间质性肺病，其过表达可抑制肺癌细胞的增殖(Hayasakaetal.，2018；Lietal.，2019)。大脑中特异高表达的PCP4与神经元的功能相关，小鼠胚胎发生过程中PCP4的过表达诱导神经元提前分化成熟，而在成年期过表达会导致小鼠学习障碍(Mouton-Ligeretal.，2014；Reneltetal.，2014)。此外我们还鉴定了不同组织共同高表达的基因。如COX1编码细胞色素c氧化酶Ⅰ，是线粒体呼吸链的末端组分，该酶可将电子从还原的细胞色素c转移到分子氧，并有助于维持线粒体内膜的电化学梯度(Baklouti-Gargourietal.，2013)，为ATP合成酶合成ATP提供电化学梯度(Dennerlein & Rehling，2015)，是细胞呼吸过程中重要的功能基因。COX1在心脏、肝脏、肾脏、大脑和骨骼肌中都是表达量前10的基因，在肺中的表达量也相对较高，证明该基因在各个组织中均发挥重要生理功能。TMSB4是一种广泛分布于真核细胞中的多功能肽，与细胞的移动性、分化能力及细胞分裂密切相关(Kuzan，2016)。TMSB4是肺、肾脏和大脑组织中共同高表达的基因，Du等(2015)发现其在大熊猫血液中也是高表达的基因之一，它可能也是维持细胞基本生理功能的重要基因。心脏和骨骼肌都属于肌肉组织，具有收缩性，本研究还鉴定到6条仅在心脏和骨骼肌中高表达转录本，这些转录本可能和肌肉组织的收缩功能密切相关。之前人类不同组织的基因表达谱研究也显示心脏和骨骼肌的表达谱相似性高(Shmuelietal.，2003)。此外肾脏和肝脏组织共同高表达的转录本有4条，这可能源于肝脏和肾脏之间有密切相关的生理功能。