大口黑鲈(Micropterus salmoides)俗名加州鲈，原产于北美洲的密西西比河水系，属于鲈形目、太阳鱼科、黑鲈属，具有生长迅速、抗病力强、肉质鲜美、耐低温和易捕捞等特点，是重要的淡水名特优养殖肉食性鱼类之一。从20世纪80年代初引种到广东省，现在全国大部分水产养殖地区均有养殖。但20多年来，由于种质的提纯复壮工作严重滞后，使得大口黑鲈养殖群体遗传多样性降低，导致出现种质严重退化的现象，降低了生产性能。

高通量RNA测序技术即RNA-seq技术，能够在没有该物种完整基因组序列的前提下，准确获取特定组织在某一状态下转录组基因的全面信息，完整记录该条件下基因表达水平、生物学过程的分子机制以及遗传标记信息等^[1-7]。当前大口黑鲈基因组序列未知，遗传背景尚不清楚，限制了对大口黑鲈相关功能基因(如与生长、繁殖和抗病等重要经济性状相关的基因)的挖掘；也造成其基因组学研究的困难。因此，丰富大口黑鲈转录组信息具有重要的科学研究价值。

鉴于此，本研究利用RNA-Seq测序技术开展大口黑鲈组织的高通量转录组测序，通过对获得的转录本数据进行拼接与组装、基因功能注释以及微卫星 (SSR)位点挖掘等研究，获得更多的大口黑鲈转录本和更为全面的转录组信息，以期为大口黑鲈重要功能基因的发掘、相关基因克隆与表达和品种改良等问题奠定理论基础。

1.   材料与方法

1.1   材料

实验材料为体表无伤健康的大口黑鲈3尾，体质量500~700 g。均采集自中国水产科学研究院珠江水产研究所养殖基地。

1.2   样品准备

取样于2016年5月13日上午10时进行，分别取大口黑鲈脾脏、肌肉、肾脏、肝脏、心脏、胃和性腺组织，取样规格100 mg，于– 80 ℃超低温冰箱中保存备用。

1.3   总RNA提取与检测

将上述大口黑鲈组织样品分别进行总RNA抽提，按照Invitrogen公司的Trizol Reagent试剂说明书方法进行操作。样品符合转录组RNA检测标准：RNA总量≥10 μg，OD_260/280为1.8~2.1，28S∶18S≥1.5∶1.0，RIN (RNA integrity number)≥8.0。通过1.2%琼脂糖凝胶电泳和Nanodrop-2000核酸蛋白测定仪检测总RNA的完整性和纯度质量，以保证使用合格的样品进行转录组测序。随后将检测合格的RNA组织样品各取10 μL混合为一个样品，准备后续测序分析。

1.4   转录组测序与分析

利用Oligo (dT)磁珠和磁分离器分离并纯化出mRNA，将mRNA打断为片段，通过RT-PCR扩增富集得到大口黑鲈转录组cDNA文库，最后将构建好的文库进行高通量测序。测序部分由联川生物科技有限公司(杭州)协助完成，并对所获得的序列进行拼接组装及SSR分析。最后将组装得到的大口黑鲈转录本数据与6个公共数据库Swiss-Prot (Swiss Prot protein database)、Nr (non-redundant protein sequences)、KEGG (Kyoto Encyclopedia of Genes and Genomes)、KOG (Eukaryotic Ortholog Groups)、Pfam (protein families database)和GO (Gene Ontology)进行Blast比对。选择阀值条件为E value<1e−10，通过序列相似性进行功能注释。

2.   结果

2.1   RNA检测

大口黑鲈混合组织的RNA样本具有完整清晰的28S、18S和5S带型，OD_260/280值为2.04，28S∶18S值为1.8，RIN值为8.5，说明所得到的RNA质量较好，符合后续转录组测序要求。

2.2   测序质量控制

