捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株miRNA-mRNA关联比较分析

王腾宇 陈昕迪 翟 帅 石雅琴 刘春霞 王文龙*

(1.内蒙古农业大学 兽医学院/农业部动物疾病临床诊疗技术重点实验室，呼和浩特 010018；2.内蒙古农业大学 生命科学院，呼和浩特 010018)

捻转血矛线虫(Haemonchuscontortus)是反刍动物最重要的寄生虫之一，给中国的畜牧业带来了巨大的经济损失，但由于缺乏有效的疫苗，捻转血矛线虫主要依靠丙硫咪唑和伊维菌素等化学药物进行防治[1]。近20年的多项研究报告表明，捻转血矛线虫在多个国家呈现多重耐药的流行趋势[2]。有关捻转血矛线虫耐药性研究国内外有很多报道且组学研究不断深入，但有关捻转血矛线虫耐药性的非编码RNA研究还处于初探阶段[3]。

微小RNA(microRNA,miRNA)是一类重要的行使基因功能但又不编码蛋白质的一类调控分子，在植物和动物的基因中普遍存在，在转录后调控中发挥重要作用，近年来逐渐成为研究的热点[4]。miRNA是一类长度约为21～24 nt的小分子单链非编码RNA(Non-coding RNA，ncRNA),能够和靶基因mRNA的3’UTR区互补或部分互补，引起靶mRNA的降解或者抑制其翻译，从而对基因的表达进行转录后水平的调控[5]。miRNA具有很强的保守性，在广泛的生物学过程中发挥着重要作用，包括发育、细胞分化、增值和凋亡等[6]。1993年Lee等[7]在秀丽隐杆线虫中发现miRNA。Reinhart等[8]发现miRNA对线虫的生长发育有重要的调节作用。赵锦燕等[9]通过miRNA表达谱芯片筛选肝癌耐药细胞株中差异表达的miRNA，发现肝癌多药耐药的发生与 miRNA 表达异常有关。但在寄生虫耐药方面的miRNA研究较人类的耐药miRNA研究还是有很大差别，对捻转血矛线虫耐药miRNA的研究也只是处于初期探索阶段。

本研究通过RNA-seq技术对miRNA和mRNA关联分析，旨在探究捻转血矛线虫敏感虫株和耐药虫株中的差异基因和miRNA，初步找到耐药相关抗性通路和可能的耐药调控机制，为捻转血矛线虫耐药性机制的进一步研究提供基础。

1 材料与方法

1.1 捻转血矛线虫丙硫咪唑耐药虫株和敏感虫株的采集

本研究用捻转血矛线虫丙硫咪唑耐药虫株采自内蒙古自治区兴安盟乌兰浩特市科尔沁右翼前旗，通过进行体内驱虫试验，虫卵减少试验和对照试验确定耐药虫株[10]，对选取的绵羊进行剖检，采集虫体并鉴定虫种，并将捻转血矛线虫在体视镜下逐条鉴别雌雄，分装入冻存管放于液氮罐中冷冻带回实验室。捻转血矛线虫敏感虫株由内蒙古农业大学寄生虫实验室通过体内和体外试验自行分离获得。

1.2 总RNA提取及质检

选取捻转血矛线虫敏感虫株雌虫和雄虫各2管，每管50条；耐药虫株雌虫和雄虫各1管，每管60条。分别称取相同质量敏感虫株和耐药虫株雌虫和雄虫进行混合，按照Trizol试剂盒说明书分别提取敏感虫株和耐药虫株的RNA并进行纯化，使用NanoDrop 2 000和Agilent 2100 RNA 6 000 Nano kit对RNA的纯度以及浓度进行检测，质量浓度高于500 ng/μL的样品RIN值≥7作为选取标准，同时根据Agilent 2 100检测RNA的完整性。

1.3 Small RNA的文库构建

