东亚钳蝎（Buthus martensii Karsch，BmK）又名马氏正钳蝎Mesobuthus martensii Karsch，属节肢动物门蛛形纲钳蝎科动物，遍布全世界各地^[1-2]。据统计，世界上大约有1 500种蝎子，其中12%是有毒的，对人类生命财产安全问题构成了极大地威胁^[3]。我国中医理论讲究“以毒攻毒”，通过传统方法炮制后的全蝎一般以整体入药，是我国名贵的中药材之一，有着息风镇痉、攻毒散结、通络镇痛等一系列功能，但中药方剂成分复杂多样，具体发挥治疗效果的成分难以明确^[4]。

转录组测序（RNA sequencing，RNA-seq）是使用高通量检测的方法对样本mRNA构建cDNA文库后进行测序，利用生物信息学对数据分析，对样本全部转录本情况进行研究，分析其基因表达水平，进而探究其具体分子生物学功能^[5-8]。和传统测序手段（如 Sanger 法测序、基因芯片）进行比较，RNA-seq的优势在于微量、分辨率高、价格低以及应用范围广等。，并且可以通过从头测序的方法对没有参考基因组的物种进行研究。随着转录组测序技术的向前发展和和进一步补充，转录组学也加快了中药现代化步伐。因此，本研究运用 Illumina 检测平台的转录组测序技术，将东亚钳蝎的前腹部、后腹部和毒液腺分别进行了转录组测序分析，筛选DEG（差异表达基因），进行差异基因的 GO 富集和 KEGG 通路分析^[9]。

1.   材料与方法

1.1   实验样品

实验动物捕捉自山东省临沂市兰陵县山区，经鉴定为东亚钳蝎，采集时间为2018年7月，选择体型健壮、完整且无病害的蝎子。

1.2   总RNA提取及文库构建

利用常规Trizol法从样品的组织中抽提总RNA，使用Nanodrop2000将所得到RNA的质量进行检测，利用琼脂糖凝胶电泳来检测RNA的完整程度，Agilent2100测定RIN值。检测符合质量标准的样品用于文库构建，将东亚钳蝎mRNA富集在含有Oligo dT的磁珠上，在Thermomixer中利用打断试剂把RNA切割成小片段，然后反转录，生成第一链cDNA，然后配置二链反应体系生成第二链，然后使用试剂盒把产物进行纯化回收，进而修复黏性末端，cDNA3’末端加碱基 A并连接接头并进行片段大小选择，最后进行PCR扩增反应。构建好的cDNA文库使用Agilent2100和ABI StepOnePlus Real Time PCR System进行质检，质检合格后使用Illumina/Hiseq-4000 RNA-seq进行上机测序^[10-11]。文库的构建与测序由深圳迪谱生物科技有限公司完成。

1.3   序列的拼接与组装

在研究没有参考基因组的转录组时，需要把所获取RNA-seq高质量测序数据进行从头组装，然后生成重叠群（contig）和单一序列（singleton），目的是为进一步进行基因产物表达、生物学功能研究等提供基础^[12]。本研究采用SeqPrep（https://github.com/jstjohn/Seqprep）对原始数据进行过滤，具体流程为：首先，将接头序列去掉，然后把3′末端质量较低的碱基修除，剔除含N百分比在10%以上的reads，舍弃adapter及质量修剪后长度小于30 bp的序列，得到clean date，使用Trinity（https://github.com/trinityrnaseq/trinityrnaseq）将所有的clean date进行denovo组装，拼接完成后使用TransRate（http://hibberdlab.com/transrate/）、CD-HIT（http://weizhongli-lab.org/cd-hit/）进行优化过滤。最后，使用BUSCO（Benchmarking universal single-copy orthologs，http://busco.ezlab.org）进行组装评估。

1.4   生物信息学分析

1.4.1   基因表达量的计算

Unigene表达量的计算使用FPKM法（Fragments Per kb per Million fragments），其计算公式为：FPKM = 10⁹C/NL

设FPKM（A）为Unigene A的表达量，则C为唯一比对到Unigene A的fragments数，N为唯一比对到所有Unigene的总fragments数，L为Unigene A的碱基数^[13-14]。

1.4.2   差异表达基因（DEGs）鉴定及功能富集分析

利用R语言（http://www.r-project.org/）对3个不同组织转录组样本进行差异基因鉴定。FDR ≤ 0.001 和 |log2Ratio| ≥ 1作为是否为差异基因的阈值。差异基因进行KO和GO功能富集分析。Blast2GO和WEGO进行GO富集，KAAS进行KO富集^[15-16]。

2.   结果