采用Illumina Hiseq 2500测序平台对大口黑鲈组织进行转录组测序。对获得的测序数据进行统计，共产生42 415 596条原始读序(reads)，总碱基数为6.36 G。序列质量评估显示，碱基Q20 (质量不低于20的碱基)为95.33%，Q30 (质量不低于30的碱基)为90.47%，GC含量平均值为48.86%，其中GC含量过高(大于80%)或过低(小于20%)的unigenes不存在，GC含量基本呈正态布(图1)，表明测序碱基组成情况良好，测序结果质量较高，能用于后续数据分析。

图  1  Unigenes GC含量分布

Figure  1.  GC content distribution of unigenes

下载: 全尺寸图片 幻灯片

2.3   测序数据组装与分析

经过De novo拼接后，获得42 264条transcripts，序列信息达到36 549 208 bp，其序列长度中位数(N50)的片段长度达到了1 331 bp，平均长度为864 bp。在transcripts数据基础上，经trinity软件对序列进行组装，最终获得35 659条unigenes，序列信息达到26 340 493 bp，序列大小201~35 516 bp，平均长度为738 bp，N50为1 052 bp (表1)。其中长度200~500 bp的transcripts有20 207条，500~1 000 bp的11 114条，1 000~2 000 bp的7 110条，≥2 000 bp的3 833条；200~500 bp的unigenes有18 863条，500~1 000 bp的9 445条，1 000~2 000 bp的5 141条，≥2 000 bp的2 210条(图2)，以上这些测序结果说明数据组装质量较高。通过GENScan软件预测unigenes编码区，有17 270条序列可被编码，占全部unigenes的48.43%。

表  1  转录组数据拼接结果统计

Table  1.  Assembly result of transcriptome data

序列种类 category	总数 total number	N50长度/bp N50 length	总长度/bp total length	最大长度/bp maximum length	中等长度/bp median length	最小长度/bp minimum length	平均长度/bp average length
transcript	42 264	1 331	36 549 208	35 516	522	201	864
unigene	35 659	1 052	26 340 493	35 516	472	201	738
注：N50表示将转录本从长到短排序，依次累加碱基数，当累计碱基数达到转录本总碱基数的50%时的转录本的长度 　Note: N50 of transcript or unigenes was calculated by ordering all sequences, then adding the lengths from longest to shortest until the summed length exceeded 50% of the total length of all sequences.

下载: 导出CSV 

| 显示表格

图  2  转录组transcript和unigenes长度分布

A. transcript的长度分布；B. unigenes的长度分布

Figure  2.  Length distribution of assembled transcripts and unigenes of transcriptome

A. size distribution of transcripts; B. size distribution of unigenes

下载: 全尺寸图片 幻灯片

2.4   微卫星位点分析

利用MISA软件对35 659条unigenes进行SSR位点搜索，共检测到含有SSR位点的序列294条，共鉴定4 030个SSR位点。SSR类型较为丰富，单核苷酸至六核苷酸重复类型均有检测到。其中出现频率最高的SSR类型为二核苷酸(2 131个)；其次为三核苷酸(1 009个)；单核苷酸SSR位点有764个；四和五核苷酸SSR位点的数量相对较少，分别有77和30个；六核苷酸SSR位点的数量最少(19个)。其中单核苷酸SSR的重复类型主要为A和T；二核苷酸重复的主要类型为GT和AC；三核苷酸重复的主要类型是TGC、GCA、GGA、TAA、GAG、CTC、TCC和AGG；四核苷酸重复的主要类型为AAGA。这些SSR位点的鉴定，可为开展大口黑鲈差异基因分析、分子标记开发和遗传图谱构建提供良好基础，并为后续大口黑鲈生物标记开发与应用提供分子水平的参考依据。

2.5   基因功能及Nr数据库注释

