王騰宇 陳昕迪 翟 帥 石雅琴 劉春霞 王文龍*

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 獸醫(yī)學(xué)院/農(nóng)業(yè)部動(dòng)物疾病臨床診療技術(shù)重點(diǎn)實(shí)驗(yàn)室，呼和浩特 010018；2.內(nèi)蒙古農(nóng)業(yè)大學(xué) 生命科學(xué)院，呼和浩特 010018)

捻轉(zhuǎn)血矛線蟲(Haemonchuscontortus)是反芻動(dòng)物最重要的寄生蟲之一，給中國(guó)的畜牧業(yè)帶來(lái)了巨大的經(jīng)濟(jì)損失，但由于缺乏有效的疫苗，捻轉(zhuǎn)血矛線蟲主要依靠丙硫咪唑和伊維菌素等化學(xué)藥物進(jìn)行防治[1]。近20年的多項(xiàng)研究報(bào)告表明，捻轉(zhuǎn)血矛線蟲在多個(gè)國(guó)家呈現(xiàn)多重耐藥的流行趨勢(shì)[2]。有關(guān)捻轉(zhuǎn)血矛線蟲耐藥性研究國(guó)內(nèi)外有很多報(bào)道且組學(xué)研究不斷深入，但有關(guān)捻轉(zhuǎn)血矛線蟲耐藥性的非編碼RNA研究還處于初探階段[3]。

微小RNA(microRNA,miRNA)是一類重要的行使基因功能但又不編碼蛋白質(zhì)的一類調(diào)控分子，在植物和動(dòng)物的基因中普遍存在，在轉(zhuǎn)錄后調(diào)控中發(fā)揮重要作用，近年來(lái)逐漸成為研究的熱點(diǎn)[4]。miRNA是一類長(zhǎng)度約為21～24 nt的小分子單鏈非編碼RNA(Non-coding RNA，ncRNA),能夠和靶基因mRNA的3’UTR區(qū)互補(bǔ)或部分互補(bǔ)，引起靶mRNA的降解或者抑制其翻譯，從而對(duì)基因的表達(dá)進(jìn)行轉(zhuǎn)錄后水平的調(diào)控[5]。miRNA具有很強(qiáng)的保守性，在廣泛的生物學(xué)過(guò)程中發(fā)揮著重要作用，包括發(fā)育、細(xì)胞分化、增值和凋亡等[6]。1993年Lee等[7]在秀麗隱桿線蟲中發(fā)現(xiàn)miRNA。Reinhart等[8]發(fā)現(xiàn)miRNA對(duì)線蟲的生長(zhǎng)發(fā)育有重要的調(diào)節(jié)作用。趙錦燕等[9]通過(guò)miRNA表達(dá)譜芯片篩選肝癌耐藥細(xì)胞株中差異表達(dá)的miRNA，發(fā)現(xiàn)肝癌多藥耐藥的發(fā)生與 miRNA 表達(dá)異常有關(guān)。但在寄生蟲耐藥方面的miRNA研究較人類的耐藥miRNA研究還是有很大差別，對(duì)捻轉(zhuǎn)血矛線蟲耐藥miRNA的研究也只是處于初期探索階段。

本研究通過(guò)RNA-seq技術(shù)對(duì)miRNA和mRNA關(guān)聯(lián)分析，旨在探究捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株中的差異基因和miRNA，初步找到耐藥相關(guān)抗性通路和可能的耐藥調(diào)控機(jī)制，為捻轉(zhuǎn)血矛線蟲耐藥性機(jī)制的進(jìn)一步研究提供基礎(chǔ)。

1 材料與方法

1.1 捻轉(zhuǎn)血矛線蟲丙硫咪唑耐藥蟲株和敏感蟲株的采集

本研究用捻轉(zhuǎn)血矛線蟲丙硫咪唑耐藥蟲株采自內(nèi)蒙古自治區(qū)興安盟烏蘭浩特市科爾沁右翼前旗，通過(guò)進(jìn)行體內(nèi)驅(qū)蟲試驗(yàn)，蟲卵減少試驗(yàn)和對(duì)照試驗(yàn)確定耐藥蟲株[10]，對(duì)選取的綿羊進(jìn)行剖檢，采集蟲體并鑒定蟲種，并將捻轉(zhuǎn)血矛線蟲在體視鏡下逐條鑒別雌雄，分裝入凍存管放于液氮罐中冷凍帶回實(shí)驗(yàn)室。捻轉(zhuǎn)血矛線蟲敏感蟲株由內(nèi)蒙古農(nóng)業(yè)大學(xué)寄生蟲實(shí)驗(yàn)室通過(guò)體內(nèi)和體外試驗(yàn)自行分離獲得。