2.1   转录组测序、组装及注释

通过对前腹部组、后腹部组及毒液腺组的转录组测序，分别获得46977080，45682534，51173286个reads，Q30均在94%以上，符合进一步数据分析标准^[17]，见表1。

表  1  测序数据的统计结果

Table  1.  Summary of transcriptome sequencing date and transcriptome assembly

样本	原始数据	有效数据	碱基数目（G）	错误率（%）	Q20（%）	Q30（%）
前腹部	46977080	46593018	6945967490	0.02	98.26	94.29
后腹部	45682534	45222188	6735593318	0.02	98.27	94.3
毒液腺	51173286	50641874	7463605213	0.02	98.26	94.36

下载: 导出CSV 

| 显示表格

2.2   差异基因统计

如图所示，前腹部与毒液腺比对鉴定出差异基因数目为10411，与前腹部相比，毒液腺上调基因3912个（37.47%），下调基因6529个（62.53%）；后腹部与毒液腺相比差异基因数目为8258，与后腹部相比，毒液腺上调基因2338个（28.31%），下调基因5920个（71.69%）。两组共同差异基因为5294个^[18]，见图1。

图  1  不同东亚钳蝎基因表达的分析

A：前腹部 vs 毒液腺；B：后腹部 vs 毒液腺。红色、绿色为显著性差异基因；灰色为非显著性差异基因。横坐标表示基因的表达倍数变化，纵坐标表示统计学上基因表达差异的显著性水平；C. Venn图。

Figure  1.  Analysis of differentially expressed genes in Buthus martensii.

下载: 全尺寸图片 幻灯片

2.3   差异表达基因的GO功能分析

将所获得的两组共同差异基因的交集，共计5294个基因，放入到GO数据库做功能富集分析^[19]，结果显示，共有913个Unigene被注释在60个功能项中，其中314个Unigene归入23个生物学途径，最多的是生物学过程，共有71个，最少的是核酸代谢过程和细胞大分子生物合成过程，均只有1个；154个 Unigene 归入细胞组成，其中最多的是细胞组分，共有31个，最少的是intracellular organelle part和organelle part，均有4个；445个Unigene归入分子功能类，占比最大的是催化活性，共有102个^[11]，见表2。

表  2  差异表达基因的 GO 功能分类

Table  2.  GO functional categories of Buthus martensii Unigenes

GO序列号	数目	功能描述	GO分类
GO:0044260	8	cellular macromolecule metabolic process	BP
GO:0009987	26	cellular process	BP
GO:0090304	1	nucleic acid metabolic process	BP
GO:0044237	22	cellular metabolic process	BP
GO:0043170	15	macromolecule metabolic process	BP
GO:0016055	7	Wnt signaling pathway	BP
GO:0034641	14	cellular nitrogen compound metabolic process	BP
GO:0006177	4	GMP biosynthetic process	BP
GO:0046037	4	GMP metabolic process	BP
GO:0006807	19	nitrogen compound metabolic process	BP
GO:0008150	71	biological_process	BP
GO:0044763	19	single-organism cellular process	BP
GO:0007275	10	multicellular organism development	BP
GO:1901360	13	organic cyclic compound metabolic process	BP
GO:0006725	13	cellular aromatic compound metabolic process	BP
GO:0009059	1	macromolecule biosynthetic process	BP
GO:0046483	13	heterocycle metabolic process	BP
GO:0006139	13	nucleobase-containing compound metabolic process	BP
GO:1901070	4	guanosine-containing compound biosynthetic process	BP
GO:0009225	4	nucleotide-sugar metabolic process	BP
GO:0034645	1	cellular macromolecule biosynthetic process	BP
GO:1901135	17	carbohydrate derivative metabolic process	BP
GO:0005975	15	carbohydrate metabolic process	BP
GO:0005615	25	extracellular space	CC
GO:0044421	29	extracellular region part	CC
GO:0044424	17	intracellular part	CC
GO:0032991	5	macromolecular complex	CC
GO:0044464	31	cell part	CC
GO:0044446	4	intracellular organelle part	CC
GO:0044422	4	organelle part	CC
GO:0044444	6	cytoplasmic part	CC
GO:0043226	8	organelle	CC
GO:0043229	8	intracellular organelle	CC
GO:0005578	6	proteinaceous extracellular matrix	CC
GO:0043234	5	protein complex	CC
GO:0031012	6	extracellular matrix	CC
GO:0030414	23	peptidase inhibitor activity	MF
GO:0061134	23	peptidase regulator activity	MF
GO:0030234	23	enzyme regulator activity	MF
GO:0004866	22	endopeptidase inhibitor activity	MF
GO:0061135	22	endopeptidase regulator activity	MF
GO:0004857	23	enzyme inhibitor activity	MF
GO:0003676	11	nucleic acid binding	MF
GO:0008238	12	exopeptidase activity	MF
GO:0098772	23	molecular function regulator	MF
GO:0097159	35	organic cyclic compound binding	MF
GO:1901363	35	heterocyclic compound binding	MF
GO:0004177	9	aminopeptidase activity	MF
GO:0004553	10	hydrolase activity，hydrolyzing O-glycosyl compounds	MF
GO:0016798	10	hydrolase activity，acting on glycosyl bonds	MF
GO:0008237	14	metallopeptidase activity	MF
GO:0003938	3	IMP dehydrogenase activity	MF
GO:0042813	3	Wnt-activated receptor activity	MF
GO:0009032	2	thymidine phosphorylase activity	MF
GO:0016154	2	pyrimidine-nucleoside phosphorylase activity	MF
GO:0005102	8	receptor binding	MF
GO:0004697	3	protein kinase C activity	MF
GO:0046914	25	transition metal ion binding	MF
GO:0008241	2	peptidyl-dipeptidase activity	MF
GO:0003824	102	catalytic activity	MF

