陳昕迪 王騰宇 劉春霞 王文龍*

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

捻轉(zhuǎn)血矛線蟲是一種寄生于反芻動(dòng)物消化道內(nèi)強(qiáng)致病性寄生蟲，可引起寄生宿主貧血和慢性消耗性臨床癥狀，嚴(yán)重可導(dǎo)致病畜死亡[1]。目前，捻轉(zhuǎn)血矛線蟲病的防治主要依靠伊維菌素和阿苯達(dá)唑等化學(xué)藥物，Dahourou等[2]對(duì)農(nóng)場(chǎng)進(jìn)行流行病學(xué)調(diào)查發(fā)現(xiàn)，在布基納法索地區(qū)的傳統(tǒng)養(yǎng)殖場(chǎng)中，綿羊胃腸道寄生蟲的總體流行率為91%。早在2000年，國(guó)際獸醫(yī)局(International epizootic office, OIE)公布的全球調(diào)查數(shù)據(jù)顯示，不科學(xué)用藥已造成全球超過(guò)80%的國(guó)家存在捻轉(zhuǎn)血矛線蟲多重耐藥的流行情況，給畜牧業(yè)造成了巨大的經(jīng)濟(jì)損失[3]。但捻轉(zhuǎn)血矛線蟲產(chǎn)生耐藥性的調(diào)控機(jī)制還處于初探階段，為了防止耐藥性的產(chǎn)生并指導(dǎo)養(yǎng)殖戶合理用藥，挖掘耐藥性候選調(diào)控因子至關(guān)重要。

隨著近年來(lái)lncRNA在生命科學(xué)領(lǐng)域被高度關(guān)注，目前在人類和小鼠等模式生物中l(wèi)ncRNA的研究已較為深入[4]，長(zhǎng)鏈非編碼RNA(Long non-coding RNA，lncRNA)是一類長(zhǎng)度大于200 nt且很少具備編碼能力的RNA[5]，廣泛參與生物體的生長(zhǎng)發(fā)育、免疫反應(yīng)以及細(xì)胞凋亡與分化等生物過(guò)程[6]。李騰等[7]在對(duì)微小隱孢子蟲的研究中發(fā)現(xiàn)，當(dāng)微小隱孢子蟲感染人結(jié)直腸癌細(xì)胞后，lncRNA NEAT1的表達(dá)顯著上調(diào)，并通過(guò)調(diào)控細(xì)胞IL-8的轉(zhuǎn)錄,參與宿主細(xì)胞抗微小隱孢子蟲免疫反應(yīng)。王振勛等[8]通RNA-蛋白沉降、Western blotting等檢測(cè)手段，對(duì)弓形蟲慢性感染的小鼠腦組織樣本研究并得出結(jié)論，lncRNA 102796通過(guò)與靶基因oprd1結(jié)合抑制小膠質(zhì)細(xì)胞的增殖，影響細(xì)胞周期，進(jìn)而導(dǎo)致神經(jīng)細(xì)胞的損傷。Gutschner等[9]在對(duì)lncRNA與癌癥發(fā)生的相關(guān)研究中發(fā)現(xiàn)，一些lncRNA表現(xiàn)出明顯的進(jìn)化保守性和調(diào)控能力，通過(guò)表達(dá)水平上的差異和某些特異性lncRNA超表達(dá)來(lái)影響腫瘤細(xì)胞增殖，逃避生長(zhǎng)抑制，進(jìn)而在癌變過(guò)程中發(fā)揮關(guān)鍵作用，此類lncRNA可作為癌癥早期診斷及預(yù)后分子標(biāo)記，并為后期藥物的開(kāi)發(fā)提供新靶點(diǎn)。Zhou等[10]通過(guò)對(duì)乳腺癌細(xì)胞mTOR通路水平的檢測(cè)發(fā)現(xiàn)，lncRNA CCAT2可以通過(guò)激活mTOR信號(hào)通路降低乳腺癌細(xì)胞對(duì)5-氟尿嘧啶的耐藥性，進(jìn)一步揭示lncRNA在解決疾病耐藥性問(wèn)題中的重要作用。隨著新發(fā)現(xiàn)的lncRNA呈現(xiàn)指數(shù)型增長(zhǎng)，人類疾病中大量耐藥性相關(guān)lncRNA被發(fā)掘出來(lái)，但在寄生蟲的耐藥性相關(guān)研究中，非編碼RNA的研究仍較為遲緩，尤其是捻轉(zhuǎn)血矛線蟲lncRNA的研究還相對(duì)滯后。

鑒于此，本研究通過(guò)對(duì)捻轉(zhuǎn)血矛線蟲伊維菌素敏感蟲株和耐藥蟲株的轉(zhuǎn)錄組測(cè)序，比較轉(zhuǎn)錄本間差異水平并根據(jù)編碼能力對(duì)lncRNA表達(dá)譜進(jìn)行預(yù)測(cè)和篩選，以期從lncRNA分子層面探究耐藥蟲體在伊維菌素藥物作用下的抵抗機(jī)制，為進(jìn)一步闡釋捻轉(zhuǎn)血矛線蟲伊維菌素耐藥調(diào)節(jié)機(jī)制及有效干擾靶點(diǎn)的開(kāi)發(fā)提供新思路。

