基于转录组测序的方斑东风螺单核苷酸多态性位点挖掘及功能注释

2021-02-03刘付柏许尤厚黄宝松王忠良

广东海洋大学学报 2021年1期

王菁，刘付柏，许尤厚，黄宝松，王忠良

（1.广东海洋大学水产学院，广东湛江 524088；2.广西北部湾海洋生物多样性养护重点实验室，广西钦州 535000）

方斑东风螺 (Babylonia areolata) 俗称花螺，隶属于软体动物门腹足纲蛾螺目，有生长速度快、养殖周期短、肉味鲜美、软体部不饱和脂肪酸含量丰富、经济价值高、便于运输等优点[1-3]。随着方斑东风螺养殖的快速发展，养殖过程中生长缓慢、病害暴发等问题日益突出，须通过遗传改良、病害防控、饲料营养优化等解决，其中遗传改良对促进东风螺养殖产业健康持续发展有现实意义[4-5]。

单核苷酸多态性(Single nucleotide polymorphism,SNP)是指因单个核苷酸变异引起的DNA 序列多态性，是一种常见的基因突变[6-7]。SNP标记为第3 代分子标记技术，与限制性片段长度多态（restriction fragment length polymerphisms，RFLP）、微卫星多态(microsatellite polymorphisms)相比，在基因组中位点丰富，代表性强，遗传稳定，可自动检测，有低成本、高效率的优点[8-10]，因而广泛应用于水产动物研究[11-18]。

随着高通量测序技术的快速发展，转录组测序成本大幅下降，利用比较转录组学方法和序列比对识别大量的SNP位点已逐渐成为一种趋势。目前，已对栉孔扇贝(Chlamys farreri)[19]、马氏珠母贝(Pinctada fucata)[20]、波纹唇鱼(Cheilinus undulatus)[21]、曼氏无针乌贼(Sepiella japonica)[22]、棘头梅童鱼(Collichthys lucidus)[23]、大菱鲆(Scophthalmus maximus)[24]、大口黑鲈(Micropterus salmoides)[25]、凡纳滨对虾(Litopenaeus vannamei)[10,26]等SNP位点进行多态特征分析，但未见关于方斑东风螺SNP位点的研究。通过转录组测序技术，结合所测物种基因组信息，更易找到与目标性状相关的SNP[27-28]。本研究通过对方斑东风螺转录组的深度测序分析，筛选出大量SNP位点，并对这些SNP所在基因进行功能注释，为方斑东风螺的抗病及育种研究提供基础数据。

1 材料和方法

1.1 材料

方斑东风螺购自广东省湛江某东风螺养殖场，平均体质量为25 g，于实验室海水桶中暂养1 周（80 L，25℃）。实验时，将方斑东风螺分为脂多糖（LPS）注射组（LPS-4h、LPS-8h）和空白对照组。实验组对方斑东风螺闭壳肌注射100 μL 0.5 mg/mL 的LPS悬浮液，对照组注射同体积的磷酸盐缓冲液（PBS），分别于刺激后4 h、8 h（分别记为LPS-4h、LPS-8h组）采集各组足组织样品，于液氮中速冻，置-80 ℃下保存备用。

1.2 转录组测序数据

样品委托广州基迪奥生物科技有限公司使用TRIzol 试剂（Invitrogen，美国）提取总RNA，经纯度（NanoDrop 2000）和浓度（Agilent 2100）检测后，同处理组10 个个体样本合并为1 例样品进行建库测序，经过Illumina/ Hiseq-2000 高通量转录组测序，删除大量低质量原始数据和适配子等，使用组装程序Trinity 对转录组数据进行序列组装(Raw Reads 数据的SRA 登录号为SRP216586)。

1.3 SNP位点的检测