对提取的Total RNA进行质检合格后，对PAGE胶电泳切胶选取18～30 nt范围的条带，对Small RNA进行回收。分别对Small RNA片段的3’和5’端加上接头进行反转录并进行PCR扩增。将PCR产物进行胶回收并纯化140 bp左右的条带，溶于EB solution，完成文库构建。将构建好的文库使用Agilent 2100以及ABI StepOnePlus Real-Time PCR System进行质量和产量的检测。最后采用Illumina Hiseq2000平台按照标准化流程完成文库上机测序。

1.4 测序数据的处理与分析

1.4.1原始数据过滤与处理

(1)将得到的Small RNA测序数据进行初步过滤，去除低质量reads和不带有3’接头以及含有5’接头的reads，过滤掉含有polyA的reads和不含插入片段的reads和插入片段长度小于18 nt或者大于30 nt的reads，得到高质量reads用于进一步分析。

(2)选取Rfam数据库对测序得到的Small RNA进行注释，利用Bowtie软件将clean reads比对捻转血矛线虫参考基因组，确定测序得到的tag来自基因组的具体位置。使用软件RepeatMasker version open-4.0.6找出repeat associate的Small RNA序列，将所有测到miRNA与miRNA权威数据库miRBase进行比对，鉴定到该物种的miRNA为已存在的miRNA,非本物种的miRNA为已知的miRNA，同时比对miRNA前体序列，过滤去除鉴定为来自mRNA降解片段，repeat区域，以及rRNA、scRNA、snRNA和tRNA等其他小RNA的tag序列，挑选出能比对上参考基因组的tag。

1.4.2数据分析

(1)由于miRNA有其特殊的二级结构，使用软件MIREAP(V0.2)预测miRNA的特殊的二级结构，并找出可能存在的新miRNA。使用RNAhybrid(V6.01)，Miranda(V3.3a)，TargetScan(V7.0) 3种方法进行靶基因预测，然后取3种方法得到的靶基因预测的结果的交集作为miRNA靶基因预测的结果。

(2)对敏感虫株和耐药虫株的miRNA进行表达分析，采用edgeR 软件对miRNA进行差异分析，表达量变化2倍以上并且P<0.05作为筛选的标准阈值。对差异miRNA对应的靶基因向GO数据库的GO terms进行映射，统计注释结果。同时进行靶基因的KEGG Pathway显著性富集分析，经过FDR校正后，选择q≤0.05的Pathway定义为在候选基因中显著富集的Pathway。

2 结果与分析

2.1 Small RNA测序数据质控及数据评估

根据得到的测序数据显示，敏感虫株和耐药虫株的样本数据过滤掉低质量序列和去接头后，分别获得14 916 974和13 888 536个clean reads，占reads总数的91.83%和94.62%，表明测序质量较高，数据准确可靠。将HC-2和HC-XA-aBZ序列注释到Rfam数据库进行比对，HC-2中rRNA占1.24%，snRNA占0.01%，snoRNA占0.00%，tRNA占0.49%；HC-XA-aBZ中rRNA占1.12%，snRNA占0.01%，snoRNA占0.00%，tRNA占0.21%(表1)，尽可能发现并去除样本中可能存在的rRNA、scRNA、snoRNA、snRNA和tRNA。

表1 敏感虫株和耐药虫株比对上Rfam中非编码RNA的tag丰度统计Table 1 Comparison on tag abundance statistics of non-coding RNA in Rfam in sensitive and resistant strains

2.2 miRNA的鉴定及分类注释

对敏感虫株和耐药虫株的miRNA进行鉴定分析，其中捻转血矛线虫的miRNA分别鉴定到132和133个，其他物种的miRNA分别鉴定到466和446个，新的miRNA分别鉴定到902和876个。对Small RNA序列长度分布进行分析，根据统计结果显示，大多数序列的长度集中在16～28 nt，20～22 nt 长度的序列频率较高，其中在22 nt出现了波峰，即序列频率最高(图1)。对敏感虫株和耐药虫株不同分类的序列丰度统计，在敏感虫株测得的序列中，miRNA占78.20%；耐药虫株测得的序列中，miRNA占83.59%(图2)。

