利用小RNA高通量测序筛查媒介昆虫携带的病原体*

张志毅 安小平 庄 璐 马麦卷 米志强 黄 勇 范 航 刘 玮 童贻刚

(军事医学科学院微生物流行病研究所，病原微生物生物安全国家重点实验室，北京 100071)

近年来，由于气候变化和环境破坏等因素，虫媒传染病疫情不时发生，其控制变得日益复杂(郑学礼，2011)。除了已熟知的虫媒传染病，20世纪70年代以来，多数年份都有一种或一种以上的新传染病被发现，新发传染病至今达40余种，而且其病原微生物种类复杂，有病毒、细菌、(包括立克次体、衣原体、螺旋体)及寄生虫等(宝福凯等，2009)。虫媒传染病病原体的本底调查，对于虫媒传染病防控有重要意义。

传统的病原体检测方法费时费力，尤其对于新发传染病，早期发现及诊断更为困难。测序方法是确定微生物种类的重要方法之一(秦楠等，2011)。传统的Sanger序列测定方法耗时长、成本高、通量低，无法满足对于未知病原的高通量分析。由于未知病原体的核酸序列未知，也不能够直接采用PCR 技术进行扩增和测序(李建彬等，2011)。近几年来，高通量测序技术，又称新一代测序技术(next-generation sequencing)的出现，使得测序费用相对于Sanger测序方法大幅下降，已广泛应用于动植物全基因组测序、基因组重测序、转录组测序、小RNAs测序和表观基因组测序等领域(岳桂东等，2012)，对于整个生命科学研究都产生巨大的影响，同时也为虫媒传染病病原体的快速检测提供了一条新的有效途径。

宏基因组(Metagenomics)高通量测序方法筛查病原体是一种新型的有效方法(Nakamuraetal.，2011)，但对于病毒病原而言，由于其基因组相对于原核生物和高等生物的基因组而言非常小，相对含量很低，其序列在测序数据分析时很容易被宿主序列所掩盖。小分子RNA(Small RNA，sRNA)是长度一般为20～30个核苷酸的非编码RNA 分子，无脊椎动物如蚕、蚊虫、线虫等应对RNA病毒入侵时可产生与其对应的小干扰RNA(siRNA)，昆虫siRNA免疫机制对入侵的RNA病毒的特定序列进行加工和放大，并利用这些小RNA分子抑制病毒基因的表达，从而达到控制病毒复制和致病的目的(Wuetal.，2010)。这种免疫机制使昆虫体内小RNA中病毒序列相对的比例更加突出。研究发现sRNA在真核细胞和原核细胞中对基因的转录后调控发挥着重要作用，因此动植物体内都有数量巨大的小分子RNA，如何从巨量的sRNA中确定病原体RNA或病原体相关sRNA，是一个值得探讨的科学问题，然而传统测序法操作复杂，花费大，测序深度有限，效率较低，难以解决这一问题，新一代高通量测序技术具有速度快、成本低、覆盖度深、产出巨大等优点，非常适合小分子RNA 测序(卫波等，2009)。已有研究报道小分子RNA高通量测序可以用来发现媒介昆虫携带病毒的报道(Janetal.，2009; Shietal.，2009；陈斌等，2011)。对于其他病原体，还未见相关报道。目前，具有传播媒介能力的节肢动物中以蚊类最多。已登记的535种虫媒病毒中，从蚊类分离到的病毒占近50%，其次是从蜱分离到的病原体(116种)，占虫媒病毒总数的21.68%(李文刚等，2011)。因此本文选取蚊和蜱为例，尝试了基于小分子RNA高通量测序筛查媒介昆虫携带病原体的方法。结果发现，小RNA病原体筛查范围其实不仅限于RNA病毒和DNA病毒，还可以发现原核以及真核生物病原体。现报道如下。

1 材料与方法

1.1 实验材料

蚊虫从云南西双版纳地区野外捕获，分别为3个种类：中华按蚊Anophelessinensis、三带喙库蚊Culextritaeniorhynchus和致倦库蚊Cx.quinquefasciatus，各100只左右。蜱采自黑龙江牡丹江市林区，包括全沟硬蜱Ixodespersulcatus18只、边缘革蜱Dermacentormarginatus19只和长角血蜱Haemaphysalislongicornis21只。

1.2 总RNA 的提取

液氮冷冻，同种的所有虫体取出研磨后各以500 μL 的Trizol (购自Invitrogen) 裂解，加入蛋白酶以65℃处理30 min后提取总RNA。具体步骤：在室温15～30℃放置5 min，加入氯仿，用力震荡约15 s，室温下放置2～3 min， 12 000 r/min(4℃) 离心15 min；取上层水相加入异丙醇，在室温下放置10 min， 12 000 r/min(2℃～8℃) 离心10 min；弃上清，加75%乙醇洗涤，涡旋混合，7 500 r/min(2～8℃) 离心5 min，弃上清; 让沉淀的RNA在室温下自然干燥5～10 min，用DEPC 水溶解。RNA 沉淀，55～60℃孵育10 min。