通过 Call snp 软件 bcftools (https://github.com/samtools/bcftools) 获得在不同处理组间有表达差异的SNPs 位点，使用SOAPsnp 对获得的SNP进行统计和分析。

1.4 SNP位点所在 unigene 的注释与功能分析

基于所得 SNP-unigenes 序列与 Nr、Nt 及Swiss-prot 数据库比对后的蛋白功能注释信息，对测序数据的KEGG 通路、差异表达基因和SNP进行分析。从转录组中筛选免疫相关的KEGG 通路，根据通路中的免疫基因筛选出差异表达的免疫防御相关基因，并进行SNP位点分析。

2 结果与分析

2.1 转录组测序及序列组装

用Illumina 高通量测序平台对方斑东风螺足组织进行转录组测序，对原始数据进行严格的质量控制，经过滤后在对照组转录组（BLANK）和实验组转录组（LPS-4h、LPS-8h 组）分别获得56 424 638、50 596 990、49 924 362 条纯净序列，其中GC 比例分别为47.8%、46.15%、46.88%，碱基质量值Q30分别为95.40%、95.61%、95.40%（表1）。说明方斑东风螺测序质量较高，转录组数据可用于后续分析。

对测序数据进行序列组装，共获得81 773 条Unigene，总长度为55 763 627 bp，平均长度为681 bp，N50 的长度为1 035 bp (表2)。转录本和Unigene 的N50 长度均远大于其平均长度，证明组装效果较佳。

2.2 转录组数据及SNP位点数据

利用在对照组转录组（BLANK）和实验组转录组（LPS-4h，LPS-8h）获得的数据，经SOAPsnp软件检测，从37 136、37 076、36 657 条unigenes中分别获得224 055、225 287、224 440 个SNP位点，所有SNP位点中，纯合SNP位点107 407 个（BLANK组35 635个，LPS-4h组35 543个，LPS-8h组36 229 个）（表3）。对于Unigene 上的SNP位点分析统计发现，BLANK 转录组中含1 个SNP位点的unigene 有9 802 条（26.39%），含2～10 个SNP位点的unigene 21 224 条（57.15%），含10 个以上SNP位点的unigene 6 110 条（16.45%）；LPS-4h 转录组中含有1 个SNP位点的unigene 有9 721 条（26.22%），含2～10 个SNP位点的unigene 21 224条（57.19%），含10 个以上SNP位点的unigene 6 110条（16.59%）；LPS-8h 转录组中含1 个SNP位点的unigene 有9 580 条（26.13%），含2～10 个SNP位点的unigene 20 921 条（57.07%），含10 个以上SNP位点的unigene 6 156 条（16.79%）（图1）。

表1 测序数据质量分析Table 1 Quality analysis of sequencing data

表2 单基因簇统计分析Table 2 Statistics analysis of Unigenes

表3 SNP位点数量概况Table 3 SNPidentified in BLANK,LPS-4h and LPS-8h transcriptomes

图2 表明，BLANK 转录组的纯合SNP位点中，颠换位点 82 447 个(36.80%)，转换位点141 608 个(63.20%)；LPS-4h 转录组的纯合SNP位点中，颠换位点82 878 个 (37.12%)，转换位点142 409 个 (63.78%)；LPS-8h 转录组的纯合SNP位点中，颠换位点82 444 个 (36.73%)，转换位点141 996 个 (63.27%)。在6 种核苷酸的变异类型中，以A/G 转换最多，分别占纯合SNP总数的31.83%（BLANK 组71 316 个）、31.85%（LPS-4h组71 706 个）和31.91%（LPS-8h 组71 630 个）。

图1 SNP的分布统计Fig.1 SNPdistribution in BLANK,LPS-4h and LPS-8h transcriptomes

图2 SNP类型分析Fig.2 SNPnumbers of different mutation types in BLANK,LPS-4h and LPS-8h transcriptomes

2.3 SNP-unigene 的注释与功能

COG 结果显示，共有16 891 条SNP-unigenes匹配相应的COG 注释信息；根据功能信息可分为26 类，其中“仅通用功能预测”和“信号转导机制”类最多，分别包含2 848 和2 814 条SNP-unigenes（图3A）。GO 分析表明，共有4 682 条SNP-unigenes 匹配到GO 条目(GO term)，GO 条目包含生物过程、细胞组分及分子功能的42 个亚类，分别在代谢过程、催化活性和细胞部分类中最为富集 (图3B)。KEGG富集分析显示，共有5 866 条SNP-unigenes 富集到298 个KEGG 子集中，其中以“内吞作用”子集中富集的SNP-Unigene 最多，共计182 条（图4）。

图3 SNP-unigenes COG 与GO 功能注释分析Fig.3 Cluster of orthologous groups (COG) and gene ontology (GO) classification of SNP-unigenes

图4 SNP-unigenes 的KEGG 信号通路富集分析Fig.4 Kyoto encyclopedia of genes and genomes (KEGG) classification of SNP-unigenes

2.4 免疫防御相关SNP-unigene的富集及特异SNP位点

根据KEGG 信号通路的富集分析，筛选到515个免疫防御相关 SNP-unigenes，注释到“Autophagy-animal”等19 条与免疫功能相关的信号通路中，其中以“Autophagy-animal”信号通路中注释的unigenes 最多(92条)，其次分别为“mTOR signaling pathway”“Wnt signaling pathway”“FoxO signaling pathway”等信号通路(表4)。

2.5 差异表达免疫防御相关基因SNP位点

根据转录组中unigene 的SNP位点分布情况及KEGG 功能注释信息筛选BLANK、LPS-4h、LPS-8h转录组中特异分布的免疫防御相关基因SNP位点，发现大量的SNP位点存在于免疫防御相关基因中（表5）。基于RPKM 标准化分析unigene 的表达水平，筛选阈值为q＜0.05，且 |log2（差异倍数）|＞1，获得大量差异表达基因。统计这些差异表达基因上的SNP位点，并对其中涉及免疫防御的基因进行KEGG 通路分析。注射LPS 4 h 后，差异表达的免疫防御相关基因主要参与胞吞作用、IL-17 信号通路、mTOR 信号通路、Wnt 信号通路。注射LPS 8 h 后，存在大量SNP位点的免疫相关基因主要参与吞噬体、溶酶体通路。比较注射4 h 与8 h，免疫相关基因则主要参与基座切除修复、谷胱甘肽代谢、药物代谢（细胞色素P450）、细胞色素代谢通路（表5）。

表4 免疫防御相关SNP-unigene 的KEGG 富集分析Table 4 KEGG enrichment analysis of immune-related SNP-unigenes

表5 BLANK、LPS-4h、LPS-8h 转录组特异分布免疫防御相关基因SNP位点分析(仅展示前10 行)Table 5 Specific immune-related SNPs identified in BLANK,LPS-4h,LPS-8h transcriptomes (Show only the first 10 rows)

3 讨论

在鱼类和贝类中，机体对病原体的抵抗力部分受基因控制[29]。开展SNPs 检测有助于了解不同个体和群体对外部刺激的反应[30]。因此，开发与方斑东风螺免疫相关的功能基因，挖掘与免疫性状连锁的分子标记是全面了解方斑东风螺免疫反应，开展分子标记辅助育种的基础。本研究基于Hiseq-2000测序技术分别从BLANK、LPS-4h、LPS-8h 转录组中获得224 055、225 287、224 440 个SNP位点。SNP碱基替换类型分为转换 (A-G、T-C) 和颠换(A-T、C-G、A-C、T-G) 两类。理论上，颠换与转换比例应为1∶2[31]。本研究中，转录组的颠换与转换SNP数比值约为1∶1.72，与先前研究中碱基转换类型发生频率高于颠换类型类似[31-33]；在6 种核苷酸变异类型中，以A/G、C/T 转换最多，约占6种核苷酸变异类型的60%，与其他物种的碱基替换频率类似[34-37]。

根据SNPs 所处位置，可将SNPs 分为基因编码区SNPs 和基因非编码区SNPs[38]，而位于基因调节区的SNP则称为调节SNPs，如其与基因相互作用影响到转录表达水平，则间接影响RNA 或蛋白质的表达数量，从而引起不同的机体反应[39]。本研究中，共有5 866 条SNP-unigenes 富集到298 个KEGG 子集中，以“内吞作用”子集中富集的SNP-Unigene 最多。进一步筛选到 515 个SNP-unigenes 注释到等19 条与免疫功能相关的信号通路，大部分SNP-unigenes 参与“自噬（动物）”“mTOR 信号通路”“Wnt 信号通路”“FoxO 信号通路”“细胞凋亡”等免疫通路。本研究表明，LPS注射4 h 后，大量差异表达的SNP-unigenes 主要参与胞吞作用、IL-17 信号通路、mTOR 信号通路/ Wnt信号通路。其中，参与胞吞作用的FGFR3基因上存在60 个SNP位点，FGFR 是酪氨酸激酶家族一员，它的激活可能导致肿瘤细胞生长和存活增加[40-41]。Wnt 信号通路是一个高度保守的系统，调控所有后生动物的复杂生物学过程。本研究发现，Wnt10基因中共存在26 个SNP位点。Wnt10可减少肿瘤细胞间互相黏连，促进肿瘤细胞间转移，进而调节干细胞行为、组织稳态和损伤修复[42]。LPS注射8 h 后，有大量SNP位点的免疫相关基因主要参与吞噬体、溶酶体通路。参与吞噬体通路的MPO上存在49 个SNP位点，参与溶酶体通路的LITAF基因上存在41 个位点。MPO 是一种血色素蛋白，是活化的中性粒细胞分泌的过氧化物酶类，可能通过各种炎症反应加速动脉粥样斑块的氧化，增强巨噬细胞吸收和泡沫细胞形成[43-44]。脂多糖诱导的肿瘤坏死因子 (Lipopolysac-charide-induced TNF-alpha factor，LITAF) 可调控TNF-α、IL-2 等细胞因子，在无脊椎动物的先天性免疫系统中扮演着重要的介质作用[45-46]。本研究LPS-4h 转录组与LPS-8h 组比较，差异表达的APEX1基因(33 个SNP位点) 和GST基因（32 个SNP位点）主要参与基座切除修复、谷胱甘肽代谢、药物代谢 (细胞色素P450)、细胞色素代谢通路。APEX1既有DNA 修复酶活性，又有氧化还原功能，在多种恶性肿瘤中参与肿瘤的侵袭与转移[47-48]。谷胱甘肽S-转移酶(GST) 作为抗氧化系统中第2 道防线主要成员，催化谷胱甘肽与亲电子外源化学物结合，将毒性物质转化成易排泄的水溶性形式[49]。推测方斑东风螺受刺激后，这些SNP-unigenes 可能广泛参与免疫反应[50-52]。