图1 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株Small RNA长度分布统计Fig.1 Statistical analysis of Small RNA length distribution in H. contortus sensitive and resistant strains of albendazole

图2 捻转血矛线虫丙硫咪唑敏感虫株(a)和耐药虫株(b)不同分类tag丰度统计Fig.2 Statistical analysis of tag abundance in different classifications of albendazole sensitive strains (a) and resistant strains (b) of H. contortus

2.3 miRNA差异表达分析

通过对HC-2和Hc-XA-aBZ鉴定得到的miRNA进行汇总，运用基于负二项分布的DEseq2，根据read counts计算每个miRNA的TPM表达量，获得全部miRNA的表达谱。根据捻转血矛线虫敏感虫株和耐药虫株间miRNA表达的差异倍数log2(Fold change)>1或log2(Fold change)<1以及显著水平筛选出差异表达的miRNA，具体以FDR<0.05为筛选条件。从edgeR的分析结果表明，HC-2和Hc-XA-aBZ在筛选条件下共发现294个差异显著的miRNA，其中上调的miRNA有113个，下调的miRNA有181个(图3)。

2.4 差异miRNA靶基因的GO富集注释

对捻转血矛线虫敏感虫株和耐药虫株差异显著的miRNA负相关靶基因进行GO功能注释。根据比对GO数据库的结果显示，显著差异表达的miRNA的靶基因被注释到了47个功能亚群，共1 430 条GO terms。其中细胞组分注释到161条GO terms，包括细胞过程、代谢过程和单生物过程等；分子功能注释到321条GO terms，包括细胞、细胞部分、膜、细胞器和大分子复合物等；生物学过程注释到948条GO terms，包括催化活性、结合和运输活动等(图4)。

2.5 差异miRNA靶基因的KEGG富集注释及分析

将miRNA预测到的差异靶基因进行KEGG富集注释，结果表明DEGs中富集到了胞吞作用(Endocytosis)，碳代谢(Carbon metabolism)，内质网中的蛋白质加工(Protein processing in endoplasmic reticulum)，吞噬体(Phagosome)等327条通路上，这些通路中包括FoxO信号通路(FoxO signaling pathway)，mTOR信号通路(mTOR signaling pathway)，PI3K-Akt信号通路(PI3K-Akt signaling pathway)等与耐药相关的一些抗性通路(图5)。结合miRNA和靶基因的差异表达结果及两者靶向关系分析结果，筛选出表达量负相关的miRNA-靶基因对(表2)。

图3 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株差异表达miRNA火山图Fig.3 miRNA volcanic plot of differential expression between sensitive and resistant strains of H. contortus

图4 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株差异miRNA靶基因的GO富集分析Fig.4 GO enrichment analysis of differential miRNA target genes in albendazole sensitive and resistant strains of H. contortus

图5 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株差异miRNA靶基因KEGG富集分析Fig.5 Enrichment analysis of differential miRNA target gene KEGG between albendazole sensitive and resistant strains of H. contortus

表2 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株差异miRNA靶基因抗性相关的KEGG通路Table 2 KEGG pathway related to differential miRNA target gene between albendazolesensitive and resistant strains of H. contortus

表2(续)

2.6 miRNA和mRNA的联合分析

在动物以及真核生物的研究当中，miRNA与其所预测到的靶基因之间大部分都是以负相关的方式进行调控。在本研究当中，将预测到的miRNA和mRNA进行关联分析，共有1 770个miRNA-mRNA对(共涉及到274个差异的miRNA和603个差异的mRNA)的表达方式是负调控(图6)。

图6 捻转血矛线虫丙硫咪唑敏感虫株和耐药虫株miRNA-mRNA负相关网络图Fig.6 Negative correlation network of miRNA-mRNA between albendazole sensitive and resistant strains of H. contortus

2.7 RT-qPCR验证差异表达的miRNA