1.2 總RNA提取及質(zhì)檢

選取捻轉(zhuǎn)血矛線蟲敏感蟲株雌蟲和雄蟲各2管，每管50條；耐藥蟲株雌蟲和雄蟲各1管，每管60條。分別稱取相同質(zhì)量敏感蟲株和耐藥蟲株雌蟲和雄蟲進(jìn)行混合，按照Trizol試劑盒說(shuō)明書分別提取敏感蟲株和耐藥蟲株的RNA并進(jìn)行純化，使用NanoDrop 2 000和Agilent 2100 RNA 6 000 Nano kit對(duì)RNA的純度以及濃度進(jìn)行檢測(cè)，質(zhì)量濃度高于500 ng/μL的樣品RIN值≥7作為選取標(biāo)準(zhǔn)，同時(shí)根據(jù)Agilent 2 100檢測(cè)RNA的完整性。

1.3 Small RNA的文庫(kù)構(gòu)建

對(duì)提取的Total RNA進(jìn)行質(zhì)檢合格后，對(duì)PAGE膠電泳切膠選取18～30 nt范圍的條帶，對(duì)Small RNA進(jìn)行回收。分別對(duì)Small RNA片段的3’和5’端加上接頭進(jìn)行反轉(zhuǎn)錄并進(jìn)行PCR擴(kuò)增。將PCR產(chǎn)物進(jìn)行膠回收并純化140 bp左右的條帶，溶于EB solution，完成文庫(kù)構(gòu)建。將構(gòu)建好的文庫(kù)使用Agilent 2100以及ABI StepOnePlus Real-Time PCR System進(jìn)行質(zhì)量和產(chǎn)量的檢測(cè)。最后采用Illumina Hiseq2000平臺(tái)按照標(biāo)準(zhǔn)化流程完成文庫(kù)上機(jī)測(cè)序。

1.4 測(cè)序數(shù)據(jù)的處理與分析

1.4.1原始數(shù)據(jù)過(guò)濾與處理

(1)將得到的Small RNA測(cè)序數(shù)據(jù)進(jìn)行初步過(guò)濾，去除低質(zhì)量reads和不帶有3’接頭以及含有5’接頭的reads，過(guò)濾掉含有polyA的reads和不含插入片段的reads和插入片段長(zhǎng)度小于18 nt或者大于30 nt的reads，得到高質(zhì)量reads用于進(jìn)一步分析。

(2)選取Rfam數(shù)據(jù)庫(kù)對(duì)測(cè)序得到的Small RNA進(jìn)行注釋，利用Bowtie軟件將clean reads比對(duì)捻轉(zhuǎn)血矛線蟲參考基因組，確定測(cè)序得到的tag來(lái)自基因組的具體位置。使用軟件RepeatMasker version open-4.0.6找出repeat associate的Small RNA序列，將所有測(cè)到miRNA與miRNA權(quán)威數(shù)據(jù)庫(kù)miRBase進(jìn)行比對(duì)，鑒定到該物種的miRNA為已存在的miRNA,非本物種的miRNA為已知的miRNA，同時(shí)比對(duì)miRNA前體序列，過(guò)濾去除鑒定為來(lái)自mRNA降解片段，repeat區(qū)域，以及rRNA、scRNA、snRNA和tRNA等其他小RNA的tag序列，挑選出能比對(duì)上參考基因組的tag。

1.4.2數(shù)據(jù)分析

(1)由于miRNA有其特殊的二級(jí)結(jié)構(gòu)，使用軟件MIREAP(V0.2)預(yù)測(cè)miRNA的特殊的二級(jí)結(jié)構(gòu)，并找出可能存在的新miRNA。使用RNAhybrid(V6.01)，Miranda(V3.3a)，TargetScan(V7.0) 3種方法進(jìn)行靶基因預(yù)測(cè)，然后取3種方法得到的靶基因預(yù)測(cè)的結(jié)果的交集作為miRNA靶基因預(yù)測(cè)的結(jié)果。

