中国公共卫生  2019, Vol. 35 Issue (7): 876-880   PDF    
长角血蜱饥饿雌成蜱转录组测序及分析
崔莹莹1,2, 刘起勇2, 张玉栋1, 栗冬梅2, 宋秀平2, 王君2, 吴海霞2, 王俊刚1    
1. 石河子大学农学院,新疆石河子 832003;
2. 中国疾病预防控制中心传染病预防控制所
摘要目的 利用IlluminaHiSeq高通量技术对长角血蜱饥饿雌成蜱进行转录组测序。方法 将测序得到的数据进行拼接组装,并利用生物信息学方法对所得序列进行基因功能注释、功能分类、代谢途径分析和简单重复序列标记等分析。结果 共获得181 246 184个clean reads数据,组装后得到107 428个unigenes序列,平均长度1 246.29。以BLAST软件将所有的unigenes序列与Nr,Nt,Pfam,KOG,Swiss-prot,KEGG,GO数据库进行比对。与Nr数据库比对发现长角血蜱基因序列与同科肩突硬蜱(Ixodes scapularis)具有较高的同源性,为55.3 %;根据GO数据库的注释,将其分为生物过程、细胞组分和分子功能3大类共56分支;通过与KOG数据库的注释划分25类;根据与KEGG数据库的分析发现含有代谢通路相关的基因有32类,其中较多的是信号转导(signal transduction)类,占12.13 %。通过SSR位点查找共鉴定出45 863个SSR;SNP分析发现发生转换的SNP数为195 369个,发生颠换的SNP数为96 780个。结论 通过对长角血蜱饥饿雌成蜱转录组水平分析,可为后续长角血蜱功能基因的挖掘及表达谱等方面的研究奠定基础。
关键词长角血蜱     转录组     基因    
Sequencing and analysis on transcriptomes of unfed female Haemaphysalis longicornis
CUI Ying-ying, LIU Qi-yong, ZHANG Yu-dong, et al     
College of Agriculture, Shihezi University, Shihezi, Xinjiang Uygur Autonomous Region 832000, China
Abstract: Objective To sequence transcriptomes of unfed female Haemaphysalis longicornis using IlluminaHiSeq high-throughput technology. Methods The data on sequences the transcriptomes were spliced and assembled, and the obtained sequences were analyzed with functional annotation, functional classification, metabolic pathway analysis and simple repeated sequence markers using bioinformatics methods. Results A total of 181 246 184 clean reads data were obtained and 107 428 unigene sequences were obtained after assembly, with an average length of 1 246.29. All unigene sequences were aligned with the Nr, Nt, Pfam, KOG, Swiss-prot, KEGG, GO databases using BLAST software. Compared with the Nr database, the long-horned blood scorpion gene sequence has a high homology (55.3%) with that of Ixodes scapularis. According to the annotation of the GO database, the functions of the all unigene sequences were divided into 3 categories (biological process, cellular component and molecular function) covering 56 functional groups; based on to the annotations of the KOG database, all the unigene sequences were assigned into 25 categories; while, according to the analysis of the KEGG database, there are 32 groups of genes involved in metabolic pathways and a major part of them (12.13%) are involved in signal transduction. A total of 45 863 simple sequence repeats (SSRs) were identified with SSR locus search. Single nucleotide polymorphism (SNP) analysis indicated that the number of SNPs for base transition was 195 369 and that for base transversion was 96 780. Conclusion The analysis on transcriptomes of unfed female Haemaphysalis longicornis lays a foundation for subsequent researches on gene expression and expression of the tick.
Key words: Haemaphysalis longicornis     transcriptome     gene    