1 材料與方法

1.1 試驗(yàn)蟲株

捻轉(zhuǎn)血矛線蟲伊維菌素耐藥蟲株采自內(nèi)蒙古自治區(qū)烏蘭浩特市科爾沁右翼前旗，通過(guò)對(duì)捻轉(zhuǎn)血矛線蟲感染強(qiáng)度EPG≥800的試驗(yàn)羊群進(jìn)行兩倍劑量的二次給藥后，經(jīng)糞便蟲卵檢查和對(duì)照試驗(yàn)調(diào)查捻轉(zhuǎn)血矛線蟲感染情況，同時(shí)采用蟲卵減少試驗(yàn)進(jìn)行復(fù)檢測(cè)定感染強(qiáng)度后確定感染蟲株為伊維菌素耐藥蟲株[11]。敏感蟲株通過(guò)幼蟲發(fā)育抑制試驗(yàn)測(cè)定敏感蟲株伊維菌素LD50敏感閾值低于9 ng/mL，且正常劑量伊維菌素體外驅(qū)蟲后的糞便蟲卵減少率大于90%，進(jìn)而確定為捻轉(zhuǎn)血矛線蟲伊維菌素敏感蟲株[12]。對(duì)選取的綿羊進(jìn)行剖檢，采集捻轉(zhuǎn)血矛線蟲成蟲并逐條鑒別雌雄，分裝至凍存管并儲(chǔ)存于液氮罐中。

1.2 成蟲cDNA文庫(kù)構(gòu)建及高通量測(cè)序

分別選取捻轉(zhuǎn)血矛線蟲伊維菌素敏感蟲株和耐藥蟲株雌蟲，雄蟲各1管，每管50條，并將雌雄蟲混合，使用Trizol(Invitrogen，美國(guó))提取蟲體的Total RNA，并對(duì)RNA純度及濃度進(jìn)行分析。采用去核糖體試劑盒(New england biolabs)對(duì)提取出的Total RNA進(jìn)行過(guò)濾去除，并將RNA隨機(jī)打斷成短片段，再以片斷化后的RNA為模板，用六堿基隨機(jī)引物(Random hexamers)合成cDNA第一鏈和第二鏈。經(jīng)過(guò)Qiaquick PCR(Qiagen)試劑盒純化洗脫并進(jìn)行末端修復(fù)、加堿基A，加測(cè)序接頭，降解第二條鏈后進(jìn)行凝膠電泳檢測(cè)。將構(gòu)建好的文庫(kù)使用Illumina HiSeqTM4 000平臺(tái)進(jìn)行雙末端測(cè)序。

1.3 轉(zhuǎn)錄組數(shù)據(jù)質(zhì)控及l(fā)ncRNA識(shí)別獲取

利用Fastp[13]軟件對(duì)Raw reads進(jìn)行質(zhì)控，過(guò)濾去除低質(zhì)量Reads。使用Bowtie2[14]將較短的Reads比對(duì)GeneBank中的核糖體數(shù)據(jù)庫(kù)，去除比對(duì)上核糖體的Reads。采用HISAT2[15]軟件將保留下來(lái)的Unmapped reads比對(duì)到參考基因組(GenBank：GCA_000469685.2)，并與RNA-seq測(cè)序數(shù)據(jù)進(jìn)行全局和局部的比對(duì)。

依據(jù)比對(duì)全部定位到基因組上的Reads(Total_mapped reads)的結(jié)果，計(jì)算參考基因組中的分布位置。將比對(duì)到基因組的區(qū)域分為外顯子區(qū)(Exon)、內(nèi)含子區(qū)(Intron)和基因間區(qū)(Intergenic)，以 Reads 在參考基因上的分布情況評(píng)價(jià) mRNA 打斷的隨機(jī)程度。使用Stringtie進(jìn)行重構(gòu)轉(zhuǎn)錄本，保留轉(zhuǎn)錄本長(zhǎng)度≥200 bp且外顯子數(shù)目≥2的Reads。

1.4 lncRNA鑒定及差異表達(dá)分析

通過(guò)CPC2[16]和CNCI[17]2個(gè)軟件對(duì)新轉(zhuǎn)錄本的編碼能力進(jìn)行分析，取交集作為候選lncRNAs。采用Stringtie軟件校正轉(zhuǎn)錄本長(zhǎng)度，計(jì)算lncRNA的FPKM(Fragments per kilobase of transcript per million mapped reads)值。使用EdgeR[18]對(duì)Reads count進(jìn)行標(biāo)準(zhǔn)化分析，計(jì)算假設(shè)檢驗(yàn)概率(P-value)并校正FDR值，基于差異分析結(jié)果篩選FDR<0.05且|log2FC|>1的lncRNA為顯著差異lncRNA。

1.5 差異表達(dá)lncRNA的靶基因預(yù)測(cè)及功能注釋