对组装得到的35 659条unigenes与已知的KOG、Nr、Pfam、Swiss-Prot、GO和KEGG数据库进行注释。共有5 617条unigenes在以上6个数据库中同时注释成功，占总数的15.75%，至少在1个数据库中注释成功的有17 253条，占总数的48.38%。其中15 024条unigenes在GO数据库中获得注释，占42.13%；15 832条unigenes在KOG数据库中获得注释，占44.40%；11 185条unigenes在KEGG数据库中获得注释，占31.37%；14 524条unigenes在Pfam数据库中获得注释，占40.73%；27 515条unigenes在Swiss-prot数据库中获得注释，占47.60%；在Nr数据库中，共有21 279条unigenes得到注释，占总数的59.67%。通过与Nr库进行比对注释，大口黑鲈转录组与深裂眶锯雀鲷(Stegastes partitus)能匹配到的相似基因序列最多，达37.4%；其次为尼罗罗非鱼(Oreochromis niloticus，12.1%)；相似性序列数量大于4%的物种有斑马拟丽鱼(Maylandia zebra，6.5%)、布氏新亮丽鲷(Neolamprologus brichardi，4.8%)、伯氏朴丽鱼(Haplochromis burtoni，4.3%)和红丽鱼(Pundamilia nyererei，4.2%)；而其他物种占了30.7%，这可能与数据库中没有大口黑鲈的全基因组序列有关。

2.6   GO功能注释

根据GO数据库注释分类信息，一共有15 024条unigenes得到注释，占42.13%。GO功能大致可分为基因参与的生物过程、所处的细胞组分和具有的分子功能3大类(图3)。这3大类别又被详细划分为50个亚类，例如代谢过程、生化过程、免疫防御生长与发育过程等。其中生物过程包含25个不同的亚类，也是3大类别中所含类别最多的一类，其中转录、DNA依赖性所占比例最高(71.13%)，其次为转录(60.47%)，而内噬作用所占比例最低，仅有56条(10.03%)；在细胞组分类别中有15个亚类，细胞核和细胞浆所占比例最高，分别为81.34%和72.48%，而高尔基体最低，仅有31条(8.65%)；在分子功能类别中有10个亚类，ATP结合和锌离子结合所占比例最高，分别为60.56%和51.45%，而蛋白质丝氨酸特异酶结合所占比例最低，仅有21条序列(7.79%)。

图  3  Unigenes GO功能注释

1. 转录、DNA依赖性；2. 转录调控、DNA依赖性；3. 蛋白质转运；4. 凋亡；5. 细胞分化；6. 细胞黏附；7. 蛋白质水解；8. 多细胞器官发育；9. 细胞周期；10. 细胞分化；11. 胞内信号转导；12. mRNA加工；13. 信号转导；14. 小G蛋白介导信号转导；15. 转运；16. 转录正调控；17. 有丝分裂；18. 染色质修饰；19. 转录负调控；20. RNA剪切；21. 胞内蛋白质运输；22. 翻译；23. DNA修复；24. 转录；25. 内噬作用；26. 细胞核；27. 细胞浆；28. 必须膜；29. 胞液；30. 细胞质膜；31. 内质网膜；32. 核仁；33. 线粒体；34. 胞外区；35. 核质；36. 隔膜；37. 高尔基体膜；38. 细胞骨架；39. 高尔基体；40. 细胞核周区；41. ATP结合；42. 锌离子结合；43. 蛋白结合；44. DNA结合；45. 金属离子结合；46. RNA结合；47. 钙离子结合；48. 特异序列DNA结合；49. 结合；50. 蛋白质丝氨酸特异酶结合

Figure  3.  GO functional annotation of unigenes

1. transcription, DNA-dependent; 2. regulation of transcription, DNA-dependent; 3. protein transport; 4. apoptosis; 5. cell division; 6. cell adhesion; 7. proteolysis; 8. multicellular organismal development; 9. cell cycle; 10. cell differentiation; 11. intracellular signal transduction; 12. mRNA procession; 13. signal transduction; 14. small GTPase mediated signal transduction; 15. transport; 16. negative regulation of transcription; 17. mitosis; 18. chromatin modification; 19. positive regulation of transcription; 20. RNA splicing; 21. intracellular protein; 22. translation; 23. DNA repair; 24. transcription; 25. endocytosis; 26. nucleus; 27. cytoplasm; 28. integral to membrane; 29. cytosol; 30. plasma membrane; 31. endoplasmic reticulum membrane; 32. nuclelous; 33. mitochondrion; 34. extracellular region; 35. nucleoplasm; 36. membrane; 37. golgi membrane; 38. cytoskeleton; 39. golgi apparatus; 40. perinuclear region of cytoplasm; 41. ATP binding; 42. zinc ion binding; 43. protein binding; 44. DNA binding; 45. metal ion binding; 46. RNA binding; 47. calcium ion binding; 48. specific DNA sequence binding; 49. binding; 50. protein serine/threonine-specific kinase