下载: 导出CSV 

| 显示表格

2.4   差异表达基因的KO富集分析

KEGG （Kyoto encyclopedia of genes and genomes）数据库可以将基因产物及其功能进行系统化分析注释^[20]。它对基因组、化学分子和生化代谢系统等方面的数据进行了整合了，包括代谢通路、天然药物、人类疾病、基因序列、生化反应及基因组等^[12]。利用该数据库有助于将基因及表达信息进行网络化研究。把获得的5294个Unigene放入KEGG数据库进行比对，显著富集的差异基因数量为1543个，为总体的29.15%。总共富集到108种代谢通路，其中Unigene最多的是钙离子信号通路，共注释到38个Unigene，占总体的2.46%；其次为补体系统途径，共注释到37个Unigene，占总体的2.40%；第三为胰液分泌途径，共注释到32个Unigene，占总体的2.07%。图2 为 Pathway富集性分析前 20 个通路的数据（图2）。为发现毒素相关基因，对转录组信息进行分析得出38个差异基因与钙离子通道相关，可能是潜在的毒素相关基因^[21]。有15个差异基因与磷脂酶代谢通路有关，可能是调节相关酶类代谢的基因。

图  2  气泡图展示了DEG富集程度最显著的前20个KEGG通路

Figure  2.  KEGG classification of Buthus martensii unigenes

下载: 全尺寸图片 幻灯片

3.   讨论

我国人民在不断地创新和实践中发展了中医药，辨证论治理论更是博大精深，但是由于传统中药成分极其复杂，具体作用机制不明确以及其他因素严重影响了中药向全世界进军的步伐^[22-23]。随着目前高通量RNA-Seq方法的不断改善和广泛应用，不仅影响了过去传统的研究方法，使转录组学，蛋白组学进入了全新的时代，而且在挖掘差异表达基因，阐释生物基因表达及调控机制中发挥了重要作用，从分子生物学角度为辨证论治理论提供强有力的基础，更加快了中药国际化的发展^[24]。

本研究通过构建东亚钳蝎前腹部、后腹部和毒液腺的转录组文库，获得大量的转录组信息。共获得高通量的转录组信息，各样品的 clean rate在 98%以上，Q 30 都在94%以上，表明测序数据质控水平较高，测序质量高，完整性好，符合接下来数据分析的标准。

为了更好地挖掘相关的目的基因，本研究将东亚钳蝎的毒液腺基因集分别于前腹部与后腹部基因集进行比较，然后取两两比较之后的交集进行生物信息学分析，这样可以更好的体现出毒液腺中基因功能。

GO（Gene ontology）数据库可以对基因和蛋白质功能进行标准化的注释，适用于不同物种，开发了三级结构的规范化语言来规范各种在线平台中基因及其产物的生物学功能注释，可将基因产物定义为相关分子功能（molecular function，MF）、生物学途径（biological process，BP）、细胞学组成（cellular component，CC）三大类^[9]。通过将DEG进行 GO 功能分类，归类到 60个大项中，含有 23个生物学途径，13个细胞组分和24个分子功能。深入对DEG进行 pathway 分析，进一步对DEG参与的主要代谢途径和信号转导途径进行研究，其中，富集最显著的前20个KEGG通路中与毒素密切相关的为钙离子信号通路、磷脂酶信号通路，更有利于明确差异表达基因在生物代谢周期中发挥的功能以及与其他基因的相互作用。