在lncRNA的順式調(diào)控中，位于上游的lncRNA可能與啟動(dòng)子或者共表達(dá)基因的其他順式調(diào)控元件有交集，進(jìn)而在轉(zhuǎn)錄或者轉(zhuǎn)錄后對(duì)基因的表達(dá)進(jìn)行調(diào)控，利用RNAplex[19]尋找2個(gè)長(zhǎng)鏈RNA之間的相互作用，預(yù)測(cè)反義lncRNA與mRNA之間的互補(bǔ)結(jié)合，同時(shí)使用該程序中的ViennaRNA包[20]，根據(jù)熱力學(xué)結(jié)構(gòu)計(jì)算最小自由能來(lái)預(yù)測(cè)最佳堿基配對(duì)關(guān)系。

將lncRNA的靶基因應(yīng)用KOBAS軟件向GO數(shù)據(jù)庫(kù)映射，同時(shí)進(jìn)行靶基因的KEGG富集分析，顯著富集的Pathway閾值為Q value≤0.05，注釋富集結(jié)果以分類二級(jí)柱狀圖和圈圖的形式展現(xiàn)。

1.6 lncRNA的ceRNA調(diào)控網(wǎng)絡(luò)構(gòu)建及連通性分析

采用Target Finder軟件將FDR<0.05且|log2FC|>1 的lncRNA和P<0.05且|log2FC|>1的miRNA依據(jù)ceRNA調(diào)控網(wǎng)絡(luò)關(guān)系，構(gòu)建lncRNA-miRNA調(diào)控網(wǎng)絡(luò)，篩選ceRNA中關(guān)鍵miRNA反應(yīng)元件(MicroRNA response elements, MRE)。

1.7 耐藥相關(guān)通路中的ceRNA調(diào)控網(wǎng)絡(luò)及連通性分析

根據(jù)Lam等[21]的研究方法對(duì)lncRNA和mRNA進(jìn)行關(guān)聯(lián)分析，結(jié)合miRNA-lncRNA-mRNA三者間的調(diào)控關(guān)系使用Cytoscape構(gòu)建ceRNA調(diào)控網(wǎng)絡(luò)。依據(jù)lncRNA、mRNA、circRNA與其存在靶向調(diào)控關(guān)系的miRNA分子數(shù)量進(jìn)行連通性分析，并對(duì)連通網(wǎng)絡(luò)進(jìn)行可視化。根據(jù)KEGG數(shù)據(jù)庫(kù)富集結(jié)果構(gòu)建Pathway-Gene網(wǎng)絡(luò)，利用Igraph對(duì)耐藥相關(guān)通路中篩選出的lncRNA-miRNA-mRNA進(jìn)行可視化分析。

1.8 差異表達(dá)lncRNA的熒光定量PCR驗(yàn)證

為保證測(cè)序數(shù)據(jù)的可靠性，從篩選出的DE lncRNA中隨機(jī)挑選10個(gè)lncRNA進(jìn)行qRT-PCR驗(yàn)證，每組樣品設(shè)置3個(gè)重復(fù)，引物序列見(jiàn)表1。選取GAPDH作為內(nèi)參基因，用Primer 6.0設(shè)計(jì)lncRNA引物序列并送至Sangon合成，采用2-ΔΔCt法計(jì)算相對(duì)表達(dá)量。

表1 lncRNA熒光定量相關(guān)引物Table 1 lncRNA fluorescence quantitative primers

2 結(jié)果與分析

2.1 RNA-seq數(shù)據(jù)質(zhì)量評(píng)估

捻轉(zhuǎn)血矛線蟲伊維菌素敏感株和耐藥株提取總RNA后，得到各樣本濃度均大于450 ng/μL，總量大于3 μg且無(wú)DNase等雜質(zhì)污染和降解，OD260/OD280的比值為1.8～2.1。電泳結(jié)果可觀察到在18 S和28 S處有條帶且在凝膠樣圖譜的相應(yīng)位置正常起峰，表明RNA完整度良好，符合建庫(kù)上機(jī)和qRT-PCR的質(zhì)量要求(圖1)。敏感株和耐藥株分別獲得12.99和11.99 G的原始Raw data，兩組數(shù)據(jù)的Q20占比分別為96.54%和97.72%，Q30占比分別為90.56%和93.44%。在進(jìn)行深度過(guò)濾后分別得到12.71 G和11.92 G的Clean data，其中比對(duì)上參考基因組唯一位置的Reads分別為73.62%和70.09%，Reads比對(duì)基因各部分均勻，且堿基組成平衡，質(zhì)量較高(圖2)，進(jìn)一步保證后續(xù)分析的準(zhǔn)確性。

圖1 捻轉(zhuǎn)血矛線蟲敏感(a)和耐藥(b)RNA凝膠樣譜圖和電泳圖Fig.1 H. contortus sensitive (a) and resistant (b) RNA gel-like spectra and electropherograms

圖2 捻轉(zhuǎn)血矛線蟲敏感蟲株(a)和耐藥蟲株(b)堿基分布圖Fig.2 Base distribution map of H. contortus sensitive (a) and resistant strains (b)