下载: 全尺寸图片 幻灯片

2.7   KOG功能注释

将在KOG数据库获得注释的13 466条unigenes进行功能分类，共获得25个功能组分，涉及了大多数的生命活动(图4)。其中，信号转导机制类所占比例最高，为17.85%(2 404条)，其后依次是一般功能预测类(17.35%，2 337条)、翻译后修饰－蛋白质周转－分子伴侣(10.02%，1 396条)、转录类(8.81%，1 187条)、功能未知类(6.73%，906条)、细胞内分泌和囊泡运输(6.61%，891条)、细胞骨架(5.94%，801条)和RNA加工和修饰(4.37%，589条)，而参与辅酶转运和代谢(0.63%，85条)及细胞运动(0.38%，51条)的unigenes最少。

图  4  Unigenes的KOG注释

A. RNA加工和修饰；B. 染色体结构和动力学；C. 能源生产和转换；D. 细胞周期调控－细胞分裂－染色体分离；E. 氨基酸转运和代谢；F. 氨基酸转运和代谢；G. 碳水化合物转运和代谢；H. 辅酶转运和代谢；I. 脂质转运和代谢；J. 翻译－核糖体结构－生物合成；K. 转录；L. 复制－重组－修复；M. 细胞壁－细胞膜合成；N. 细胞运动；O. 翻译后修饰－蛋白质周转－分子伴侣；P. 无机离子转运与代谢；Q. 次生代谢产物生物合成、转运和分解代谢；R. 一般功能预测；S. 功能未知；T. 信号转导机制；U. 细胞内分泌和囊泡运输；V. 防御机制；W. 胞外结构；Y. 核结构；Z. 细胞骨架

Figure  4.  KOG annotation of unigenes

A. RNA processing and modification; B. chromatin structure and dynamics; C. energy production and conversion; D. cell cycle control, cell division, chromosome partitioning; E. amino acid transport and metabolism; F. nucleotide transport and metabolism; G. carbohydrate transport and metabolism; H. coenzyme transport and metabolism; I. lipid transport and metabolism; J. translation, ribosomal structure and biogenesis; K. transcription; L. replication, recombination and repair; M. cell wall/membrane/envelope biogenesis; N. cell motility; O. posttranslational modification, protein turnover, chaperones; P. inorganic ion transport and metabolism; Q. secondary metabolites biosynthesis, transport and catabolism; R. general function prediction only; S. function unknown; T. signal transduction mechanisms; U. intracellular trafficking, secretion, and vesicular transport; V. defense mechanisms; W. extracellular structure; Y. nuclear structure; Z. cytoskeleton

下载: 全尺寸图片 幻灯片

2.8   KEGG功能分类

采用KEGG数据库对大口黑鲈可能参与的生理生化反应途径进行注释，有11 185条(31.37%) unigenes在KEGG数据库中得到注释。这些unigenes共参与了细胞过程、环境信息处理、遗传信息处理、代谢和有机系统5大类共29小类的267个功能通路(图5)。参与信号转导通路的unigenes较多，共有1 349条，其中较多的几个分别是MAPK信号通路(471条)、Calcium信号通路(241条)和ErbB信号通路(161条)，这几个代谢通路与环境信息处理大类中信号转导相关；其次是免疫系统通路(1 002条)；参与萜类和酮类化合物代谢的unigenes数量最少，仅有49条。

图  5  Unigenes的KEGG注释

Figure  5.  KEGG annotation of unigenes

下载: 全尺寸图片 幻灯片

3.   讨论