(2)對(duì)敏感蟲株和耐藥蟲株的miRNA進(jìn)行表達(dá)分析，采用edgeR 軟件對(duì)miRNA進(jìn)行差異分析，表達(dá)量變化2倍以上并且P<0.05作為篩選的標(biāo)準(zhǔn)閾值。對(duì)差異miRNA對(duì)應(yīng)的靶基因向GO數(shù)據(jù)庫(kù)的GO terms進(jìn)行映射，統(tǒng)計(jì)注釋結(jié)果。同時(shí)進(jìn)行靶基因的KEGG Pathway顯著性富集分析，經(jīng)過(guò)FDR校正后，選擇q≤0.05的Pathway定義為在候選基因中顯著富集的Pathway。

2 結(jié)果與分析

2.1 Small RNA測(cè)序數(shù)據(jù)質(zhì)控及數(shù)據(jù)評(píng)估

根據(jù)得到的測(cè)序數(shù)據(jù)顯示，敏感蟲株和耐藥蟲株的樣本數(shù)據(jù)過(guò)濾掉低質(zhì)量序列和去接頭后，分別獲得14 916 974和13 888 536個(gè)clean reads，占reads總數(shù)的91.83%和94.62%，表明測(cè)序質(zhì)量較高，數(shù)據(jù)準(zhǔn)確可靠。將HC-2和HC-XA-aBZ序列注釋到Rfam數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，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)，盡可能發(fā)現(xiàn)并去除樣本中可能存在的rRNA、scRNA、snoRNA、snRNA和tRNA。

表1 敏感蟲株和耐藥蟲株比對(duì)上Rfam中非編碼RNA的tag豐度統(tǒng)計(jì)Table 1 Comparison on tag abundance statistics of non-coding RNA in Rfam in sensitive and resistant strains

2.2 miRNA的鑒定及分類注釋

對(duì)敏感蟲株和耐藥蟲株的miRNA進(jìn)行鑒定分析，其中捻轉(zhuǎn)血矛線蟲的miRNA分別鑒定到132和133個(gè)，其他物種的miRNA分別鑒定到466和446個(gè)，新的miRNA分別鑒定到902和876個(gè)。對(duì)Small RNA序列長(zhǎng)度分布進(jìn)行分析，根據(jù)統(tǒng)計(jì)結(jié)果顯示，大多數(shù)序列的長(zhǎng)度集中在16～28 nt，20～22 nt 長(zhǎng)度的序列頻率較高，其中在22 nt出現(xiàn)了波峰，即序列頻率最高(圖1)。對(duì)敏感蟲株和耐藥蟲株不同分類的序列豐度統(tǒng)計(jì)，在敏感蟲株測(cè)得的序列中，miRNA占78.20%；耐藥蟲株測(cè)得的序列中，miRNA占83.59%(圖2)。

圖1 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株Small RNA長(zhǎng)度分布統(tǒng)計(jì)Fig.1 Statistical analysis of Small RNA length distribution in H. contortus sensitive and resistant strains of albendazole

圖2 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株(a)和耐藥蟲株(b)不同分類tag豐度統(tǒng)計(jì)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差異表達(dá)分析

通過(guò)對(duì)HC-2和Hc-XA-aBZ鑒定得到的miRNA進(jìn)行匯總，運(yùn)用基于負(fù)二項(xiàng)分布的DEseq2，根據(jù)read counts計(jì)算每個(gè)miRNA的TPM表達(dá)量，獲得全部miRNA的表達(dá)譜。根據(jù)捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株間miRNA表達(dá)的差異倍數(shù)log2(Fold change)>1或log2(Fold change)<1以及顯著水平篩選出差異表達(dá)的miRNA，具體以FDR<0.05為篩選條件。從edgeR的分析結(jié)果表明，HC-2和Hc-XA-aBZ在篩選條件下共發(fā)現(xiàn)294個(gè)差異顯著的miRNA，其中上調(diào)的miRNA有113個(gè)，下調(diào)的miRNA有181個(gè)(圖3)。

2.4 差異miRNA靶基因的GO富集注釋