2.2 lncRNA及其靶基因的預(yù)測(cè)與鑒定

基于Stringtie重構(gòu)的5 374條轉(zhuǎn)錄本，使用CPC2和CNCI 2個(gè)軟件分別預(yù)測(cè)到1 141條和1 180 條lncRNA，共靶向到4 377個(gè)mRNA。兩款軟件預(yù)測(cè)到的沒(méi)有編碼潛能的轉(zhuǎn)錄本交集篩選出858個(gè)新lncRNA(圖3)，根據(jù)ceRNA競(jìng)爭(zhēng)調(diào)控機(jī)制共靶向到474個(gè)mRNA。結(jié)合新預(yù)測(cè)到的lncRNA在基因組上的位置，統(tǒng)計(jì)到基因間長(zhǎng)鏈非編碼RNA(Intergenic lncRNAs)400個(gè)、雙向長(zhǎng)鏈非編碼RNA(Bidirectional lncRNAs)21個(gè)、內(nèi)含子長(zhǎng)鏈非編碼RNA(Intronic lncRNAs)11個(gè)、反義長(zhǎng)鏈非編碼RNA(Antisense lncRNAs)63個(gè)和正義長(zhǎng)鏈非編碼RNA(Sense overlapping lncRNAs)245個(gè)。

圖3 捻轉(zhuǎn)血矛線蟲蟲體lncRNA預(yù)測(cè)韋恩圖Fig.3 Venn diagram of lncRNA prediction of H. contortus

2.3 顯著差異表達(dá)的lncRNA篩選

根據(jù)RNA-seq獲得的數(shù)據(jù)，將所有比對(duì)預(yù)測(cè)到的lncRNA進(jìn)行FPKM值和Count值統(tǒng)計(jì)后發(fā)現(xiàn)，在敏感蟲株和耐藥蟲株中分別有65.85%和81.95%的lncRNA的FPKM值≥1(圖4)，說(shuō)明lncRNA在捻轉(zhuǎn)血矛線蟲耐藥性的產(chǎn)生過(guò)程中呈現(xiàn)中高豐度表達(dá)且處于上升趨勢(shì)，且在所有鑒定到的lncRNA中MSTRG.9400.1的表達(dá)量最高。

圖4 捻轉(zhuǎn)血矛線蟲敏感株和耐藥株蟲體 lncRNA表達(dá)豐度小提琴圖Fig.4 Violin plot of lncRNA expression abundance ofH. contortus sensitive and resistant strains

伊維菌素敏感蟲株和耐藥蟲株預(yù)測(cè)到的858個(gè)新lnRNA均存在不同程度的差異表達(dá)，將FDR<0.05且|log2FC|≥1作為進(jìn)一步顯著性篩選閾值，共獲得205個(gè)顯著差異表達(dá)的lncRNA，其中102個(gè)lncRNA上調(diào)，最顯著上調(diào)的3個(gè)lncRNA分別是MSTRG.13692.1、MSTRG.8350.1和MSTRG.3950.2；103個(gè)lncRNA下調(diào)，其中表達(dá)量最顯著下調(diào)3個(gè)lncRNA是MSTRG.10930.1、MSTRG.7170.1和MSTRG.3121.1(圖5)。

圖5 捻轉(zhuǎn)血矛線蟲敏感株和抗藥株之間 顯著差異表達(dá)lncRNA火山圖Fig.5 Volcano map of lncRNA differentially expressed between H. contortus sensitive and resistant strains

2.4 差異lncRNA的順式調(diào)控基因富集分析

將lncRNA與mRNA進(jìn)行關(guān)聯(lián)分析，其中723個(gè)lncRNA通過(guò)順式調(diào)控作用對(duì)1 320個(gè)上游和下游基因進(jìn)行潛在調(diào)控，共匹配到了1 566個(gè)調(diào)控關(guān)系。將lncRNA通過(guò)順式調(diào)控的基因向GO數(shù)據(jù)庫(kù)映射，得到每個(gè)GO功能富集情況的基因列表，應(yīng)用超幾何檢驗(yàn)，獲得顯著富集的GO條目。敏感蟲株和耐藥蟲株的GO富集結(jié)果顯示：lncRNA所靶向到的基因在生物學(xué)過(guò)程中注釋到了代謝過(guò)程(GO:0008152)95個(gè)，細(xì)胞過(guò)程(GO:0009987)89個(gè)，對(duì)刺激的反應(yīng)(GO:0050896)24個(gè)共計(jì)18條GO terms，表明靶基因在受到藥物的刺激下廣泛參與到了各類化合物的代謝；在細(xì)胞組分中注釋到了細(xì)胞膜結(jié)構(gòu)(GO:0016020)47個(gè)，細(xì)胞器(GO:0043226)45個(gè)和細(xì)胞部分(GO:0005623)61個(gè)等9條GO terms(圖6)，細(xì)胞膜等膜結(jié)構(gòu)被大量注釋，說(shuō)明膜組織在捻轉(zhuǎn)血矛線蟲耐藥性的產(chǎn)生過(guò)程中發(fā)揮重要作用。在分子功能中被大量注釋到催化活性(GO:0003824)106個(gè)，各類結(jié)合(GO:0005488)79個(gè)，轉(zhuǎn)運(yùn)活動(dòng)(GO:0005215)14個(gè)等8條GO terms；在分子功能注釋到的一些lncRNA靶基因中，很多基因具備離子通道催化活性劑、轉(zhuǎn)移酶活性等分子功能，大量的功能性物質(zhì)被轉(zhuǎn)運(yùn)，廣泛參與藥物脅迫應(yīng)激和化合物代謝。