1.3 样品的制备与Solexa 测序

siRNA 分离与高通量测序由深圳华大基因公司完成，具体流程如下：首先用聚丙烯酰胺电泳分离siRNA，切胶回收18～30 nt 大小的片段，两端分别连接RNA 接头，然后进行反转录和PCR 扩增，最后将扩增产物上样进行Solexa 法测序。

1.4 数据过滤

初始数据采用软件Solexa pipeline进行以下几步过滤：1)去掉低质量的reads；2)去掉5′接头和3′接头序列；3)过滤无插入接头序列；4)过滤ployA reads。

1.5 判别方法

1.5.1序列比对。使用NCBI的Blast(2.22)软件中的blastn作核酸比对，下载NCBI核酸数据库(nt库2013年4月)建立本地化核酸数据库，以过滤后的测序reads作为查询项，进行全面比对搜索。只给出得分在前10位的结果。

1.5.2blast结果过滤。通过以下几步对结果进行过滤：每个结果只保留得分最高的，如果有并列最高的两项属于不同物种则认为属于不同物种间保守序列而予以忽略；去掉核糖体RNA(因为核糖体RNA有较高保守性且数量较多，因此特异性差)；去掉匹配长度小于21 bp的结果(因为reads序列较短，匹配短的很可能是随机匹配)。

1.5.3物种分析。因为NCBI只给出了gi号与种属分类编号的对应文件，从nt库的标题可以得到gi号与具体数据库编号(如ref、emb、dbj等)的对应关系。而blast结果只给出具体数据库编号，所以由具体数据库编号得到对应gi号，再得到种属分类编号，进行统计得到与各属匹配的reads数和匹配总长度，最后转化为属的名称。全部过程用自行编写的python程序完成。

1.5.4校正排序。因为小RNA序列较短，比对搜索的结果中必然有较多的非特异的随机匹配。考虑到nt库中搜集的某一种属的序列越多，被随机匹配到的可能性越大，因此我们又将最后得到每个种属的匹配总长度除以该种属在nt库中的总长度，得到检出比值，按照这个比值从大到小排序。这个比值越大的，表明该种属有更大的检出可能性。

2 结果

2.1 高通量测序统计结果

6只样品分别经过RNA提取，小RNA分离、纯化和测序，并经过数据过滤，得到有效小RNA序列的数量和测序的平均覆盖倍数估计值在表1中列出。

表1 各虫种测序reads数和平均覆盖倍数Tab.1 Number of reads and coverage of the 6 vector insects

2.2 Blast及结果过滤

将各样品的有效reads作为查询项，输入到nt库做Blast比对搜索。只保留得分在前10的结果。然后将每一条比对结果中得分最高的一项汇集起来进行统计。因为每株昆虫样品中，核糖体RNA的数量都占有较高比例，所以首先将结果中的核糖体RNA的条目去掉。这一步可去除11%～18%的匹配项。

小RNA长度较短，从表1可以看出，22 nt长的reads占了大多数。而nt库含有1 600多万条序列，总碱基长度达到37 594 899 229 nt，有许多研究详尽的物种如人、鼠、果蝇等物种的序列，基因组很大且收录得有冗余，因此搜索时，当长度小于等于20时，得到随机匹配的几率非常大。我们曾做过测试，设计了10 000条22 bp的随机序列，去和nt库比对，发现匹配在18～20 bp的条目占了接近60%的比例。因此，选取Blast结果中匹配碱基数大于20的条目能够进一步去除非特异随机匹配，使其不至于干扰真实结果。结果表明，可进一步去掉约2/3的条目。

2.3 病原体筛查结果

对剩下的Blast结果进行统计，首先根据每一条的数据库编号找到对应的GenBank序号(gi号)，并记下其匹配碱基数，这一步完成后，将gi号相同的条目的匹配碱基数相加，得到每个gi号对应匹配总长度之和。而后将gi号变为对应的物种分类编号，就得到所包含的每个物种的匹配总长度。结果中可以发现，由于人、小鼠等物种在nt库中总长度很大，所以其随机匹配几率很高，在结果中排名很靠前。为了排除这类非随机匹配的干扰，用每个种属的匹配长度除以该种属在nt库中的总长度的比值大小来排序，就可以使检出比超出背景值的物种排在前面。