蜱是陆生脊椎动物的专性、非永久性体外寄生虫[1],是仅次于蚊虫的第二大人畜共患病的传播媒介和宿主。长角血蜱(Haemaphysalis longicornis)隶属于硬蜱科,血蜱属,是我国流行最广泛和最重要的蜱种之一,主要分布于东北、华北、西北等地的17个省市,寄生于马、狗、兔、牛、羊、鼠等多种动物体表,可携带斑点热立克次体、犬吉氏巴贝虫、牛卵圆巴贝虫、瑟氏泰勒虫等多种病原[25],引起家畜和人的疾病,对畜牧业发展和人类的健康造成重大经济损失。

近年来,随着分子生物和信息技术的逐渐发展,对长角血蜱的研究聚集在功能基因的发掘上,这大大补充和丰富了长角血蜱基因资源,对进一步了解其基因功能、结构及代谢具有重大意义。其中孙文英通过构建长角血蜱孤雌生殖种群cDNA表达文库,克隆了3条防御素基因[6];赵华渠利用蛋白质组学方法筛选出影响长角血蜱吸血的驱动蛋白和动力蛋白[7];刘建刚等也通过构建长角血蜱雄性全蜱抑制消减杂交文库,获得了5个有潜在生物学功能的表达序列标签(expressed sequence tag, EST)[8];以及Mulenga等人基于基因的保守区从长角血蜱中扩增出丝氨酸蛋白酶、半胱氨酸蛋白酶cDNA等[9]。但这些技术不能对物种基因进行全面、系统的分析,也不能对表达基因进行准确的定量研究。随着新一代高通量转录组测序Illumina技术的运用,可以在较短时间内获得某一物种大量的数据,能更全面地揭示生物个体在特定时期和特定组织的全局基因的表达情况,可用于研究功能基因、结构、新转录本预测等[10]。是一个通量高、成本低,可获得低丰度的表达基因,适用于未知基因组序列的物种发掘功能基因的重要途径[11]。本文通过对长角血蜱饥饿雌成蜱进行转录组测序,并对转录组水平进行系统分析,为后续对长角血蜱进行功能基因、蛋白质、代谢物等方面的深入研究提供依据,并在分子层面上为蜱病控制、生物学特性、传病机制等提供参考。

1 材料与方法 1.1 材料

选取室内温度(25 ± 1)℃、相对湿度(70 ± 10)%,实验室连续多代培养的世代整齐的长角血蜱饥饿雌成蜱为研究材料。

1.2 样品处理及转录组数据库分析流程

每组(共3组)取5只长角血蜱饥饿雌成蜱,利用QIAGEN试剂盒提取总RNA,并用带有Oligo(dT)的磁珠富集mRNA。随后加入fragmentation buffer将mRNA打断成短片段,以打断后的mRNA为模板合成第一链cDNA,然后加入缓冲液、dNTPs和DNA polymerase I和RNase H合成第二链cDNA,再用AMPure XP beads纯化双链cDNA,进行末端修复、加A尾并连接测序接头,再用AMPure XP beads进行片段大小选择。最后进行PCR扩增后构建cDNA文库。文库构建完成后,使用Qubit2.0进行初步定量,随后使用Agilent 2100和Q-PCR对文库进行质检。最后用IlluminaHiSeq进行测序,并将序列数据进行拼接组装及注释分析。

1.3 测序原始序列数据处理与过滤

测序得到的原始图像数据文件经CASAVA碱基识别(base calling)分析转化为原始测序序列(sequenced reads),即raw data或raw reads,其中包含测序序列(reads)的序列信息以及其对应的测序质量信息。为了保证信息分析质量,对raw reads过滤,得到clean reads,数据处理的步骤包括:①去除带接头(adapter)的reads;②去除N(N表示无法确定碱基信息)的比例大于10 % 的reads;③去除低质量reads(质量值Qphred ≤ 20的碱基数占整个reads的50 % 以上的reads)。

1.4 转录本拼接及聚类

本文通过对长角血蜱雌成蜱无参考基因组的测序,获得clean reads后,采用Trinity[12]对clean reads进行拼接。以Corset[13]层次聚类后得到最长转录本(unigene)进行后续的分析。