對(duì)捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株差異顯著的miRNA負(fù)相關(guān)靶基因進(jìn)行GO功能注釋。根據(jù)比對(duì)GO數(shù)據(jù)庫(kù)的結(jié)果顯示，顯著差異表達(dá)的miRNA的靶基因被注釋到了47個(gè)功能亞群，共1 430 條GO terms。其中細(xì)胞組分注釋到161條GO terms，包括細(xì)胞過(guò)程、代謝過(guò)程和單生物過(guò)程等；分子功能注釋到321條GO terms，包括細(xì)胞、細(xì)胞部分、膜、細(xì)胞器和大分子復(fù)合物等；生物學(xué)過(guò)程注釋到948條GO terms，包括催化活性、結(jié)合和運(yùn)輸活動(dòng)等(圖4)。

2.5 差異miRNA靶基因的KEGG富集注釋及分析

將miRNA預(yù)測(cè)到的差異靶基因進(jìn)行KEGG富集注釋，結(jié)果表明DEGs中富集到了胞吞作用(Endocytosis)，碳代謝(Carbon metabolism)，內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工(Protein processing in endoplasmic reticulum)，吞噬體(Phagosome)等327條通路上，這些通路中包括FoxO信號(hào)通路(FoxO signaling pathway)，mTOR信號(hào)通路(mTOR signaling pathway)，PI3K-Akt信號(hào)通路(PI3K-Akt signaling pathway)等與耐藥相關(guān)的一些抗性通路(圖5)。結(jié)合miRNA和靶基因的差異表達(dá)結(jié)果及兩者靶向關(guān)系分析結(jié)果，篩選出表達(dá)量負(fù)相關(guān)的miRNA-靶基因?qū)?表2)。

圖3 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株差異表達(dá)miRNA火山圖Fig.3 miRNA volcanic plot of differential expression between sensitive and resistant strains of H. contortus

圖4 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株差異miRNA靶基因的GO富集分析Fig.4 GO enrichment analysis of differential miRNA target genes in albendazole sensitive and resistant strains of H. contortus

圖5 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株差異miRNA靶基因KEGG富集分析Fig.5 Enrichment analysis of differential miRNA target gene KEGG between albendazole sensitive and resistant strains of H. contortus

表2 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株差異miRNA靶基因抗性相關(guān)的KEGG通路Table 2 KEGG pathway related to differential miRNA target gene between albendazolesensitive and resistant strains of H. contortus

表2(續(xù))

2.6 miRNA和mRNA的聯(lián)合分析

在動(dòng)物以及真核生物的研究當(dāng)中，miRNA與其所預(yù)測(cè)到的靶基因之間大部分都是以負(fù)相關(guān)的方式進(jìn)行調(diào)控。在本研究當(dāng)中，將預(yù)測(cè)到的miRNA和mRNA進(jìn)行關(guān)聯(lián)分析，共有1 770個(gè)miRNA-mRNA對(duì)(共涉及到274個(gè)差異的miRNA和603個(gè)差異的mRNA)的表達(dá)方式是負(fù)調(diào)控(圖6)。

圖6 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株miRNA-mRNA負(fù)相關(guān)網(wǎng)絡(luò)圖Fig.6 Negative correlation network of miRNA-mRNA between albendazole sensitive and resistant strains of H. contortus

2.7 RT-qPCR驗(yàn)證差異表達(dá)的miRNA

從測(cè)序結(jié)果中隨機(jī)挑選出6個(gè)差異表達(dá)的且參與靶向調(diào)控的miRNA，以U6作為內(nèi)參基因進(jìn)行RT-qPCR驗(yàn)證。驗(yàn)證結(jié)果顯示，hco-miR-5 896，hco-miR-5 948，hco-miR-86，hco-miR-55，novel-m1 032-3p，novel-m0 887-3p的相對(duì)表達(dá)水平與測(cè)序結(jié)果表達(dá)趨勢(shì)一致，具有較好的一致性(圖7)。

圖7 RT-qPCR驗(yàn)證RNA-Seq試驗(yàn)結(jié)果Fig.7 Results of RNA-Seq verification by RT-qPCR

3 討 論