圖6 lncRNA cis調(diào)控的上下游基因GO功能注釋Fig.6 GO functional annotation of lncRNA cis-regulated upstream and downstream genes

KEGG通路富集顯示，預(yù)測(cè)得到的lncRNA所靶向的上下游基因共注釋到291條通路。其中注釋基因數(shù)最多的前10位分別為代謝途徑(KO01100)、光傳導(dǎo)(KO04745)、黏著力(KO04510)、肌動(dòng)蛋白細(xì)胞骨架的調(diào)節(jié)(KO04810)、內(nèi)吞作用(KO04144)、血小板活化(KO04611)、卵母細(xì)胞減數(shù)分裂(KO04114)、甘油磷脂代謝(KO00564)、黏附連接(KO04520)、多巴胺能突觸(KO04728)和孕酮介導(dǎo)的卵母細(xì)胞成熟(KO04914)，由KEGG數(shù)據(jù)庫(kù)通路富集圈圖(圖7)可以看出代謝相關(guān)通路被大量富集，其中細(xì)胞色素P450及對(duì)異生物質(zhì)的代謝，視黃醇代謝和甘油磷脂代謝等耐藥相關(guān)代謝通路也被大量富集，由此進(jìn)一步闡述在耐藥發(fā)生過(guò)程中捻轉(zhuǎn)血矛線蟲產(chǎn)生較大的代謝功能障礙。

圖7 lncRNA cis調(diào)控的上下游基因KEGG數(shù)據(jù)庫(kù)富集圈圖Fig.7 KEGG database enrichment circle map of lncRNA cis-regulated upstream and downstream genes

2.5 lncRNA的ceRNA調(diào)控網(wǎng)絡(luò)構(gòu)建

根據(jù)ceRNA競(jìng)爭(zhēng)性調(diào)節(jié)機(jī)制，蛋白編碼基因mRNA、長(zhǎng)鏈非編碼RNA(lncRNA)、假基因(Pseudogenes)、環(huán)狀RNA(circRNA)均是通過(guò)共同的MREs結(jié)合調(diào)控相同的miRNA，由此構(gòu)建lncRNA-miRNA調(diào)控網(wǎng)絡(luò)(圖8)共預(yù)測(cè)到205個(gè)lncRNA靶向結(jié)合371個(gè)miRNA，組成3 227組靶向調(diào)控關(guān)系組合。同一個(gè)miRNA可以靶向結(jié)合多個(gè)lncRNA，在眾多的調(diào)控組合中novel-m0251-3p能夠結(jié)合14個(gè)lncRNA且數(shù)量最多。miR-203-y、miR-2184-x、miR-219-z、miR-32-x和miR-92-y均能結(jié)合到12個(gè)lncRNA，具有較強(qiáng)的結(jié)合能力，但有約74.25%的miRNA靶向到的lncRNA小于6個(gè)。lncRNA中MSTRG.11468.4、MSTRG.11468.5、MSTRG.11468.8和MSTRG.11468.9可以靶向結(jié)合到的miRNA數(shù)量較多，分別為322 422和22個(gè)miRNA，其中89.54%的lncRNA能夠結(jié)合到miRNA的數(shù)量均小于10個(gè)。

圖8 捻轉(zhuǎn)血矛線蟲敏感蟲株和耐藥蟲株lncRNA-miRNA調(diào)控網(wǎng)絡(luò)Fig.8 lncRNA-miRNA regulatory network of H. contortus sensitive and resistant strains

2.6 耐藥相關(guān)調(diào)控信號(hào)通路的ceRNA調(diào)控網(wǎng)絡(luò)

