許雅靜， 孫明會， 劉佳美， 郭意龍， 胡 穎， 張佳欣，趙 蕭， 張凱遙， 熊翠玲， 陳大福,2， 郭 睿,2

(1.福建農(nóng)林大學動物科學學院 蜂學學院， 福州 350002； 2.福建農(nóng)林大學蜂療研究所， 福州 350002)

1 引 言

意大利蜜蜂(Apismelliferaligustica)作為西方蜜蜂(Apismellifera)的亞種之一，在我國和其他國家的養(yǎng)蜂生產(chǎn)中廣泛使用，具有重要的生態(tài)和經(jīng)濟價值[1].蜜蜂是一種真社會性昆蟲，因群居性而易遭受細菌、真菌和病毒等病原微生物的侵襲[1].東方蜜蜂微孢子蟲(Nosemaceranae)是一種專性侵染成年蜜蜂中腸上皮細胞的真菌病原，對蜂王、工蜂、雄蜂及蜜蜂幼蟲均具有侵染性，導(dǎo)致的蜜蜂微孢子蟲病嚴重危害蜜蜂健康并影響?zhàn)B蜂業(yè)的可持續(xù)發(fā)展[2].東方蜜蜂微孢子蟲侵染可對蜜蜂宿主造成諸多不良影響，包括能量脅迫、免疫抑制、細胞凋亡抑制、中腸上皮細胞結(jié)構(gòu)破壞及壽命縮短等[3].

非編碼RNA(non-coding RNA, ncRNA)起初被認為不具有蛋白編碼能力和生物學功能，是基因轉(zhuǎn)錄過程的“噪音”和“副產(chǎn)物”，但隨著技術(shù)手段的不斷進步和相關(guān)研究的持續(xù)深入，ncRNA被證實在許多生物學過程中起到關(guān)鍵作用[4].根據(jù)長度，ncRNA可分為長鏈非編碼RNA(long non-coding RNA, lncRNA)和小ncRNA，后者又包括微小RNA(microRNA, miRNA)、小干擾RNA(Small interfering RNA, siRNA)和PIWI蛋白互作RNA(PIWI-interacting RNA, piRNA)等.成熟piRNA的長度一般介于26～30 nt，與miRNA(20～24 nt)和siRNA(21～25 nt)的長度相近[5].在真核生物中，piRNA能參與控制轉(zhuǎn)座子的活性及維持基因組的穩(wěn)定性[6].此外，有研究表明piRNA通過與靶mRNA的堿基互補配對發(fā)揮基因表達調(diào)控作用[7].PIWI蛋白最早在果蠅(Drosophilamelanogaster)體內(nèi)被發(fā)現(xiàn)，并被證實參與了果蠅種系干細胞的維持和自我更新[8].隨著深度測序技術(shù)的興起與發(fā)展，人們在哺乳動物[9]、昆蟲[10]和魚類[11]等物種中均陸續(xù)鑒定出piRNA.與人和果蠅等少數(shù)模式生物相比，昆蟲的piRNA研究較為滯后，迄今僅在果蠅[12]、白紋伊蚊(Aedesalbopictus)[13]和蜜蜂[5]等少數(shù)昆蟲中有過相關(guān)報道.Jones等人[14]研究表明piRNA途徑可能在維持果蠅基因組完整性、組織功能、代謝穩(wěn)態(tài)和免疫防御方面起到一定作用.Guo等通過過表達和敲減piRNA-3312發(fā)現(xiàn)piRNA-3312可靶向結(jié)合腸道酯酶1基因進而負調(diào)控淡色庫蚊Culexpipienspallens對擬除蟲菊酯的抗性[15].

目前，關(guān)于piRNA參與蜜蜂宿主與病原互作的研究未見報道；東方蜜蜂微孢子蟲脅迫下意蜂工蜂的piRNA應(yīng)答研究仍然缺失.本研究擬基于已獲得的small RNA(sRNA)組學數(shù)據(jù)對東方蜜蜂微孢子蟲脅迫下意蜂工蜂中腸的DEpiRNA進行篩選、分析和驗證，進而對DEpiRNA進行靶向預(yù)測及分析，以期解析意蜂工蜂中腸響應(yīng)東方蜜蜂微孢子蟲脅迫的piRNA差異表達譜及調(diào)控網(wǎng)絡(luò)，并揭示DEpiRNA調(diào)控宿主應(yīng)答的潛在作用.