在结果中，排在最前面的是所测昆虫的物种及其近缘物种，其数量占有绝对优势。此外一些环境微生物、植物、真菌也占有相当比例。除去此类一般意义的物种外，一些重要的病原体是我们最关注的。

根据以上算法，我们得到的初步结果中(表2)，中华按蚊可能携带有乙型脑炎病毒。乙型脑炎又称“日本脑炎”，主要通过蚊虫叮咬传播，库蚊、伊蚊、按蚊中均有一些种可以传播此病，主要流行区域为东南亚、西太平洋地区，我国除新疆、西藏、青海外均为乙脑疫区(郭杨等，2008)。三带喙库蚊则可能携带的病原体是库蚊浓核病毒(Brevidensovirus)。蚊浓核病毒属于细小病毒科，特异性感染蚊类，引发特异性的细胞核致密增生病变并可最终导致宿主发病或死亡(顾金保等，2008)。根据顾金保的统计，2007年以前发现的蚊浓核病毒，有两株是在中国发现的，分别是淡色库蚊浓核病毒(Culexpipienspallensdensovirus, CpDNV)和C6/36浓核病毒(C6/36 DNV)。我们样品中的病毒应该更接近于前者。库蚊黄病毒为只感染蚊虫的黄病毒(Kimetal.，2009)，有研究报道该病毒可与黄病毒属的其他病毒共同感染三带喙库蚊(Kentetal.，2010)，但是否相互影响还有待研究。致倦库蚊中发现的可能病原体则是果蝇X病毒，该病毒由两条核酸链组成，A链序列于1996年测定(Chungetal.，1996)，B链于2002年测定(Shwedetal.，2002)，但该病毒在研究昆虫的病毒免疫机制中应用较多，对人的致病性还待深入研究。全沟硬蜱最可能携带的病原体有：peacockii立克次体(Rickettsiapeacockii，Rp)。边缘革蜱最可能携带的病原体有：立氏立克次体(Rickettsiarickettsii，Rr)。长角血蜱最可能携带的病原体是：伯氏柯克斯体Coxiellaburnetii。R.peacockii和R.rickettsii都属于斑点热群立克次体，目前NCBI物种分类数据库收录的斑点热群立克次体不少于几十种。2005的一项调查研究，总结了我国之前引起立克次体病的10种病原体(张丽娟等，2005)，其中只有伯氏柯克斯体出现血蜱中，R.peacockii和R.rickettsii在国内较少有报道。伯氏柯克斯体习惯上称为Q热立克次体，它所引起的Q热是一种人兽共患病。

表2 蜱与蚊虫中可能含有的病原体Tab.2 Possible pathogens contained in ticks and mosquitos

3 讨论

常规病原体高通量检测需要对不同的病原体采用不同的纯化方案和测序方案，如细菌和原虫作为细胞生物，需要按照细胞生物的分离方法将它们与人体血液和组织成分分离，并提取DNA进行测序；而DNA病毒及支原体、衣原体需要按照可滤过性病原体进行分离和DNA样品提取；RNA病毒则需要采用滤过性提取方法提取和进行RNA测序。对于未知病原体，上述3种方案需要分别考虑，选择使用，最可靠的办法是同时采用上述3种方案进行高通量测序分析，这些操作在时间、人力和经济上都会造成很大负担。如果采用RNA高通量测序技术，则可以用一种测序方案检测所有类型的病原体，这不仅能够节省很多的人工和试剂费用，更重要的是可以节省大量时间，尤其是在有暴发疫情的时候，病原体筛查的速度将显得十分关键。小RNA高通量测序筛选媒介昆虫携带病毒已经被证明是一种有效的策略(Jan，2009; Wuetal.，2010)。本研究结果表明，用这一方法也不仅可以检测媒介昆虫携带的病毒，而且对细菌及原虫等病原体的检测也具有很大的潜力。因此用小RNA高通量测序筛查病原体，只需一次测序，一次分析就得到样品可能含有的病毒、原核及真核的各类病原体，具有快速、简便、节省费用的突出优点。

小RNA由于其在调控、免疫等方面的重要作用，当前受到格外关注，各类研究已积累了大量的小RNA高通量测序数据。因此用本文所述的方法对这些小RNA数据进行分析，可获得样品中可能含有的各类物种的信息，这对于充分利用实验数据获取有用信息具有一定意义，也可能对疾病的发生机制、发展进程与预后提供一些新的线索。

通过上述的生物信息学分析，大大缩小了后期实验验证的范围。但由于小分子RNA核苷酸序列长度较短，因此在对大的核酸数据库比对搜索时，非特异性的随机匹配可能性较大，因此有必要进一步优化生物信息学分析参数来尽可能排除掉比对结果中的各种非特异匹配，并需进一步做实验验证。