通過(guò)對(duì)伊維菌素敏感蟲株和耐藥蟲株差異mRNA與miRNA，miRNA與lncRNA的靶向調(diào)控關(guān)系預(yù)測(cè)結(jié)果進(jìn)行l(wèi)ncRNA-circRNA-miRNA-mRNA連通性分析(圖9)，結(jié)果顯示在連通性前5的特定連通性節(jié)點(diǎn)上分別有118、62、33、22和9條調(diào)控網(wǎng)絡(luò)，且MSTRG.11468、MSTRG.14365.2和MSTRG.5594.2與circRNA-miRNA-mRNA網(wǎng)絡(luò)具有較強(qiáng)的連通性，進(jìn)一步表明其具有更強(qiáng)的潛在調(diào)控能力。通過(guò)KEGG富集結(jié)果對(duì)FoxO signaling pathway、cAMP signaling pathway、MAPK signaling pathway、AMPK signaling pathway和PI3K-Akt signaling pathway等耐藥相關(guān)信號(hào)通路與連通性高的節(jié)點(diǎn)網(wǎng)絡(luò)進(jìn)行關(guān)聯(lián)分析，同時(shí)篩選出部分假定蛋白(表2)，根據(jù)上述分析結(jié)果構(gòu)建的耐藥性相關(guān)通路lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)(圖10)可以看出lncRNA所靶向到的mRNA單一性較高，結(jié)合miRNA多靶效應(yīng)具有較強(qiáng)指向性的特點(diǎn)進(jìn)一步篩選出ceRNA調(diào)控網(wǎng)絡(luò)中的關(guān)鍵調(diào)節(jié)因子。

圖9 連通性前10的ceRNA連通性桑基圖Fig.9 Sankey diagram of the top ten ceRNA connectivity

表2 捻轉(zhuǎn)血矛線蟲敏感蟲株與耐藥蟲株差異lncRNA序列互補(bǔ)靶標(biāo)的功能注釋Table 2 Functional annotations of complementary targets with sequences of lncRNA differentially expressed between H. contortus sensitive and resistant strains

圖10 耐藥相關(guān)通路中顯著差異基因的lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)關(guān)系Fig.10 The lncRNA-miRNA-mRNA regulatory network relationship of significantly differentially expressed genes in drug resistance-related pathways

2.7 lncRNA的表達(dá)驗(yàn)證

利用qRT-PCR技術(shù)，在伊維菌素敏感蟲株和耐藥蟲株比較組中隨機(jī)挑選10個(gè)(上)下調(diào)表達(dá)量差異水平較大的DE lncRNA進(jìn)行表達(dá)趨勢(shì)的驗(yàn)證。如圖11所示，10個(gè)DE lncRNA在比較組中的變化情況均通過(guò)qRT-PCR結(jié)果在組內(nèi)真實(shí)表達(dá)出來(lái)，且與RNA-seq測(cè)序結(jié)果趨勢(shì)一致，進(jìn)一步證實(shí)轉(zhuǎn)錄組測(cè)序結(jié)果的穩(wěn)定性和可靠性。

圖11 差異表達(dá)lncRNA qRT-PCR與RNA-seq比較分析Fig.11 Comparative analysis between qRT-PCR and RNA-seq of differentially expressed lncRNAs

3 討 論

反芻動(dòng)物感染消化道線蟲特別是捻轉(zhuǎn)血矛線蟲給畜牧業(yè)生產(chǎn)造成較為嚴(yán)重的影響，由于無(wú)法及時(shí)準(zhǔn)確篩查出已感染的羊群，大多數(shù)養(yǎng)殖戶選擇大劑量群體給藥進(jìn)行驅(qū)蟲和預(yù)防[11]，致使羊群對(duì)伊維菌素和阿苯達(dá)唑等藥物產(chǎn)生耐藥性，此現(xiàn)象受到國(guó)際社會(huì)的廣泛關(guān)注[22]。相較于秀麗隱桿線蟲等模式生物，捻轉(zhuǎn)血矛線蟲的組學(xué)數(shù)據(jù)較為匱乏且研究相對(duì)滯后。目前，lncRNA已被證實(shí)在擬南芥(Arabidopsisthaliana)、蜜蜂球囊菌(Ascosphaeraapis)和羅伯茨綠僵菌(Metarhiziumrobertsii)[23-24]等真核生物中均參與到代謝和免疫等過(guò)程并調(diào)控特定蛋白的表達(dá)水平。有研究證實(shí)耐藥性的產(chǎn)生可能與藥物代謝的增加、藥物靶點(diǎn)的突變以及細(xì)胞凋亡和自噬存在密切聯(lián)系，同時(shí)信號(hào)通路的異常激活也可能導(dǎo)致耐藥強(qiáng)度的改變。非編碼 RNA 能夠極大地增加調(diào)控的復(fù)雜性并起到微調(diào)轉(zhuǎn)錄組的作用，并能快速調(diào)整蛋白質(zhì)組以響應(yīng)刺激[25]。