近年来，高通量测序技术在揭示转录组差异基因的表达、功能基因的发掘、分子标记开发、相关标记的分析和非编码RNA鉴定与调控等方面起着重要作用，已广泛应用于动植物的转录组和基因发掘与功能研究中^[8-12]。这种高通量测序技术在一些鱼类的转录组研究中也有报道，如虹鳟(Oncorhynchus mykiss)^[13]，欧洲鳗鲡(Anguilla anguilla)^[14]，大菱鲆(Scophthalmus maximuz)^[15]、达氏鳇(Huso dauricus)^[16]、大黄鱼(Larimichthys crocea)^[17]、斑点叉尾鮰(Ictalurus punctatus)^[18]、大西洋鲑(Salmo salar)^[19]和草鱼(Ctenopharyngodon idella)^[20]等。通过这些研究获得了大量的鱼类转录组信息，如赵刚等^[21]利于该技术对岩原鲤(Pocypris rabaudi)全组织转录组进行了测序，结果获得了83 252条unigenes，平均长度为787 bp，有37 157条unigenes得到注释；许建等^[22]利用RNA-Seq技术对鲮(Cirrhina molitorella)进行了转录组分析，经过拼接与组装，最终获得了19 962条unigenes，平均长度为1 269 bp，有10 577个特异蛋白通过基因功能得到注释，鉴定了具有完整开放阅读框 (ORF)的全长cDNA序列共5 931条。这些信息的获得能为后续水产动物功能基因发掘和转录调控研究提供有价值的数据参考。目前，有关大口黑鲈转录组的研究还很少，已有的文献中，仅见到对大口黑鲈肝脏和肌肉混合组织通过转录组测序进行SNPs筛选与生长关联分析的研究，而没有进一步对转录组数据进行功能注释信息分析^[23]。本研究利用拼接得到的unigenes序列在Pfam、Swiss-prot、KEGG、Nr、GO和KOG 6个公共数据库中进行注释，结果发现仅有5 617条unigenes在公共数据库中均能获得注释。造成这种结果的原因可能有3个方面：1)数据库中大口黑鲈基因序列信息少，尚无大口黑鲈基因组测序研究的报道；2)转录组测序也可能有一些新出现的基因是目前数据库里所没有的，因鱼种类繁多、种间差异巨大可能导致其同源序列难以被发现；3)测序技术的局限性。基于从头拼接测序技术的准确度和拼接长度极大依赖于测序深度和参考基因组序列，而组装测序中序列长度较短也会影响后期数据质量，如过短的unigenes序列在数据库中无法获得有效的比对信息。

Nr数据库物种注释中，有37.4%的unigenes在深裂眶锯雀鲷基因中获得注释，在所注释物种中比例最高，可能是由于这2个物种同属鲈形目，在亲缘关系上可能较近；而有30.7%的unigenes未被注释到其他物种中，以上结果也说明，数据库中已有的大口黑鲈序列信息还很少，难以将大口黑鲈本身特有的基因信息与其他物种的基因信息分离开来。本研究获得的unigenes进一步丰富了目前数据库中已有的鲈形目基因资源。KEGG代谢通路分析表明11 185条大口黑鲈unigenes共参与到267个代谢通路中。其中与信号转导通路相关的unigenes数量最多，主要有MAPK、Calcium和ErbB等信号通路，这些信号通路都与疾病发生、细胞离子交换和细胞凋亡相关。其次是参与免疫系统通路的unigenes，有1 002条。这些与免疫通路相关的基因包括抗病毒免疫因子，如炎症细胞因子和受体、补体成分、模式识别受体、B细胞、T细胞抗原活化分子和干扰素等。SSRs作为一种分子标记在动物中研究很广泛。利用MISA软件查找测序序列，共发现了4 030个SSR位点，而且发现SSR重复碱基以二核苷酸类型最多，有2 131个分子标记，占所有SSR位点的52.88%，这为进一步利用分子标记研究大口黑鲈功能基因定位、群体遗传多样性、遗传连锁图谱构建等奠定了基础。