光倒刺鲃（Spinibarbus hollandi）属鲤形目（Cypriniformes）鲤科（Cyprinidae）鲃亚科（Barbinae）倒刺鲃属（Spinibarbus），广泛分布于长江以南各水系，是重要的经济鱼种之一，其抗病性较强，肉质鲜美，营养价值较高[1-2]。近年来，有关光倒刺鲃的研究大多集中在养殖技术[3]、性别标记[4]、功能基因[5-6]和线粒体基因组[7] 等方面。
转录组测序是一种能全面快速获得细胞或组织在某种状态下所有转录本的序列信息和基因表达信息的技术[8]。通过分析转录组数据，可以发掘未知和稀有基因，从而准确研究基因表达差异、基因组成变异、分子标记筛选等问题[9-12]。目前，转录组测序技术因其数量大、覆盖率高、成本低、操作方便等特点，在鱼类研究中广泛应用。Hao等[13]对草鱼进行转录组测序分析，评估了草鱼肌肉的肉质和肌肉品质相关基因的表达情况，发现氨基酸代谢途径与鱼肉品质密切相关；石立冬等[14]对红鳍东方鲀肝脏转录组进行分析，筛选出与牛磺酸抗应激功能有关的基因；邵嘉棋等[15]对大口黑鲈进行RNA-seq技术测序，获得64个与驯食性状相关的候选基因和1个单核苷酸多态性（Single nucleotide polymorphism，SNP）标记，为其食性驯化遗传改良提供参考。
目前，光倒刺鲃基因组研究较少，本研究对该鱼类进行转录组测序与分析，以构建光倒刺鲃转录组文库，探讨其基因组结构，为解决光倒刺鲃基因进化、遗传育种和资源恢复等问题提供参考，对揭示其生长机制、抗逆能力，选育优势养殖品种以及种质资源保护具有重大意义。
1 材料与方法
1.1 试验材料
试验所用光倒刺鲃样品采自广东省珠江水产研究所。选取约200 g的雌雄个体各3尾，刮去表面鳞片，剪取各尾背部肌肉约20 mg。
1.2 RNA提取、文库构建与测序
利用Trizol试剂盒提取光倒刺鲃肌肉组织的RNA，采用Thermo Nanodrop 2000c、Agilent 2100分别测定RNA纯度、浓度及完整性。用含有Oligo（dT）的磁珠富集mRNA，然后加入裂解缓冲剂将mRNA消化成小片段，再以片段化的mRNA为模板合成cDNA，合成的cDNA用试剂盒进行纯化，在cDNA 3’末端加polyA，并连接测序接头，利用Hieff NGS® DNA Selection Beads磁珠进行纯化目标片段，然后进行PCR文库扩增，经质检后用Illumina HiSeq™ 2000测序仪对合格文库进行测序。
1.3 原始数据的组装拼接
对原始序列进行处理，去除接头序列和低质量序列，以获得Clean Data。然后对Clean Data进行统计，以评估数据的成分和质量，包括产量、Q20、N含量、GC含量等。最后用组装软件Trinity [16]将短reads组装成Contig和Unigene，并统计组装质量。
1.4 转录组注释、分析
1.5 微卫星查找和SNP检测
利用MicroSAtellite软件对拼接得到的Unigene进行微卫星定位。在MISA结果的基础上，保留前后序列均不小于150 bp的微卫星，对微卫星进行分型，用Primer primier 6.0软件的默认设置设计引物。最后使用SOAPsnp检测样品的SNP。
2 结果与分析
2.1 高通量测序产量统计
测序结果表明，总计产出13 525 001 683 nt数据，过滤后得到的clean reads数目为95 231 190个，Q20含量为96.40%，N含量为0，GC含量为49.43%。拼接后得到91 467个Unigene，总长为73 474 804 nt，Unigene长度分布如图1所示，其中序列长度均大于200 nt，且大于3 000 nt的有4 163个，最大、最小长度分别为40 474 和201 nt，平均长度为803 nt，N50达1 568个。
2.2 Unigene功能注释
利用Blast工具将91 467 个Unigene序列与NR、KEGG、KOG和GO数据库（E-value<10-5）进行比对，NR、GO、KOG、KEGG数据库分别注释了48 912、17 577、28 503、25 927个Unigene，共有49 061（53.64%）个Unigene注释到蛋白数据库。
2.2.1 NR注释
2.2.2 GO功能注释及分类
与GO数据库进行比对注释，共有17 557个Unigene分别注释到生物过程、细胞组分和分子功能3个主要生物类别的49个功能亚类中（图3）。在生物过程中，最多的是与细胞过程相关的基因，共8 793个，其次是与单生物过程相关基因，共有8 508个；在细胞组分中，与细胞和细胞组分有关的基因均有6 071个；在分子功能中，主要是与结合和催化活性功能相关的基因，分别有8 880个和5 758个。
2.2.3 KOG功能分类
将Unigene和KOG数据库进行比对（图4），结果显示，28 503个Unigene被分为25个KOG聚类。其中，排在前三的聚簇有信号转导机制，一般功能预测和翻译后修饰、蛋白质转化、分子伴侣，分别有14 872，10 527和4 846个Unigene，而Unigene数量最少的为细胞运动，仅181个。
2.2.4 KEGG代谢通路分析
2.3 微卫星分析
由表1可知，使用MISA软件共发现9 036条微卫星序列，含微卫星的Unigenes序列共有7 283条，占总Unigene（91 467条）的7.96%，其中，复合型的微卫星序列共有835条，其余为以不间断重复方式组成的单一型微卫星。由表2可知，在光倒刺鲃转录组中，二核苷酸重复序列最多，有5 204条，其次是三核苷酸重复序列和四核苷酸重复序列，分别有2 523和965条，五、六核苷酸重复序列共344条；此外，不同微卫星标记位点的重复次数也有差异，除了单碱基微卫星位点，最丰富的重复次数为6次，有2 718个位点，占总数的30.08%，其次是重复次数为5次，有1 605个位点，占总数的17.76%。由图6可知，在二核苷酸重复序列中，所占比例最高的是AC/GT，占总数的36.8%，其次是AG/CT和AT/AT，分别为14.6%和5.9%；在三核苷酸重复中，所占比例最高的是ATC/CTG，占总数的5.6%；四核苷酸重复序列中，比例最高的重复单元为AAAT/ATTT，占总数的3%；五、六核苷酸重复序列中，各重复序列所占比例均较低。
表1 微卫星位点检测结果 |
| 类型 | 数量/个 |
|---|---|
| Unigene数 | 91 467 |
| 发现微卫星位点数 | 9 036 |
| 包含微卫星位点的Unigene数 | 7 283 |
| 含有1个以上微卫星位点的Unigene数 | 1 240 |
| 复合型的微卫星序列 | 835 |
表2 微卫星不同重复单元的分布 |
| 重复类型 | 二核苷酸 | 三核苷酸 | 四核苷酸 | 五核苷酸 | 六核苷酸 | 总计 | 比例/% |
|---|---|---|---|---|---|---|---|
| 4 | 0 | 0 | 620 | 150 | 96 | 866 | 9.58 |
| 5 | 0 | 1 356 | 187 | 42 | 20 | 1 605 | 17.76 |
| 6 | 2 039 | 617 | 46 | 15 | 1 | 2 718 | 30.08 |
| 7 | 1 025 | 290 | 43 | 3 | 0 | 1 361 | 15.06 |
| 8 | 616 | 75 | 14 | 3 | 2 | 710 | 7.86 |
| 9 | 442 | 27 | 10 | 1 | 0 | 480 | 5.31 |
| ≥10 | 1082 | 158 | 45 | 10 | 1 | 1 296 | 14.35 |
| 总计 | 5 204 | 2 523 | 965 | 224 | 120 | 9 036 | 100 |
3 结论与讨论
高通量测序是一种高质量、低成本、高效率的转录组测序分析法[20]。本研究对光倒刺鲃转录组基因进行测序、组装和拼接，再通过生物信息学方法分析得到光倒刺鲃Unigene序列，对其进行注释、功能分类和微卫星分子标记开发，结果表明，拼接得到91 467个Unigene，基因序列平均长度为803 nt。Franssen等[21]认为，较长的测序读长可有效减少拼接错误，提高拼接重叠群长度，对于无参考基因组的物种来说具有重要意义。91 467个Unigene共有49 061（53.64%）个Unigene可注释到蛋白数据库，有48 912个Unigene可与NR数据库匹配；未与任何已知蛋白匹配的Unigene，推断其可能为新的基因。在这些匹配的Unigene中，与犀角金线鲃的Unigene匹配率最高，达到26%，其次是安水金线鲃和金线鲃，分别为23%和18%，说明在数据库中，光倒刺鲃的基因序列信息较少，而本次获得的光倒刺鲃Unigene可以丰富数据库中光倒刺鲃的基因信息。
53.64%的Unigene被注释到KOG、GO和KEGG数据库，接近50%基因未搜索到同源基因序列，说明目前基因数据库中收录的鱼类基因丰富度有待进一步提高。这3个数据库的注释结果有助于进一步了解光倒刺鲃基因的分子功能、所处细胞位置、参与生物过程、代谢途径和信号通路等，为光倒刺鲃的功能基因、相关基因的信号通路、基因所处的位置及生理功能等的开发利用提供基础信息。其中一些与细胞免疫功能相关的基因可以用来制作基因表达芯片，用于光倒刺鲃免疫水平检测，还可作为光倒刺鲃抗病品系选育的生化指标。
在光倒刺鲃转录组中，二核苷酸重复序列的微卫星含量最高，其次是三核苷酸重复序列，这与大部分真核生物的研究结果一致；四、五、六核苷酸重复序列的微卫星丰富度明显低于二、三核苷酸，该结果也在其他真核生物中得到了验证[22-23]。在光倒刺鲃转录组的二核苷酸重复序列中，以AC/GT的数量最多，占36.8%，该结果与绝大多数脊椎动物的基因组序列研究结果相一致[24-26]。其次是AG/CT和AT/AT，而CG/CG的含量最少，极少的CG微卫星重复单元数量可能是因为在基因组中胞嘧啶C容易被甲基化脱氨基转化为胸腺嘧啶T。染色体间有频繁的滑链错配，具有较高的突变率。较高的突变率与高多态性相关联，序列越长，其在转录组中的多态率越高。本研究中，大多数的微卫星序列属于单一型，与大部分真核生物通过转录组检测到的微卫星序列结果相一致[27-28]。可能是因为转录组研究的是在RNA水平上的基因表达情况，其主要位于CDS区，基因序列的保守度高，变异小，因此微卫星重复数较少。
综上，本研究对光倒刺鲃转录组进行高通量测序并对其结果进行分析，在光倒刺鲃转录组的高通量测序中，共得到91 467个Unigene，其中49 061个Unigene注释到蛋白数据库，NR、GO、KOG、KEGG数据库分别注释了48 912、17 577、28 503、25 927个Unigene；在光倒刺鲃转录组中，二核苷酸重复序列的微卫星含量最高，有5 204条，其次是三核苷酸重复序列，有2 523条。研究结果为光倒刺鲃的功能基因、相关基因的信号通路及生理功能等的开发利用提供参考。