本研究以伊維菌素敏感蟲株和耐藥蟲株為研究對(duì)象，通過(guò)高通量測(cè)序技術(shù)將獲得的Raw data進(jìn)行重新組裝，經(jīng)預(yù)測(cè)發(fā)現(xiàn)捻轉(zhuǎn)血矛線蟲中不存在已知和已存在的lncRNA，根據(jù)序列在基因組的位置信息，長(zhǎng)度大小和編碼能力鑒定到858個(gè)新lncRNA，可以與474個(gè)mRNA靶向結(jié)合，鑒定到的lncRNA以|log2FC|≥1且FDR<0.05為顯著性閾值，存在205個(gè)顯著差異lncRNA。本研究通過(guò)表達(dá)量分析發(fā)現(xiàn)，預(yù)測(cè)到的所有新lncRNA在敏感蟲株與耐藥蟲株中均存在表達(dá)，但顯著差異的lncRNA數(shù)目相較于總數(shù)目而言占比較少，且高表達(dá)基因差異水平未呈現(xiàn)斷崖式變化，由此推測(cè)捻轉(zhuǎn)血矛線蟲在對(duì)伊維菌素產(chǎn)生抗性過(guò)程中僅通過(guò)較少數(shù)量的調(diào)控網(wǎng)絡(luò)便可對(duì)抗性水平進(jìn)行調(diào)節(jié)。依據(jù)GC含量高的基因組間區(qū)存在更多的編碼基因的假說(shuō)[26]，對(duì)轉(zhuǎn)錄本數(shù)據(jù)中獲得的大量lncRNA座進(jìn)行預(yù)測(cè)發(fā)現(xiàn)，lncRNA中的GC含量明顯低于mRNA，本研究結(jié)果與該假說(shuō)猜想一致。Azlan等[27]通過(guò)對(duì)埃及伊蚊(Aedesaegypti)lncRNA表達(dá)譜鑒定發(fā)現(xiàn)，lncRNA呈現(xiàn)低表達(dá)、低 GC 含量、長(zhǎng)度短和低保守性等特點(diǎn)，在本研究所有預(yù)測(cè)到的新lncRNA中，有262個(gè)lncRNA的FPKM值小于1，其余大部分均呈現(xiàn)出中高豐度表達(dá)，但相較于mRNA的表達(dá)情況而言，lncRNA依舊處于低表達(dá)水平，與前人研究結(jié)果一致，同時(shí)也說(shuō)明lncRNA的表達(dá)情況往往與鄰近的反義蛋白編碼基因相關(guān)。在lncRNA-mRNA的調(diào)控網(wǎng)絡(luò)中有356組調(diào)控關(guān)系是上調(diào)lncRNA抑制mRNA下調(diào)表達(dá)，378組下調(diào)lncRNA從而上調(diào)mRNA表達(dá)，這種反向調(diào)控模式占所有l(wèi)ncRNA-mRNA調(diào)控網(wǎng)絡(luò)的46.9%，lncRNA與編碼基因間通過(guò)正向和反向調(diào)控構(gòu)成了遠(yuǎn)超本身數(shù)量的復(fù)雜機(jī)體微調(diào)網(wǎng)絡(luò)，再次說(shuō)明同一共有l(wèi)ncRNA廣泛參與到捻轉(zhuǎn)血矛線蟲伊維菌素耐藥性的產(chǎn)生過(guò)程。

KEGG富集結(jié)果顯示，被富集到較多的通路是代謝和轉(zhuǎn)運(yùn)相關(guān)通路，其中顯著差異基因約有三分之二富集到代謝過(guò)程。Siddiqui等[28]發(fā)現(xiàn)在大多數(shù)給藥后或耐藥生物的通路富集結(jié)果中其他藥物代謝酶、藥物代謝細(xì)胞色素 P450、谷胱甘肽代謝、甘油磷脂代謝、糖酵解/糖異生和異生物質(zhì)代謝等通路均存在富集情況，這些代謝途徑與信號(hào)通路存在緊密的聯(lián)系，不同的信號(hào)通路都發(fā)揮著不同的調(diào)控作用。在腫瘤細(xì)胞的研究中耐藥細(xì)胞與藥物敏感細(xì)胞相比，許多脂質(zhì)都發(fā)生了顯著變化[29]，本研究中甘油磷脂代謝等相關(guān)脂質(zhì)通路被顯著富集，甘油磷脂代謝作為重要的藥物代謝通路已被證實(shí)在蠕蟲中藥物會(huì)通過(guò)氨基酸和甘油磷脂產(chǎn)生的小分子進(jìn)行代謝[30]。根據(jù)功能預(yù)測(cè)分析，甘油磷脂代謝通路中的 F09G2.8作為一類磷脂水解酶酶和磷脂酶D病毒包膜結(jié)構(gòu)域蛋白在MSTRG.14103.3的下調(diào)影響下低豐度表達(dá)，由此推測(cè)MSTRG.14103.3可靶向調(diào)節(jié)脂類小分子蛋白參與捻轉(zhuǎn)血矛線蟲伊維菌素的代謝過(guò)程從而影響耐藥性。通過(guò)對(duì)lncRNA富集情況進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)cAMP信號(hào)通路和MAPK信號(hào)通路同時(shí)富集到MSTRG.3141.1及其靶基因PRKACA。Moody等[31]發(fā)現(xiàn)除了MAPK和PI3K-AKt等信號(hào)通路外，PRKACA在曲妥珠單抗耐藥的乳腺癌樣本中上調(diào)表達(dá)，并在人表皮生長(zhǎng)因子受體2的治療中拯救了細(xì)胞，從而使細(xì)胞對(duì)拉帕替尼治療敏感，其中PRKACA的表達(dá)情況與本研究中的結(jié)果一致。在人類肺癌對(duì)鉑類藥物的抗藥性研究中表明，PRKACA作為一類cAMP依賴性蛋白，ATP轉(zhuǎn)化為的環(huán)AMP可以與cAMP依賴蛋白結(jié)合調(diào)節(jié)細(xì)胞轉(zhuǎn)運(yùn)激酶、鳥(niǎo)嘌呤交換因子和感覺(jué)信號(hào)向環(huán)苷酸門控通道轉(zhuǎn)導(dǎo)信號(hào)，進(jìn)而提高細(xì)胞對(duì)藥物的敏感性[32]。同樣被注釋到的MAPK信號(hào)通路可以將細(xì)胞外信號(hào)傳導(dǎo)至細(xì)胞核內(nèi)，在脅迫應(yīng)答和細(xì)胞增殖與凋亡等方面起到關(guān)鍵作用[33]，由此可以推測(cè)在伊維菌素耐藥性產(chǎn)生過(guò)程中MSTRG.3141.1的下調(diào)可能激活下游轉(zhuǎn)錄因子與cAMP依賴性蛋白啟動(dòng)子的近端區(qū)域相互作用，并將信號(hào)傳遞到細(xì)胞內(nèi)以促進(jìn)靶基因PRKACA表達(dá)。

