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

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

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

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

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

1 材料與方法

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

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

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

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與分析

2.1 Small RNA測序數(shù)據(jù)質(zhì)控及數(shù)據(jù)評估

根據(jù)得到的測序數(shù)據(jù)顯示，敏感蟲株和耐藥蟲株的樣本數(shù)據(jù)過濾掉低質(zhì)量序列和去接頭后，分別獲得14 916 974和13 888 536個clean reads，占reads總數(shù)的91.83%和94.62%，表明測序質(zhì)量較高，數(shù)據(jù)準(zhǔn)確可靠。將HC-2和HC-XA-aBZ序列注釋到Rfam數(shù)據(jù)庫進(jìn)行比對，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 敏感蟲株和耐藥蟲株比對上Rfam中非編碼RNA的tag豐度統(tǒng)計Table 1 Comparison on tag abundance statistics of non-coding RNA in Rfam in sensitive and resistant strains

2.2 miRNA的鑒定及分類注釋

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

圖1 捻轉(zhuǎn)血矛線蟲丙硫咪唑敏感蟲株和耐藥蟲株Small RNA長度分布統(tǒng)計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)計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á)分析

通過對HC-2和Hc-XA-aBZ鑒定得到的miRNA進(jìn)行匯總，運用基于負(fù)二項分布的DEseq2，根據(jù)read counts計算每個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個差異顯著的miRNA，其中上調(diào)的miRNA有113個，下調(diào)的miRNA有181個(圖3)。

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

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

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

將miRNA預(yù)測到的差異靶基因進(jìn)行KEGG富集注釋，結(jié)果表明DEGs中富集到了胞吞作用(Endocytosis)，碳代謝(Carbon metabolism)，內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工(Protein processing in endoplasmic reticulum)，吞噬體(Phagosome)等327條通路上，這些通路中包括FoxO信號通路(FoxO signaling pathway)，mTOR信號通路(mTOR signaling pathway)，PI3K-Akt信號通路(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)中，miRNA與其所預(yù)測到的靶基因之間大部分都是以負(fù)相關(guān)的方式進(jìn)行調(diào)控。在本研究當(dāng)中，將預(yù)測到的miRNA和mRNA進(jìn)行關(guān)聯(lián)分析，共有1 770個miRNA-mRNA對(共涉及到274個差異的miRNA和603個差異的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驗證差異表達(dá)的miRNA

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

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

3 討 論

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

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

Wang等[25]在對KLFs家族基因的研究中發(fā)現(xiàn)，KLF2的異位表達(dá)抑制了胃癌細(xì)胞的增殖、遷移和侵襲，在過表達(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提供了很有價值的信息。

4 結(jié) 論

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