2 材料與方法

2.1 材 料

2.1.1 供試蜜蜂及微孢子蟲 意蜂工蜂取自福建農(nóng)林大學動物科學學院(蜂學學院)教學蜂場的實驗蜂群.東方蜜蜂微孢子蟲孢子由福建農(nóng)林大學動物科學學院(蜂學學院)蜜蜂保護實驗室純化和保存[16].

2.1.2 sRNA組學數(shù)據(jù)來源 前期研究中，筆者所在團隊已利用small RNA-seq(sRNA-seq)技術(shù)對東方蜜蜂微孢子蟲脅迫7 d和10 d的意蜂工蜂中腸(AmT1和AmT2)及未受脅迫的工蜂中腸(AmCK1和AmCK2)進行測序，共測得16 589 574條原始讀段(raw reads)，過濾后得到132 028 990條有效讀段(clean reads)，各組中clean reads占比均達到77.04%及以上；AmCK1、AmT1、AmCK2和AmT2組內(nèi)3個不同生物學重復(fù)之間的Pearson相關(guān)性系數(shù)分別達到99.93%、87.92%、99.67%和91.14%[17].高質(zhì)量的sRNA-seq組學數(shù)據(jù)可為本研究中piRNA的差異表達及調(diào)控網(wǎng)絡(luò)分析提供可靠的數(shù)據(jù)基礎(chǔ).

2.2 方 法

2.2.1 DEpiRNA的篩選與分析 根據(jù)edgeR軟件(http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)對AmCK1 vs AmT1和AmCK2 vs AmT2比較組中的DEpiRNA進行篩選，篩選標準為|log2(Fold change)|≥1(|log2FC|≥1)且P≤0.05，其余采用軟件默認參數(shù).利用基迪奧在線云平臺(www.omicshare.com)的相關(guān)工具對上述兩個比較組中的DEpiRNA進行Venn分析，采用默認參數(shù).