根據(jù)本研究前期對(duì)捻轉(zhuǎn)血矛線蟲耐藥性流行病學(xué)調(diào)查結(jié)果來(lái)看，耐藥性已經(jīng)普遍存在[11]。早期研究表明，丙硫咪唑耐藥性的產(chǎn)生是由于捻轉(zhuǎn)血矛線蟲的I型β微管蛋白基因的單核苷酸突變引起的[12]。但近年來(lái)的研究發(fā)現(xiàn)，影響耐藥性的不僅僅是單一基因所導(dǎo)致[13]。研究表明，miRNA可以調(diào)控細(xì)胞周期，參與細(xì)胞凋亡和遷移，因此許多疾病都與異常miRNA的表達(dá)有關(guān)[14]。而miRNA主要通過(guò)與靶基因的3′UTR 結(jié)合抑制蛋白或mRNA的表達(dá)從而發(fā)揮調(diào)控的作用[15]。本研究主要以捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株為研究對(duì)象，將得到的miRNA進(jìn)行差異篩選并分析與其負(fù)調(diào)控的靶基因，發(fā)現(xiàn)多個(gè)在藥物代謝過(guò)程中起著重要調(diào)控作用的基因，如下調(diào)基因KLF4、ELO、PDHB和上調(diào)基因UGP,其靶向miRNA在人或小鼠等物種的耐藥性研究中均有報(bào)道[16-19]，本研究中的敏感蟲株和耐藥蟲株中也出現(xiàn)顯著差異表達(dá)，并通過(guò)GO分析顯示，在前10的富集條目中大多數(shù)與轉(zhuǎn)運(yùn)相關(guān)，在KEGG富集到的通路中，其中6條通路與代謝相關(guān)，這暗示著某些藥物代謝和轉(zhuǎn)運(yùn)通路對(duì)耐藥基因的產(chǎn)生具有很重要的作用。

Xu等[20]研究發(fā)現(xiàn)miR-22過(guò)表達(dá)可通過(guò)PTEN/PI3K/AKT信號(hào)通路增加順鉑在人胃腸道間質(zhì)瘤細(xì)胞中的化學(xué)敏感性，降低細(xì)胞對(duì)順鉑的耐藥性。本研究發(fā)現(xiàn)在敏感蟲株和耐藥蟲株比較組中miR-22存在較大差異，miR-22的靶向基因FZD4[21]已在其他物種的耐藥性研究中發(fā)現(xiàn)具有抑制細(xì)胞調(diào)亡等作用，與本研究結(jié)果一致。Shi等[22]報(bào)道m(xù)iR-29a在結(jié)腸癌細(xì)胞中也具有調(diào)控耐藥基因的功能，miR-29a通過(guò)調(diào)控PI3K/AKT通路的PTEN抑制細(xì)胞增殖和凋亡從而使耐藥性降低。本研究發(fā)現(xiàn)在耐藥蟲株中miR-29a的相對(duì)表達(dá)量較敏感蟲株差異顯著，有研究表明，miR-29a靶向的基因SEH1[23]和RAC1[24]具有抑制細(xì)胞分化和凋亡等作用，由此可以推測(cè)，F(xiàn)ZD4、SEH1和RAC1在捻轉(zhuǎn)血矛線蟲對(duì)丙硫咪唑耐藥作用中也具有抗性調(diào)節(jié)的作用。

Wang等[25]在對(duì)KLFs家族基因的研究中發(fā)現(xiàn)，KLF2的異位表達(dá)抑制了胃癌細(xì)胞的增殖、遷移和侵襲，在過(guò)表達(dá)后顯著促進(jìn)細(xì)胞凋亡。本研究發(fā)現(xiàn)耐藥蟲株中KLF2的表達(dá)量較敏感蟲株中KLF2的表達(dá)量明顯降低，KLF2基因靶向?qū)?yīng)到了hco-miR-83，hco-miR-87c，hco-miR-133，hco-miR-5910，hco-miR-5979等miRNA，并且hco-miR-83表達(dá)量明顯升高，基于miRNA-mRNA之間的負(fù)調(diào)控關(guān)系，由此可以推斷hco-miR-83在耐藥調(diào)控中發(fā)揮關(guān)鍵作用，為后續(xù)篩選耐藥相關(guān)基因和找到與其靶向調(diào)節(jié)負(fù)相關(guān)的miRNA提供了很有價(jià)值的信息。

4 結(jié) 論

本研究利用RNA-seq技術(shù)對(duì)捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株轉(zhuǎn)錄組和miRNA進(jìn)行測(cè)序分析，獲得miRNA表達(dá)譜，篩選出差異的miRNA和與之調(diào)控的差異基因，為從miRNA角度分析捻轉(zhuǎn)血矛線蟲耐藥機(jī)制提供了豐富的數(shù)據(jù)基礎(chǔ)和重要參考。