1.5 基因功能注释

为获得全面的基因功能信息,本文进行了七大数据库的基因功能注释,包括:Nr,Nt,Pfam,KOG,Swiss-prot,KEGG,GO[14]

2 结 果 2.1 转录组测序产量及组装质量分析(表1
表 1 转录组测序产量统计

长角血蜱饥饿雌成蜱3个样品(HaeA1、HaeA2、HaeA3)的转录组测序结果一共得到186 226 470个total raw reads,经过滤,剔除低质量的数据后,剩余181 246 184个clean reads。每个样品超过了8.86 G的clean bases。测序结果的GC含量正常在59.14 %~59.41 % 之间。不同样品的质量值Q30大于88.49 %,说明本次测序准确度较高(Q30是碱基识别出错的概率的整数映射,Q30越高表明碱基识别越可靠,准确度越高)。

组装后得到107 428个unigenes序列,长度 > 500 bp的数目为70 890个,平均长度为1 246.29,N50的长度为2, 087;对于transcripts共得到143 809个序列,长度 > 500 bp的数目为71 440个,平均长度为1 012.32,N50的长度为1 910,组装完整性较高(将序列按照长度递减累加,当累加之和刚好大于总长度的一半时,最后被累加的那条序列长度,即为N50,数值越大说明组装的质量越好)。

2.2 基因功能注释

通过BLAST软件将组装得到的所有unigenes序列与Nr,Nt,Pfam,KOG,Swiss-prot,KEGG,GO数据库进行比对,获得unigenes的注释信息。本次实验在以上7个数据库中至少1个数据库中注释成功的unigenes数目占总unigenes数量的52.66 % ;在以上7个数据库中都注释成功的unigenes数目占总unigenes数的比例为6.35 %。其中基因注释率高的数据库分别为占比41.61 % 的Nr数据库和占比41.57 % 的GO数据库,其次为Pfam数据库,基因注释占比为39.83 %。KOG、KEGG、Swiss-prot数据库注释基因占比分别为18.37 %,18.21 %。本次测序结果中还有较大一部分unigenes未获得注释,因此,对于长角血蜱及其相近物种在分子方面需进一步研究。

2.2.1 Nr功能注释

通过与Nr库比对注释,我们可以获取长角血蜱基因序列与近缘种基因序列的相似性。其中,72.7 % 的注释基因与已知蛋白有60 % 以上的相似性;通过E-value值分析,53 % 的注释基因具有很强的同源性(1e-60~1e-45),而17.3 % 的注释基因表现出很低的同源性(1e-15~1e-5);从物种比对分布来看,长角血蜱与肩突硬蜱(Ixodes scapularis)的比对重复率最高,达到55.3 %。因此,转录组序列中必然存在很大一部分蜱类特异的基因模型和基因序列,其相似度与蜱类以外的其他物种相比不是很高,但在蜱类物种内具有一定的相似程度。

2.2.2 GO注释与分类

根据Nr注释信息得到GO功能注释。该数据库包括分子功能(molecular function)、生物学过程(biological process)和细胞组分(cellular component)三大类,分别用来描述基因编码的产物所参与的生物过程、所具有的分子功能及所处的细胞环境。

在本次转录组unigene注释是对预测的长角血蜱unigene的功能进行分类。基于序列同源性,44 665个unigenes被分成三大主要类别共有56个功能群(图1)。统计每一类的基因数目发现,在生物学过程大类的26个分支中涉及本转录组unigene 124 758条,在细胞组分大类的20分支中涉及本转录组unigene 74 387个,在分子功能大类的10个分支中涉及本转录组unigene 55 381个。这56个功能群在二级分类中注释基因最多依次是结合功能(binding)25 892个;细胞过程(cellular process)25 546个;代谢过程(metabolic process)23 273个;单一生物过程(single-organism process)20 466个;催化活性(catalytic activity)20 106个。注释基因最少的是细胞外基质成分(extracellular matrix component)6个;拟核(nucleoid)5个;细胞聚集(cell aggregation)3个。

图 1 长角血蜱unigenes GO分类图

2.2.3 KOG注释和分类

KOG数据库通过对多种生物的蛋白质序列大量比较而来,是对基因进行同源分类的数据库,共有 25 个分类。在本次unigene注释结果中(图2),获得KOG注释的基因序列有 19 745 个,这些序列主要聚集在以下几个分类中:一般功能预测基因(general function prediction)3 238 个,信号转导机制(signal transduction mechanisms)2 731个,翻译后修饰及蛋白质转换(posttranslational modification,protein turnover,chaperones)1 933个,转录(transcription)1 349条,分别占的比例是16.40 %,13.83 %,9.79 %,6.83 %。核结构(nuclear structure)和细胞运动(cell motility)注释到的最少,仅有89个和36个,分别占0.45 %,0.18 %。但未知功能(function unknown)高达1 396个,为7.07 %,因此,开发这些未知基因的功能尤为重要。

图 2 长角血蜱Unigenes KOG分类图

2.2.4 KEGG分类

KEGG是系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库。KEGG已建立了一套完整KO注释的系统,可完成新测序物种的基因组或转录组的功能注释。根据参与的KEGG代谢通路可分为5个分支:细胞过程(cellular processes),环境信息处理(environmental information processing),遗传信息处理(genetic information processing),代谢(metabolism),有机系统(organismal systems)。本文将107 428个unigene映射到KEGG中,注释为代谢通路相关的序列共有19 570个(图3)。具有2373个基因的信号转导(signal transduction)是环境信息处理分支中丰度最高的亚类,占12.13 %,这可能与蜱参与气味物质的结合和转运有关。其次是内分泌系统(endocrine system)1 366个,占6.98 %。

图 3 长角血蜱unigenes KEGG分类图

2.3 长角血蜱雌成蜱转录组差异表达分析

通过对3个样品进行差异表达基因数目统计,其中HlonA1与HlonA2之间上调的基因为1 513个,下调的基因为1 270个;HlonA1与HlonA3之间上调的基因为1 482个,下调的基因为1 325个;而HlonA2与HlonA3之间差异表达的基因较少,其中上调的为59个,下调的为89个。因此,在长角血蜱不同个体之间也会出现不同程度的基因差异表达。其中一些基因可能与嗅觉相关,这将为后续蜱嗅觉基因相关研究提供线索。

对于长角血蜱饥饿雌成蜱转录组总长度135 820 526 bp的107 428个基因序列进行SSR分析,在27 419个基因中共鉴定出45 863个SSRs,有1个以上SSR序列动物基因有9 899个,其中化合物形成中存在的SSR数为4 014个。对于长角血蜱雌成蜱转录组SNP分析,发生转换的SNP数为195 369个,其中C/T数量为97 145个,出现的频率为0.72,A/G数量为98 224个,出现的频率也为0.72;而发生颠换的SNP数为96 780个,其中A/T数量为16 923个,出现的频率为0.12,A/C数量为25 301个,出现的频率为0.19,T/G数量为26 079个,出现的频率为0.19,C/G数量为28 477个,出现的频率为0.21。

3 讨 论

Illumina高通量测序的数据量大,效率高,试验成本低,适用于没有开展全基因组测序的物种开展转录组测序研究。本研究采用了Illumina高通量测序技术对长角血蜱饥饿雌成蜱进行了测序,序列拼接后得到的 107 428个unigenes序列,采用生物信息学分析方法,对基因进行注释,其中 56 576 个unigenes得到功能注释,所占比例为 52.66 %,这一结果显示采用高通量测序可以挖掘长角血蜱饥饿雌成蜱表达的大量基因,也为进一步研究长角血蜱基因克隆及验证提供参考。

简单重复序列标记(simple sequence repeat,SSR),又称为短串联重复序列或微卫星标记。广泛均匀分布于真核生物基因组中。是一类由几个核苷酸(1~6个)为重复单位组成的长达几十个核苷酸的重复序列,长度较短。本文于长角血蜱雌成蜱转录组中共发现了 45 863 个SSR位点,并利用Primer3设计了SSR引物。对SNP(single nucleotide polymorphisms)遗传标记,在转录组数据中发生转换和颠换的SNP数分别为195 369个和96 780个。可为后续长角血蜱多态性分析及开发分子标记提供丰富的数据。

本研究仅初步对长角血蜱饥饿雌成蜱进行了转录组分析,后续还将对长角血蜱雄成蜱,以及不同组织进行转录组测序分析,进一步完善长角血蜱不同时期,不同组织基因功能的注释及分析。

志谢 感谢中国疾病预防控制中心传染病预防控制所信息室张雯老师的帮助

参考文献
[1] 余增莹. 长角血蜱和嗜群血蜱ITS-2、COI和COⅡ基因的序列分析研究[D]. 四川农业大学, 2008.
[2] Jiang XL, Wang XJ, Li JD, et al. Isolation, identification and characterization of SFTS bunyavirus from ticks collected on the surface of domestic animals[J]. Chinese Journal of Virology, 2012, 28(3): 252–257.
[3] Park SW, Song BG, Shin EH, et al. Prevalence of severe fever with thrombocytopenia syndrome virus in Haemaphysalis longicornisticks in South Korea [J]. Ticks and Tick-Borne Diseases, 2014, 5(6): 975–977. DOI:10.1016/j.ttbdis.2014.07.020
[4] 张臣臣, 张仪. 我国常见蜱的生物学特性及其传播疾病研究进展[J]. 国际医学寄生虫病杂志, 2015, 42(3): 184–188. DOI:10.3760/cma.j.issn.1673-4122.2015.03.015
[5] Chen X, Xu S, Yu Z, et al. Multiple lines of evidence on the genetic relatedness of the parthenogenetic and bisexual Haema-physalis longicornis (Acari: Ixodidae) [J]. Infection, Genetics and Evolution, 2014, 21: 308–314. DOI:10.1016/j.meegid.2013.12.002
[6] 孙文英. 长角血蜱孤雌生殖种群防御素基因克隆、原核表达及功能分析[D]. 河北师范大学, 2016.
[7] 赵华渠. 长角血蜱唾液腺功能蛋白发掘与吸血适应机制研究[D]. 河北师范大学, 2017.
[8] 刘建刚, 刘光远, 田占成, 等. 长角血蜱雄蜱抑制消减杂交cDNA文库的构建与分析[J]. 中国兽医科学, 2009, 39(8): 659–664.
[9] da Silva Vaz I Jr, Torino Lermen T, Michelon A, et al. Effect of acaricides on the activity of a Boophilus microplus glutathione S-transferase [J]. Veterinary Parasitology, 2004, 119(2): 237–245.
[10] Josek T, Walden KKO, Allan BF, et al. A foreleg transcriptome for Ixodes scapularis ticks: candidates for chemoreceptors and binding proteins that might be expressed in the sensory Haller′s organ[J]. Ticks and Tick-borne Diseases, 2018.
[11] Maher CA, Kumar-Sinha C, Cao XH, et al. Transcriptome sequencing to detect gene fusions in cancer[J]. Nature, 2009, 458(7234): 97–101. DOI:10.1038/nature07638
[12] Grabherr MG, Haas BJ, Yassour M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nature Biotechnology, 2011, 29(7): 644. DOI:10.1038/nbt.1883
[13] Davidson NM, Oshlack A. Corset: enabling differential gene expression analysis for de novoassembled transcriptomes[J]. Genome Biology, 2014, 15(7): 410.
[14] Dillies MA, Rau A, Aubert J, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[J]. Briefings in Bioinformatics, 2013, 14(6): 671–683. DOI:10.1093/bib/bbs046