2.2.2 DEpiRNA的靶mRNA預(yù)測與分析 參照杜宇等[18]的方法，利用TargetFinder軟件[19]對DEpiRNA進行靶向預(yù)測.通過BLAST工具將靶mRNA比對到GO數(shù)據(jù)庫(https://www.geneontology.org)和KEGG數(shù)據(jù)庫(https://www.genome.jp/kegg/)以獲得相應(yīng)的功能和通路注釋.最后利用Cytoscape軟件[20]分別對上述兩個比較組中的DEpiRNA及其靶mRNA、共有DEpiRNA及其靶mRNA的調(diào)控網(wǎng)絡(luò)進行可視化.均采用默認參數(shù).

2.2.3 DEpiRNA的Stem-loop RT-PCR與qPCR驗證 利用FastPure?Cell/Tissue Total RNA Isolation Kit V2(Vazyme公司，中國)分別提取意蜂7日齡和10日齡工蜂中腸樣品的總RNA.利用超微量分光光度儀(Thermo Fisher，美國)測定RNA質(zhì)量.從上述兩個比較組中隨機選取4個DEpiRNA(piR-ame-967173、piR-ame-459976、piR-ame-779990、piR-ame-1044429)進行驗證.根據(jù)上述piRNA的核酸序列設(shè)計相應(yīng)的特異性Stem-loop引物和上游引物(F)以及通用下游引物(R)，委托上海生物工程有限公司合成引物(附表1).參照HiScript?1st Strand cDNA Synthesis Kit試劑盒說明書利用Stem-loop引物反轉(zhuǎn)錄合成cDNA作為piRNA的PCR模板；利用隨機引物和oligo (dT)引物進行常規(guī)反轉(zhuǎn)錄，合成的cDNA作為內(nèi)參基因snRNA U6的PCR模板.PCR體系(20 μL)包括：10 μL PCR mix(南京諾唯贊，中國)，1 μL上游引物(F)、1 μL通用下游引物(R)、1 μL cDNA模板以及7 μL DEPC處理水.PCR反應(yīng)條件為：95 ℃ 5 min，95 ℃ 10 s，55 ℃ 30 s，40個循環(huán).產(chǎn)物經(jīng)1.8%的瓊脂糖凝膠電泳檢測.采用凝膠圖像成像分析系統(tǒng)(上海培清，中國)進行凝膠觀察和拍照.

參照筆者所在團隊前期已建立的方法[18]對上述4個DEpiRNA進行qPCR驗證.利用Stem-loop引物進行反轉(zhuǎn)錄得到相應(yīng)的cDNA，作為模板進行qPCR反應(yīng).反應(yīng)按照SYBR Green Dye試劑盒(Vazyme公司，中國)操作說明書進行.反應(yīng)體系(20 μL)包含：SYBR Green Dye 10 μL，cDNA模板1.3 μL，正、反向引物各1 μL，DEPC水6.7 μL.反應(yīng)在ABI QuantStudio 3熒光定量PCR儀(ABI，美國)上進行，程序設(shè)置為：95 ℃預(yù)變性5 min；95 ℃變性10 s，60 ℃退火延伸30 s，共40個循環(huán).每個反應(yīng)進行3次技術(shù)重復(fù).采用2-ΔΔCt法計算DEpiRNA相對表達量.最后利用Graph Prism 7軟件對qPCR數(shù)據(jù)進行student’st檢驗，并計算組間差異顯著性及繪圖，采用默認參數(shù).

3 結(jié)果與分析

3.1 東方蜜蜂微孢子蟲脅迫下意蜂工蜂中腸的piRNA差異表達譜

差異表達分析結(jié)果顯示，AmCK1 vs AmT1比較組包含50個DEpiRNA，上調(diào)piRNA和下調(diào)piRNA的數(shù)量分別為46和4個；其中上調(diào)倍數(shù)最高的是piR-ame-1207641(log2FC=16.81,P=6.37E-9)，其次是piR-ame-500944(log2FC=13.13,P=4.45E-6)和piR-ame-748814(log2FC=12.49,P=5.14E-7)；下調(diào)倍數(shù)最高的是piR-ame-1221425(log2FC=-19.97,P=0.0001)，其次是piR-ame-1182394(log2FC=-19.65,P=0.0001)和piR-ame-1219424(log2FC=-17.26,P=0.0002).AmCK2 vs AmT2比較組包含207個DEpiRNA，上調(diào)piRNA和下調(diào)piRNA的數(shù)量分別為56和151個；其中上調(diào)倍數(shù)最高的是piR-ame-1008436(log2FC=15.01,P=3.11E-7)，其次是piR-ame-771408(log2FC=15.01,P=3.11E-7)和piR-ame-716466(log2FC=14.95,P=5.36E-7)；下調(diào)倍數(shù)最高的是piR-ame-706008(log2FC=-14.46,P=7.49E-6)，其次是piR-ame-634882(log2FC=-14.09,P=9.88E-7)和piR-ame-1253391(log2FC=-13.96,P=3.08E-5).Venn分析結(jié)果顯示，上述兩個比較組的共有DEpiRNA為10個(表1)，特有DEpiRNA分別為40和197個.

表1 東方蜜蜂微孢子蟲脅迫7 d和10 d的意蜂工蜂中腸共有DEpiRNA的詳細信息

3.2 東方蜜蜂微孢子蟲脅迫下意蜂工蜂中腸DEpiRNA的靶mRNA預(yù)測及注釋

靶向預(yù)測結(jié)果顯示，AmCK1 vs AmT1比較組中piR-ame-748814、piR-ame-762269、piR-ame-1128833、piR-ame-175077可分別靶向2035、927、3021和1273條mRNA；AmCK2 vs AmT2比較組中piR-ame-1048676、piR-ame-1055898、piR-ame-11093、piR-ame-1128833、piR-ame-165683、piR-ame-190949、piR-ame-202265、piR-ame-471566、piR-ame-503081、piR-ame-630674、piR-ame-706008、piR-ame-742536、piR-ame-797627、piR-ame-932156可分別靶向1290、631、1195、3021、3549、4040、1152、2583、512、4935、1297、3191、593、1163條mRNA.

GO數(shù)據(jù)庫注釋結(jié)果顯示，AmCK1 vs AmT1比較組中DEpiRNA的靶mRNA可注釋到細胞進程(1381)和代謝進程(1214)等18個生物學進程相關(guān)條目，細胞器(503)和細胞(687)等16個細胞組分相關(guān)條目，結(jié)合(1509)和轉(zhuǎn)運活性(161)等11個分子功能相關(guān)條目(圖1a)；AmCK2 vs AmT2比較組中DEpiRNA的靶mRNA可注釋到細胞進程(3744)和生物學調(diào)控(1087)等19個生物學進程相關(guān)條目，細胞膜組件(1529)和細胞膜(1581)等17個細胞組分相關(guān)條目，信號轉(zhuǎn)導(dǎo)活性(417)和結(jié)構(gòu)分子活性(76)等11個分子功能相關(guān)條目(圖1b).括號內(nèi)的數(shù)字代表注釋在該條目的靶mRNA數(shù)目.

KEGG數(shù)據(jù)庫注釋結(jié)果顯示，AmCK1 vs AmT1比較組中DEpiRNA的靶mRNA共注釋到117條代謝通路，包括Jak-STAT(6)、磷脂酰肌醇信號系統(tǒng)(46)、Wnt信號通路(42)、mTOR信號通路(9)、氨基酸代謝(10)、碳水化合物代謝(14)及脂質(zhì)代謝(10)等(圖2a)；AmCK2 vs AmT2比較組中DEpiRNA的靶mRNA共注釋到140條代謝通路，包括Wnt信號通路(117)、FoxO信號通路(82)、MAPK信號通路(26)、氨基酸代謝(13)、碳水化合物代謝(14)及輔助因子和維生素代謝(12)等(圖2b).括號內(nèi)的數(shù)字代表注釋在該通路的靶mRNA數(shù)目.

3.3 東方蜜蜂微孢子蟲脅迫下意蜂工蜂中腸DEpiRNA及其靶mRNA的調(diào)控網(wǎng)絡(luò)

調(diào)控網(wǎng)絡(luò)分析結(jié)果顯示，AmCK1 vs AmT1比較組中同一個DEpiRNA可靶向多個mRNA，如piR-ame-748814可靶向多達2035個mRNA；此外，部分mRNA可被1～2個DEpiRNA靶向結(jié)合，如NM_001014992.1可被piR-ame-1128833和piR-ame-175077共同靶向結(jié)合，NM_001040233.1可被piR-ame-762269和piR-ame-175077共同靶向結(jié)合(圖3a).類似地，AmCK2 vs AmT2比較組中也存在同一個DEpiRNA靶向多個mRNA、部分mRNA被多個DEpiRNA靶向結(jié)合的現(xiàn)象(圖3b).

圖1 AmCK1 vs AmT1 (a)和AmCK2 vs AmT2 (b)比較組DEpiRNA的靶mRNA的GO數(shù)據(jù)庫注釋Fig.1 GO database annotation of mRNAs targeted by DEpiRNAs in AmCK1 vs AmT1 (a) and AmCK2 vs AmT2 (b) comparison groups

圖2 AmCK1 vs AmT1 (a)和AmCK2 vs AmT2 (b)比較組DEpiRNA的靶mRNA的KEGG數(shù)據(jù)庫注釋

圖3 東方蜜蜂微孢子蟲脅迫7 d和10 d的意蜂工蜂中腸中DEpiRNA及其靶mRNA的調(diào)控網(wǎng)絡(luò)Fig.3 Regulation network of DEpiRNAs and their target mRNAs in A. m. ligustica workers’ midguts at 7 d and 10 d post N.ceranae stress

進一步分析發(fā)現(xiàn)，上述兩個比較組共有的DEpiRNA piR-ame-1128833靶向的mRNA中有30個涉及Hippo、Notch、FoxO、Wnt和Hedgehog等發(fā)育相關(guān)的信號通路，有54個涉及內(nèi)吞作用、MAPK信號通路、Jak-STAT信號通路及泛素介導(dǎo)的蛋白水解等免疫途徑.

3.4 DEpiRNA的Stem-loop RT-PCR與RT-qPCR驗證

利用Stem-loop RT-PCR對隨機選取4個DEpiRNA進行擴增，產(chǎn)物的電泳結(jié)果顯示均能擴增出符合預(yù)期大小(約80 bp)的目的片段(圖4).進一步的qPCR結(jié)果顯示上述4個DEpiRNA的表達趨勢與測序數(shù)據(jù)中的趨勢一致.以上結(jié)果證實了本研究中piRNA表達趨勢的真實性和測序數(shù)據(jù)的可靠性.

4 討 論

本研究選擇東方蜜蜂微孢子蟲脅迫7 d和10 d兩個時間點的意蜂工蜂中腸主要基于三點依據(jù)：(1) 東方蜜蜂微孢子蟲孢子接種西方蜜蜂工蜂后，病原孢子載量在宿主細胞可持續(xù)增長到 20 d[21]；(2) 東方蜜蜂微孢子蟲孢子接種意蜂工蜂后，接種組工蜂的累積死亡率在1～10 dpi階段均高于未接種組工蜂，但兩組累計死亡率只在7 dpi和10 dpi兩個時間點具有顯著性差異[22]；(3) 筆者團隊前期

圖4 DEpiRNA的Stem-loop RT-PCR(a)和qPCR驗證(b～e)Fig.4 Stem-loop RT-PCR (a) and qPCR(b～e) validation of DEpiRNAs

已基于東方蜜蜂微孢子蟲接種7 d和10 d的意蜂工蜂中腸及未接種工蜂中腸轉(zhuǎn)錄組數(shù)據(jù)解析了宿主響應(yīng)脅迫的免疫基因應(yīng)答[23]，分析和探討了miRNA和lncRNA在宿主的脅迫應(yīng)答中的調(diào)控作用[17,22]，本研究是在前期基礎(chǔ)上進一步探討東方蜜蜂微孢子蟲脅迫下意蜂工蜂的piRNA差異表達譜及潛在功能，以期全面深入地闡釋意蜂工蜂響應(yīng)東方蜜蜂微孢子蟲脅迫的作用機制.相比于未受脅迫的意蜂工蜂中腸，我們在東方蜜蜂微孢子蟲脅迫7 d和10 d的工蜂中腸中分別鑒定到50和207個DEpiRNA，說明這些DEpiRNA參與了宿主的脅迫應(yīng)答.AmCK1 vs AmT1和AmCK2 vs AmT2比較組共有的DEpiRNA為10個，特有的DEpiRNA分別為40和197個，推測上述共有DEpiRNA在宿主應(yīng)答過程的2個時間點均發(fā)揮作用，而多數(shù)DEpiRNA在宿主應(yīng)答的不同時間點發(fā)揮特定作用.

piRNA參與了許多后生動物的生殖系發(fā)育和配子發(fā)生且與人類的癌癥高度相關(guān)[24].piRNA也被證實能通過靶向mRNA發(fā)揮基因表達調(diào)控的作用[25].本研究中，AmCK1 vs AmT1中僅有4個DEpiRNA可靶向7256條mRNA，AmCK2 vs AmT2中有14個DEpiRNA可靶向29152條mRNA，說明少數(shù)DEpiRNA參與宿主應(yīng)答過程中基因表達調(diào)控，而多數(shù)DEpiRNA可能主要發(fā)揮保持基因組穩(wěn)定性的作用.AmCK1 vs AmT1中DEpiRNA的靶mRNA可注釋到細胞進程等45個功能條目(圖1a)；還注釋到磷脂酰肌醇信號系統(tǒng)等117條代謝通路(圖2a).類似地，AmCK2 vs AmT2中DEpiRNA的靶mRNA可注釋到細胞進程等47個功能條目(圖1b)以及MAPK信號通路等140條代謝通路(圖2b).以上結(jié)果表明DEpiRNA通過調(diào)節(jié)諸多生物學過程潛在影響宿主對東方蜜蜂微孢子蟲脅迫的應(yīng)答.miRNA作為基因表達的關(guān)鍵調(diào)控因子，其“種子區(qū)”根據(jù)堿基互補配對原則靶向結(jié)合mRNA的3’UTR特定區(qū)域，同一個miRNA能夠靶向多個mRNA，反之亦然[26].本研究也發(fā)現(xiàn)，同一個DEpiRNA可靶向結(jié)合多個mRNA，一些mRNA可同時被多個DEpiRNA靶向結(jié)合，說明DEpiRNA能以一種類似于miRNA的方式參與宿主應(yīng)答中的基因表達調(diào)控[27].

信號轉(zhuǎn)導(dǎo)途徑在真核生物的環(huán)境適應(yīng)、細胞活動、生長發(fā)育、新陳代謝和免疫應(yīng)答等生物學過程中發(fā)揮舉足輕重的作用[28].Hippo信號通路通過調(diào)節(jié)細胞增殖和凋亡參與器官大小控制和組織再生[29].Notch信號通路可調(diào)節(jié)果蠅中腸上皮細胞分化及細胞凋亡[30].FoxO信號通路參與調(diào)控碳水化合物代謝和能量代謝等重要代謝過程[31].同時，F(xiàn)oxO和Wnt信號通路在調(diào)節(jié)細胞生長、增殖和凋亡及新陳代謝等方面均起到關(guān)鍵作用[32-34].Hedgehog信號通路與細胞的生長和發(fā)育密切相關(guān)[35].上述兩個比較組共有的DEpiRNA僅piR-ame-1128833可靶向3021條mRNA，暗示piR-ame-1128833在意蜂工蜂響應(yīng)東方蜜蜂微孢子蟲脅迫的過程中發(fā)揮重要作用.本研究發(fā)現(xiàn)，piR-ame-1128833的靶mRNA中分別有11、1、2、11和5個可注釋到Hippo、Notch、FoxO、Wnt和Hedgehog等信號轉(zhuǎn)導(dǎo)途徑，表明piR-ame-1128833通過調(diào)控生長、發(fā)育和物質(zhì)能量代謝通路潛在影響意蜂工蜂對東方蜜蜂微孢子蟲脅迫的應(yīng)答.

為抵御病原體和寄生蟲的入侵，昆蟲進化出了高效的先天免疫系統(tǒng)，包括以內(nèi)吞和吞噬作用為代表的細胞免疫[36]和以Jak-STAT和MAPK等信號通路及抗菌肽釋放為代表的體液免疫[37].媒介昆蟲的中腸腸壁是病毒入侵的主要屏障，研究表明依賴于網(wǎng)格蛋白的內(nèi)吞作用在促進病毒進入煙粉虱(Bemisiatabaci)中腸細胞方面起到重要作用[38].MAPK信號通路可以通過級聯(lián)反應(yīng)調(diào)控中腸基因表達進而引起小菜蛾(Plutellaxylostella)的抗病毒反應(yīng)[39].JAK/STAT信號通路參與調(diào)節(jié)細胞生長、分化、凋亡及炎癥反應(yīng)，并與埃及伊蚊(Aedesaegypti)抵御真菌侵染密切相關(guān)[40].泛素介導(dǎo)的蛋白水解系統(tǒng)作為消除錯誤折疊或損傷細胞的主要機制，是許多生理過程如信號轉(zhuǎn)導(dǎo)、細胞周期及免疫細胞功能的控制器[41].本研究發(fā)現(xiàn)，piR-ame-1128833的靶mRNA中分別有32、2、1和19個可注釋到內(nèi)吞作用、MAPK信號通路、Jak-STAT通路及泛素介導(dǎo)的蛋白水解系統(tǒng)，說明piR-ame-1128833可能通過調(diào)節(jié)上述細胞和體液免疫途徑參與宿主對東方蜜蜂微孢子蟲脅迫的應(yīng)答，但仍需要進一步實驗驗證.目前，參照miRNA的過表達和敲減方法對真核生物組織或細胞中的piRNA進行功能研究已見諸報道[42,43].下一步我們將根據(jù)piR-ame-1128833的序列設(shè)計合成相應(yīng)的模擬物和抑制物，并通過飼喂法對東方蜜蜂微孢子蟲感染的意蜂工蜂內(nèi)的piR-ame-1128833進行過表達和敲減，以明確該piRNA在宿主的脅迫應(yīng)答中的功能.

綜上，東方蜜蜂微孢子蟲脅迫可引起意蜂工蜂中腸piRNA表達譜的改變；DEpiRNA可通過靶向mRNA潛在調(diào)控宿主響應(yīng)脅迫中的多個生物學過程；piR-ame-1128833可通過靶向調(diào)控相關(guān)mRNA潛在調(diào)節(jié)Hippo和Wnt等信號通路及內(nèi)吞作用和MAPK信號通路等免疫途徑.