从测序结果中随机挑选出6个差异表达的且参与靶向调控的miRNA，以U6作为内参基因进行RT-qPCR验证。验证结果显示，hco-miR-5 896，hco-miR-5 948，hco-miR-86，hco-miR-55，novel-m1 032-3p，novel-m0 887-3p的相对表达水平与测序结果表达趋势一致，具有较好的一致性(图7)。

图7 RT-qPCR验证RNA-Seq试验结果Fig.7 Results of RNA-Seq verification by RT-qPCR

3 讨 论

根据本研究前期对捻转血矛线虫耐药性流行病学调查结果来看，耐药性已经普遍存在[11]。早期研究表明，丙硫咪唑耐药性的产生是由于捻转血矛线虫的I型β微管蛋白基因的单核苷酸突变引起的[12]。但近年来的研究发现，影响耐药性的不仅仅是单一基因所导致[13]。研究表明，miRNA可以调控细胞周期，参与细胞凋亡和迁移，因此许多疾病都与异常miRNA的表达有关[14]。而miRNA主要通过与靶基因的3′UTR 结合抑制蛋白或mRNA的表达从而发挥调控的作用[15]。本研究主要以捻转血矛线虫敏感虫株和耐药虫株为研究对象，将得到的miRNA进行差异筛选并分析与其负调控的靶基因，发现多个在药物代谢过程中起着重要调控作用的基因，如下调基因KLF4、ELO、PDHB和上调基因UGP,其靶向miRNA在人或小鼠等物种的耐药性研究中均有报道[16-19]，本研究中的敏感虫株和耐药虫株中也出现显著差异表达，并通过GO分析显示，在前10的富集条目中大多数与转运相关，在KEGG富集到的通路中，其中6条通路与代谢相关，这暗示着某些药物代谢和转运通路对耐药基因的产生具有很重要的作用。

Xu等[20]研究发现miR-22过表达可通过PTEN/PI3K/AKT信号通路增加顺铂在人胃肠道间质瘤细胞中的化学敏感性，降低细胞对顺铂的耐药性。本研究发现在敏感虫株和耐药虫株比较组中miR-22存在较大差异，miR-22的靶向基因FZD4[21]已在其他物种的耐药性研究中发现具有抑制细胞调亡等作用，与本研究结果一致。Shi等[22]报道miR-29a在结肠癌细胞中也具有调控耐药基因的功能，miR-29a通过调控PI3K/AKT通路的PTEN抑制细胞增殖和凋亡从而使耐药性降低。本研究发现在耐药虫株中miR-29a的相对表达量较敏感虫株差异显著，有研究表明，miR-29a靶向的基因SEH1[23]和RAC1[24]具有抑制细胞分化和凋亡等作用，由此可以推测，FZD4、SEH1和RAC1在捻转血矛线虫对丙硫咪唑耐药作用中也具有抗性调节的作用。

Wang等[25]在对KLFs家族基因的研究中发现，KLF2的异位表达抑制了胃癌细胞的增殖、迁移和侵袭，在过表达后显著促进细胞凋亡。本研究发现耐药虫株中KLF2的表达量较敏感虫株中KLF2的表达量明显降低，KLF2基因靶向对应到了hco-miR-83，hco-miR-87c，hco-miR-133，hco-miR-5910，hco-miR-5979等miRNA，并且hco-miR-83表达量明显升高，基于miRNA-mRNA之间的负调控关系，由此可以推断hco-miR-83在耐药调控中发挥关键作用，为后续筛选耐药相关基因和找到与其靶向调节负相关的miRNA提供了很有价值的信息。

4 结 论

本研究利用RNA-seq技术对捻转血矛线虫敏感虫株和耐药虫株转录组和miRNA进行测序分析，获得miRNA表达谱，筛选出差异的miRNA和与之调控的差异基因，为从miRNA角度分析捻转血矛线虫耐药机制提供了丰富的数据基础和重要参考。