ABC轉(zhuǎn)運(yùn)蛋白是一類專一性外排轉(zhuǎn)運(yùn)蛋白，通過(guò)跨細(xì)胞膜轉(zhuǎn)運(yùn)分子逆濃度梯度調(diào)節(jié)內(nèi)源性和外源性物質(zhì)的細(xì)胞水平，并參與多種生物過(guò)程。ABC轉(zhuǎn)運(yùn)蛋白的過(guò)度表達(dá)與大環(huán)內(nèi)酯類藥物的抗性存在密切聯(lián)系[34]。P-糖蛋白(Pgp)是ABC轉(zhuǎn)運(yùn)蛋白超家族中重要的功能蛋白，在人類的多種疾病和寄生蟲的耐藥性研究中均出現(xiàn)明顯異常表達(dá)，是重要的耐藥相關(guān)蛋白[35-36]。在捻轉(zhuǎn)血矛線蟲耐伊維菌素的研究中發(fā)現(xiàn)，Pgp基因均出現(xiàn)表達(dá)情況的異常改變及SNP位點(diǎn)的突變[37]，本研究中ABC轉(zhuǎn)運(yùn)途徑中富集到Pgp-1和ced-7等多個(gè)基因，但通過(guò)ceRNA結(jié)合能力預(yù)測(cè)Pgp-1并未靶向到lncRNA和miRNA，而ced-7卻靶向到大量miRNA，由此推測(cè)機(jī)體可能通過(guò)ABC轉(zhuǎn)運(yùn)蛋白家族中的其他轉(zhuǎn)運(yùn)分子進(jìn)行ceRNA網(wǎng)絡(luò)調(diào)節(jié)進(jìn)而達(dá)到藥物逃逸的作用。lncRNA作為研究起步較晚的一類非編碼RNA，在捻轉(zhuǎn)血矛線蟲伊維菌素耐藥相關(guān)的研究中尚未見(jiàn)相關(guān)報(bào)道，由于各物種間lncRNA序列的保守性較差，無(wú)法確定lncRNA的轉(zhuǎn)錄確切位置。本研究通過(guò)編碼潛能預(yù)測(cè)獲取的lncRNA表達(dá)譜并與靶基因和miRNA進(jìn)行ceRNA網(wǎng)絡(luò)構(gòu)建，盡可能去除無(wú)效的非編碼RNA分子。對(duì)于篩選出的耐藥相關(guān)且特異性表達(dá)的lncRNA及其調(diào)控組合，后續(xù)還需對(duì)蟲體進(jìn)行干擾和過(guò)表達(dá)等方式進(jìn)行生物學(xué)功能驗(yàn)證，以期進(jìn)一步闡述非編碼RNA在捻轉(zhuǎn)血矛線蟲耐藥性產(chǎn)生過(guò)程中發(fā)揮的重要作用。

4 結(jié) 論

通過(guò)RNA-seq技術(shù)，獲取到捻轉(zhuǎn)血矛線蟲伊維菌素敏感蟲株和耐藥蟲株lncRNA表達(dá)譜并進(jìn)行全面分析，解析在耐藥性產(chǎn)生過(guò)程中l(wèi)ncRNA的數(shù)量、種類以及表達(dá)譜間存在的明顯差異，揭示部分lncRNA可通過(guò)順式調(diào)控作用參與ceRNA競(jìng)爭(zhēng)調(diào)節(jié)機(jī)制，并廣泛參與到藥物代謝和轉(zhuǎn)運(yùn)等生物學(xué)過(guò)程。研究結(jié)果可為耐藥性關(guān)鍵調(diào)控分子的篩選和功能性研究提供重要參考，也為深度開(kāi)發(fā)捻轉(zhuǎn)血矛線蟲耐藥性關(guān)鍵lncRNA及其潛在ceRNA調(diào)控元件提供數(shù)據(jù)支持。