意大利蜜蜂工蜂中腸發(fā)育過(guò)程中長(zhǎng)鏈非編碼RNA的差異表達(dá)分析

郭睿，耿四海，熊翠玲，鄭燕珍，付中民，王海朋，杜宇，童新宇，趙紅霞，陳大福

?

意大利蜜蜂工蜂中腸發(fā)育過(guò)程中長(zhǎng)鏈非編碼RNA的差異表達(dá)分析

郭睿1，耿四海1，熊翠玲1，鄭燕珍1，付中民1，王海朋1，杜宇1，童新宇1，趙紅霞2，陳大福1

（1福建農(nóng)林大學(xué)蜂學(xué)學(xué)院，福州 350002；2廣東省生物資源應(yīng)用研究所，廣州 510260）

【目的】長(zhǎng)鏈非編碼RNA（lncRNA）在真核生物的基因表達(dá)、表觀遺傳和細(xì)胞周期調(diào)控等方面發(fā)揮重要功能。本研究旨在探究意大利蜜蜂（, 簡(jiǎn)稱意蜂）工蜂中腸發(fā)育過(guò)程中l(wèi)ncRNA的表達(dá)譜及其作用。【方法】利用RNA-seq技術(shù)和鏈特異性建庫(kù)方法對(duì)意蜂7和10日齡工蜂中腸（Am7、Am10）進(jìn)行深度測(cè)序，下機(jī)的原始數(shù)據(jù)經(jīng)過(guò)Perl腳本過(guò)濾得到高質(zhì)量有效讀段。利用bowtie工具將有效讀段比對(duì)核糖體數(shù)據(jù)庫(kù)，進(jìn)一步利用TopHat2軟件將未比對(duì)到核糖體數(shù)據(jù)庫(kù)上的數(shù)據(jù)比對(duì)到參考基因組。利用CPC和CNCI軟件對(duì)轉(zhuǎn)錄本的編碼能力進(jìn)行預(yù)測(cè)。通過(guò)RT-PCR對(duì)部分lncRNA進(jìn)行鑒定。利用edgeR軟件進(jìn)行差異表達(dá)lncRNA（DElncRNA）分析，進(jìn)而預(yù)測(cè)lncRNA的上下游基因，并對(duì)上下游基因進(jìn)行GO及KEGG代謝通路富集分析。聯(lián)用RNAhybrid、Miranda和TargetScan軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA及miRNA靶向結(jié)合的靶基因，并通過(guò)Cytoscape軟件構(gòu)建DElncRNAs-miRNAs-mRNAs的調(diào)控網(wǎng)絡(luò)。最后，通過(guò)RT-qPCR驗(yàn)證測(cè)序數(shù)據(jù)的可靠性?！窘Y(jié)果】Am7和Am10的深度測(cè)序分別獲得134 802 058和147 051 470條原始讀段，經(jīng)嚴(yán)格過(guò)濾分別得到134 166 157和146 293 288條有效讀段；共得到3 890個(gè)DElncRNA，包括2 005個(gè)上調(diào)lncRNA與1 885個(gè)下調(diào)lncRNA。RT-PCR驗(yàn)證結(jié)果顯示共有8個(gè)lncRNA能擴(kuò)增出符合預(yù)期的目的片段，表明預(yù)測(cè)出的lncRNA真實(shí)存在。DElncRNA的上下游基因數(shù)為1 793個(gè)，它們涉及42個(gè)GO條目，包括代謝進(jìn)程、發(fā)育進(jìn)程、細(xì)胞進(jìn)程、應(yīng)激反應(yīng)和免疫系統(tǒng)進(jìn)程等；這些上下游基因還涉及251條代謝通路，包括碳代謝、嘌呤代謝和脂肪酸的生物合成等物質(zhì)代謝通路，硫代謝、甲烷代謝和氧化磷酸化等能量代謝通路，Hippo信號(hào)通路、Wnt信號(hào)通路和Notch信號(hào)通路等信號(hào)通路，溶酶體、內(nèi)吞作用和泛素介導(dǎo)的蛋白水解等細(xì)胞免疫通路，以及MAPK信號(hào)通路、Jak-STAT信號(hào)通路和NF-kappa B信號(hào)通路等體液免疫通路，上述結(jié)果表明DElncRNA在意蜂中腸發(fā)育過(guò)程中參與物質(zhì)和能量代謝、細(xì)胞生命活動(dòng)和免疫調(diào)控。進(jìn)一步分析發(fā)現(xiàn)TCONS_00020918可通過(guò)調(diào)控西方蜜蜂王漿主蛋白1編碼基因在意蜂工蜂中腸的營(yíng)養(yǎng)吸收、級(jí)型分化中發(fā)揮功能。DElncRNA的調(diào)控網(wǎng)絡(luò)分析結(jié)果顯示DElncRNA與miRNA、mRNA間存在復(fù)雜的調(diào)控關(guān)系，部分DElncRNA處于調(diào)控網(wǎng)絡(luò)的中心位置且能靶向結(jié)合較多的miRNA，也有部分miRNA可被多個(gè)DElncRNA共同靶向，表明這些DElncRNA可能在中腸發(fā)育中發(fā)揮重要作用。隨機(jī)挑取5個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證，結(jié)果顯示它們的表達(dá)量變化趨勢(shì)與測(cè)序結(jié)果一致，證實(shí)了本研究測(cè)序數(shù)據(jù)的可靠性?！窘Y(jié)論】差異表達(dá)長(zhǎng)鏈非編碼RNA（DElncRNA）廣泛參與意蜂工蜂中腸的新陳代謝、細(xì)胞活動(dòng)和免疫調(diào)控并作為競(jìng)爭(zhēng)性內(nèi)源RNA（ceRNA）發(fā)揮作用，研究結(jié)果為關(guān)鍵lncRNA的篩選和功能研究提供了必要的數(shù)據(jù)支持。

意大利蜜蜂；中腸；發(fā)育；長(zhǎng)鏈非編碼RNA；上下游基因

0 引言

【研究意義】蜜蜂作為社會(huì)學(xué)模式昆蟲(chóng)，在發(fā)育學(xué)、行為學(xué)和神經(jīng)生物學(xué)等方面具有重要價(jià)值[1]。同時(shí)，蜜蜂也是最重要的授粉昆蟲(chóng)，具有不可替代的經(jīng)濟(jì)價(jià)值和生態(tài)價(jià)值[2]。意大利蜜蜂（，簡(jiǎn)稱意蜂）具有優(yōu)越的采集能力、造脾能力和分泌蜂王漿能力[3]，在世界各地的養(yǎng)蜂生產(chǎn)中廣泛使用。目前，有關(guān)蜜蜂中腸發(fā)育機(jī)理及調(diào)控機(jī)制的研究極為滯后，非編碼RNA（non-coding RNA，ncRNA）在中腸發(fā)育過(guò)程中的作用的相關(guān)信息也十分有限。在轉(zhuǎn)錄組水平對(duì)長(zhǎng)鏈非編碼RNA（long non-coding RNA，lncRNA）在意蜂工蜂中腸發(fā)育中的作用進(jìn)行研究，可為深入解析意蜂工蜂中腸發(fā)育的分子機(jī)理提供新的思路和線索?！厩叭搜芯窟M(jìn)展】隨著測(cè)序技術(shù)的發(fā)展，起初被認(rèn)為基因轉(zhuǎn)錄“噪音”的ncRNA成為近幾年的研究熱點(diǎn)。在人類基因組中，編碼蛋白序列所占比例不到2%，超過(guò)98%的序列是不編碼蛋白質(zhì)的[4]。ncRNA根據(jù)其長(zhǎng)短和形狀分為小RNA、lncRNA和環(huán)狀RNA等。其中，lncRNA已被證明在劑量補(bǔ)償效應(yīng)[5]、表觀遺傳[6]、細(xì)胞周期[7]和細(xì)胞分化調(diào)控[8]等眾多生命活動(dòng)中發(fā)揮重要作用，因此成為生命科學(xué)領(lǐng)域的研究熱點(diǎn)。lncRNA的作用方式多樣，包括作為信號(hào)分子調(diào)控上下游基因轉(zhuǎn)錄[9]，作為誘導(dǎo)分子充當(dāng)分子阻斷劑[10]，作為引導(dǎo)分子與蛋白結(jié)合[11]，作為支架分子發(fā)揮“腳手架”作用[12]，作為小RNA的前體[13]，以及作為競(jìng)爭(zhēng)內(nèi)源RNA結(jié)合miRNA對(duì)mRNA起剪接調(diào)控作用[14]。近年來(lái)，有關(guān)lncRNA在腫瘤發(fā)生機(jī)制方面作用的研究取得了顯著進(jìn)展，在脂肪代謝、肌肉發(fā)育及免疫抗病等機(jī)理方面作用的研究也取得了許多突破性進(jìn)展[15-16]。然而，lncRNA在昆蟲(chóng)領(lǐng)域的相關(guān)研究較為滯后，僅在果蠅和家蠶等極少數(shù)模式昆蟲(chóng)中有較少的lncRNA相關(guān)研究報(bào)道[17]。此前，有研究報(bào)道lncRNA在蜜蜂的級(jí)型分化[18]、卵巢發(fā)育[19-20]和抵抗病毒入侵[21]中起重要的調(diào)控作用，但總體上有關(guān)蜜蜂lncRNA的研究仍非常有限。腸道是昆蟲(chóng)的重要消化器官和免疫器官。蜜蜂腸道的相關(guān)研究主要集中在腸道微生物方面[22-27]，如Engel等[26]利用二代測(cè)序技術(shù)對(duì)意蜂工蜂后腸進(jìn)行宏基因組測(cè)序，基于數(shù)據(jù)分析發(fā)現(xiàn)意蜂后腸內(nèi)的少量細(xì)菌種類存在極大的遺傳多樣性，進(jìn)一步的比較分析發(fā)現(xiàn)不同的細(xì)菌種類涉及不同的功能，包括宿主互作、菌膜形成及碳水化合物水解等。【本研究切入點(diǎn)】蜜蜂中腸發(fā)育機(jī)理及調(diào)控機(jī)制仍不明確，關(guān)于lncRNA在蜜蜂中腸發(fā)育過(guò)程中作用的研究未見(jiàn)報(bào)道?！緮M解決的關(guān)鍵問(wèn)題】結(jié)合RNA-seq技術(shù)和鏈特異性cDNA建庫(kù)方法對(duì)意蜂7和10日齡的工蜂中腸進(jìn)行測(cè)序，對(duì)中腸發(fā)育過(guò)程的lncRNA進(jìn)行差異表達(dá)分析，進(jìn)一步通過(guò)上下游基因分析和lncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)分析DElncRNA在意蜂工蜂中腸發(fā)育過(guò)程中的作用。

1 材料與方法

試驗(yàn)于2017年在福建農(nóng)林大學(xué)蜂學(xué)學(xué)院蜜蜂保護(hù)實(shí)驗(yàn)室進(jìn)行。

1.1 生物材料

供試意大利蜜蜂工蜂取自福建農(nóng)林大學(xué)蜂學(xué)學(xué)院教學(xué)蜂場(chǎng)。

1.2 樣品制備與Illumina測(cè)序

從群勢(shì)較強(qiáng)的意蜂蜂群中選取有封蓋子的巢脾，快速提至實(shí)驗(yàn)室，放入（34±0.5）℃培養(yǎng)箱中培養(yǎng)至工蜂出房，將剛出房的工蜂（記為0 d）放入干凈的四周打孔以通風(fēng)的塑料盒中，盒子上方插入一支裝有50%（w/v）無(wú)菌糖水的飼喂器（圖1），（34±0.5）℃條件下恒溫培養(yǎng)，每日檢查工蜂存活情況，及時(shí)清理死亡的工蜂。進(jìn)行3次生物學(xué)重復(fù)。

圖1 意蜂工蜂的人工飼養(yǎng)

選取意蜂7和10日齡工蜂中腸作為測(cè)序材料。分別快速拉取意蜂7日齡工蜂中腸（Am7）和10日齡工蜂中腸（Am10），液氮速凍后置于-80℃超低溫冰箱保存?zhèn)溆?，每只中腸的取樣時(shí)間嚴(yán)格控制在15 s以內(nèi)。首先用RNA抽提試劑盒（AxyPrepTMMultisource Total RNA Miniprep Kit）（TaKaRa公司，日本）抽提蜜蜂中腸樣品的總RNA，為最大限度地保留所有非編碼RNA（ncRNA），去除核糖體RNA后的mRNA和ncRNA用裂解緩沖液隨機(jī)打斷為小片段，作為模板用六堿基隨機(jī)引物、緩沖液、dNTPs、RNase H和DNA polymerase I合成cDNA第二鏈。經(jīng)過(guò)QiaQuick PCR試劑盒（Qiagen公司，德國(guó)）純化并加EB緩沖液洗脫經(jīng)末端修復(fù)、加堿基A，加測(cè)序接頭，然后通過(guò)尿嘧啶-N-糖基化酶（UNG）降解cDNA第二鏈。消化產(chǎn)物經(jīng)瓊脂糖凝膠電泳和PCR擴(kuò)增，PCR產(chǎn)物經(jīng)高碘酸鈉（Sigma 公司，美國(guó)）處理。委托廣州基迪奧生物科技有限公司對(duì)上述cDNA文庫(kù)進(jìn)行雙端（paired-end）測(cè)序，測(cè)序平臺(tái)為Illumina HiSeq 4000。測(cè)序數(shù)據(jù)已上傳NCBI SRA數(shù)據(jù)庫(kù)，BioProject號(hào)：PRJNA406998。

1.3 高通量數(shù)據(jù)分析

1.3.1 LncRNA的生物信息學(xué)預(yù)測(cè) 對(duì)于下機(jī)數(shù)據(jù)，利用Perl腳本去除含有adaptor、未知核苷酸比例＞10%和低質(zhì)量reads，獲得有效讀段（clean reads）。使用短reads比對(duì)工具bowtie[28]分別將樣品Am7-1、Am7-2、Am7-3、Am10-1、Am10-2和Am10-3的clean reads比對(duì)（mapping）到核糖體數(shù)據(jù)庫(kù)（最多允許5個(gè)錯(cuò)配），去除mapping上核糖體的reads，利用TopHat2軟件[29]將保留下來(lái)的數(shù)據(jù)進(jìn)一步比對(duì)西方蜜蜂的參考基因組（Amel_4.5）[30]。利用FPKM（Fragments Per Kilobase of transcript per Million mapped reads）法計(jì)算基因表達(dá)量。利用R軟件（version 2.16.2）計(jì)算各樣品之間的相關(guān)性系數(shù)。

lncRNA與mRNA相比在序列保守型、ORF長(zhǎng)度等特點(diǎn)上具有一定差異，利用軟件CPC[31]和CNCI[32]對(duì)新轉(zhuǎn)錄本的編碼能力進(jìn)行預(yù)測(cè)，選取二者的交集作為可靠的預(yù)測(cè)結(jié)果。

1.3.2 DElncRNA及其上下游基因分析 利用edgeR軟件[33]進(jìn)行l(wèi)ncRNA的差異分析，得到DElncRNA。LncRNA的功能與其編碼基因位置臨近的蛋白編碼基因關(guān)系密切，位于上游的lncRNA可能與啟動(dòng)子或共表達(dá)基因的其他順式作用元件存在交集，從而在轉(zhuǎn)錄或轉(zhuǎn)錄后水平進(jìn)行基因表達(dá)調(diào)控；位于3′ UTR或基因下游的lncRNA可能參與其他調(diào)控作用。因此，對(duì)lncRNA進(jìn)行注釋，如果其位于一個(gè)基因的上游或下游，這些lncRNA有可能與順式作用元件所在區(qū)域有交集，從而參與轉(zhuǎn)錄調(diào)控。利用Omicshare平臺(tái)（http://www.omicshare.com/tools/index.php/）對(duì)DElncRNA的上下游基因進(jìn)行GO（Gene Ontology）分類和KEGG pathway富集分析。

1.3.3 DElncRNA的調(diào)控網(wǎng)絡(luò)構(gòu)建 LncRNA可以競(jìng)爭(zhēng)結(jié)合miRNA，從而減少miRNA結(jié)合mRNA，表明lncRNA可以通過(guò)miRNA調(diào)控mRNA的表達(dá)。利用RNAhybrid[34]、Miranda[35]和TargetScan[36]軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA，以及miRNA結(jié)合的mRNA，根據(jù)上述靶向結(jié)合關(guān)系構(gòu)建lncRNA-miRNA- mRNA的調(diào)控網(wǎng)絡(luò)。利用Cytoscape軟件[37]對(duì)上述調(diào)控網(wǎng)絡(luò)進(jìn)行可視化。

1.4 lncRNA的RT-PCR驗(yàn)證

隨機(jī)選取9個(gè)lncRNA，根據(jù)它們的序列利用DNAMAN軟件設(shè)計(jì)相應(yīng)的特異性引物。委托上海生工生物工程有限公司合成引物。利用RNA抽提試劑盒（TaKaRa，中國(guó)）提取Am7和Am10的總RNA，作為模板進(jìn)行反轉(zhuǎn)錄，得到的cDNA作為模板進(jìn)行PCR。PCR反應(yīng)體系（20 μL）包括cDNA模板1 μL，上游引物1 μL，下游引物1 μL，Mixture 10 μL，無(wú)菌水補(bǔ)至20 μL。PCR程序：94℃預(yù)變性5 min；94℃變性50 s，55℃退火30 s，72℃延伸50 s，共34個(gè)循環(huán)；72℃再延伸10 min。PCR產(chǎn)物經(jīng)1.5%瓊脂糖凝膠電泳和凝膠成像儀（上海培清，中國(guó)）檢測(cè)。

1.5 DElncRNA的實(shí)時(shí)熒光定量PCR（RT-qPCR）驗(yàn)證

為了驗(yàn)證RNA-seq數(shù)據(jù)的可靠性，隨機(jī)選取5個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證。利用DNAMAN軟件設(shè)計(jì)相應(yīng)的特異性引物。委托上海生工生物工程有限公司合成引物。利用RNA抽提試劑盒（TaKaRa，中國(guó)）提取Am7和Am10的總RNA，作為模板進(jìn)行反轉(zhuǎn)錄，得到的cDNA作為模板進(jìn)行PCR。RT-qPCR反應(yīng)按照SYBR Green Dye試劑盒（Vazyme公司，中國(guó)）操作說(shuō)明書(shū)進(jìn)行，每個(gè)反應(yīng)進(jìn)行3次重復(fù)。反應(yīng)體系（20 μL）包含正、反向引物（10.0 μmol·L-1）各1 μL，cDNA模板DNA 1 μL，SYBR Green Dye 10 μL，DEPC水7 μL。qRT-PCR反應(yīng)在ABI 7500熒光定量PCR儀（ABI公司，美國(guó)）上進(jìn)行，反應(yīng)條件：95℃預(yù)變性1 min，95℃變性15 s，60℃延伸30 s，共40個(gè)循環(huán)，最后72℃延伸45 s。利用2-ΔΔCt法對(duì)上述基因的相對(duì)表達(dá)量進(jìn)行計(jì)算。

2 結(jié)果

2.1 數(shù)據(jù)質(zhì)控與評(píng)估

在實(shí)驗(yàn)室條件下對(duì)意蜂工蜂進(jìn)行人工飼養(yǎng)（圖1），Am7和Am10中腸樣品的測(cè)序得到raw reads平均分別為134 802 058和147 051 470條，過(guò)濾后得到clean reads平均分別為134 166 157和146 293 288條。各樣品的平均Q20和Q30分別為97.34%和93.86%（表1）。Am7與Am10的組內(nèi)各生物學(xué)重復(fù)之間的Pearson相關(guān)系數(shù)均在0.97以上，說(shuō)明各樣品的重復(fù)性較好（圖2）。上述結(jié)果說(shuō)明本研究的測(cè)序數(shù)據(jù)質(zhì)量良好，可用于進(jìn)一步分析。

隨機(jī)挑選9個(gè)lncRNA進(jìn)行RT-PCR驗(yàn)證，電泳結(jié)果顯示其中有8個(gè)lncRNA均能擴(kuò)增出符合預(yù)期的目的片段，說(shuō)明本研究預(yù)測(cè)出的絕大多數(shù)的lncRNA真實(shí)存在（圖3）。相關(guān)引物信息詳見(jiàn)表2。進(jìn)一步對(duì)上述8個(gè)驗(yàn)證的lncRNA的上下游基因進(jìn)行預(yù)測(cè)，其中TCONS_00020918可調(diào)控西方蜜蜂王漿主蛋白1編碼基因，TCONS_00021005可調(diào)控西方蜜蜂未知蛋白LOC725PPP，TCONS_00025221可調(diào)控西方蜜蜂未知蛋白LOC726321、LOC100577788和LOC102656594。

橫縱坐標(biāo)均表示基因表達(dá)量（FPKM）The horizontal and vertical coordinates represent gene expression level (FPKM)

表1 RNA-seq數(shù)據(jù)概覽

M: DNA marker; 1: TCONS_00020918; 2: TCONS_00021005; 3: TCONS_00019675; 4: TCONS_00019678; 5: TCONS_00025221; 6: TCONS_00025232; 7: TCONS_00025235; 8: TCONS_00025236

表2 RT-PCR與RT-qPCR引物信息

2.2 DElncRNA及其上下游基因分析

LncRNA的差異表達(dá)分析結(jié)果顯示，在Am7 vs Am10比較組中共有3 890個(gè)DElncRNA，包括2 005個(gè)上調(diào)lncRNA和1885個(gè)下調(diào)lncRNA。lncRNA可通過(guò)影響其編碼基因的上下游基因而發(fā)揮作用。共預(yù)測(cè)出1 698個(gè)DElncRNA的上下游基因，其中上游基因、下游基因和重疊基因分別為485、535和593個(gè)。對(duì)DElncRNA編碼基因的上下游基因進(jìn)行GO分類，結(jié)果顯示它們分布在生物學(xué)進(jìn)程、細(xì)胞組分和分子功能3大類，涉及42個(gè)GO term，富集基因數(shù)最多的前10個(gè)GO term分別為結(jié)合（355 genes）、細(xì)胞進(jìn)程（315 genes）、代謝進(jìn)程（299 genes）、單細(xì)胞進(jìn)程（262 genes）、催化活性（259 genes）、細(xì)胞（147 genes）、細(xì)胞組分（147 genes）、細(xì)胞膜（138 genes）、細(xì)胞膜組分（131 genes）和細(xì)胞器（117 genes）（圖4）。

進(jìn)一步對(duì)DElncRNA編碼基因的上下游基因進(jìn)行KEGG pathway富集分析，結(jié)果顯示它們富集在251個(gè)pathway，富集基因數(shù)最多的前10位pathway是Hippo信號(hào)通路（36 genes）、Wnt信號(hào)通路（19 genes）、癌癥microRNA（16 genes）、癌癥通路（16 genes）、剪接體（15 genes）、PI3K-AKt信號(hào)通路（14 genes）、碳代謝（12 genes）、單純皰疹感染（12 genes）、癌癥蛋白聚糖（12 genes）和神經(jīng)營(yíng)養(yǎng)因子信號(hào)通路（11 genes）（圖5-A）。Hippo信號(hào)通路的概貌如圖5-B所示。

2.3 DElncRNA的調(diào)控網(wǎng)絡(luò)分析

LncRNA可以作為一種競(jìng)爭(zhēng)性內(nèi)源RNA，可通過(guò)結(jié)合miRNA減少靶向結(jié)合mRNA的miRNA數(shù)，從而減輕miRNA對(duì)mRNA的抑制作用[38]。利用軟件預(yù)測(cè)DElncRNA靶向結(jié)合的miRNA，以及miRNA靶向結(jié)合的mRNA，其中靶向結(jié)合miRNA最多的DElncRNA為TCONS_00004891（36 miRNAs）、XR_001704571.1（28 miRNAs）和TCONS_00024939（27 miRNAs），靶向結(jié)合mRNA最多的miRNA為ame-miR-6001-3p（98 mRNAs）、mir-136-y（61 mRNAs）和mir-410-y（38 mRNAs）。進(jìn)一步利用Cytoscape軟件對(duì)DElncRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)進(jìn)行可視化，結(jié)果顯示上調(diào)、下調(diào)lncRNA與miRNA和mRNA間均形成復(fù)雜的調(diào)控網(wǎng)絡(luò)（圖6），共有125個(gè)lncRNA預(yù)測(cè)到靶向結(jié)合的miRNA，其中僅有6個(gè)lncRNA靶向結(jié)合1個(gè)miRNA，其余的119個(gè)lncRNA均靶向結(jié)合2個(gè)或更多個(gè)miRNA。

圖4 DElncRNA上下游基因的GO分類

2.4 DElncRNA的qRT-PCR驗(yàn)證

隨機(jī)挑取6個(gè)DElncRNA進(jìn)行RT-qPCR驗(yàn)證，結(jié)果顯示其中有5個(gè)DElncRNA表達(dá)水平的變化趨勢(shì)和轉(zhuǎn)錄組數(shù)據(jù)中相應(yīng)DElncRNA表達(dá)水平的變化趨勢(shì)一致（圖7），說(shuō)明本研究中的測(cè)序數(shù)據(jù)真實(shí)可靠。相關(guān)引物序列信息詳見(jiàn)表2。對(duì)驗(yàn)證的5個(gè)DElncRNA進(jìn)行上下游基因分析，發(fā)現(xiàn)TCONS_00029692可調(diào)控肥大細(xì)胞脫粒肽編碼基因，TCONS_00038012可調(diào)控未知蛋白LOC100576733編碼基因。

3 討論

此前，蜜蜂腸道的相關(guān)研究集中在腸道微生物方面，有關(guān)蜜蜂腸道的發(fā)育機(jī)理、ncRNA在腸道發(fā)育過(guò)程中作用的研究未見(jiàn)報(bào)道。前人研究結(jié)果表明lncRNA與生物的神經(jīng)發(fā)育、細(xì)胞組織發(fā)育、細(xì)胞周期調(diào)控、干細(xì)胞分化、免疫應(yīng)答和癌癥發(fā)生等過(guò)程關(guān)系密切。近年來(lái)，隨著二代測(cè)序技術(shù)的迅速發(fā)展和廣泛應(yīng)用，越來(lái)越多的lncRNA在動(dòng)物、植物和微生物中被鑒定出來(lái)，但蜜蜂lncRNA的研究較為滯后，相關(guān)信息極為有限。為探究意蜂工蜂中腸的發(fā)育機(jī)理及調(diào)控機(jī)制，本研究結(jié)合RNA-seq技術(shù)和鏈特異性建庫(kù)方法對(duì)意蜂7和10日齡工蜂中腸進(jìn)行測(cè)序，基于高質(zhì)量的測(cè)序數(shù)據(jù)共預(yù)測(cè)出6 353個(gè)lncRNA，RT-PCR驗(yàn)證的擴(kuò)增成功率高達(dá)88.9%，表明本研究預(yù)測(cè)出的絕大多數(shù)lncRNA真實(shí)存在。相對(duì)于非特異性建庫(kù)，特異性建庫(kù)的方法在合成cDNA第二鏈時(shí)，將dTTP替換為dUTP，加上接頭后再用UNG酶降解第二鏈，因而該方法可以確定轉(zhuǎn)錄本來(lái)自正鏈或負(fù)鏈，故能更準(zhǔn)確地獲得基因的結(jié)構(gòu)和基因表達(dá)信息。lncRNA的表達(dá)具有時(shí)間和組織特異性，本研究的測(cè)序?qū)ο笫且夥涔し渲心c組織，它們?cè)谄渌M織中是否表達(dá)以及表達(dá)水平有待于進(jìn)一步研究。本研究中人工飼養(yǎng)工蜂的過(guò)程中僅飼喂糖水，工蜂的腸道內(nèi)會(huì)有少量流體且多聚集在后腸，為排除雜質(zhì)對(duì)測(cè)序結(jié)果的影響，舍棄了后腸，將中腸組織用于測(cè)序。測(cè)得的原始數(shù)據(jù)包含工蜂本身的數(shù)據(jù)和少量的中腸內(nèi)微生物的數(shù)據(jù)，為排除后者的干擾，對(duì)原始數(shù)據(jù)進(jìn)行了嚴(yán)格的過(guò)濾和質(zhì)控，先比對(duì)核糖體數(shù)據(jù)庫(kù)以去除rRNA數(shù)據(jù)，再比對(duì)西方蜜蜂基因組去除未比對(duì)上的數(shù)據(jù)，保留下來(lái)的能比對(duì)上的數(shù)據(jù)用于后續(xù)的生物信息學(xué)分析，考慮到昆蟲(chóng)和微生物的物種親緣關(guān)系較遠(yuǎn)，二者的基因保守性很低，因而比對(duì)上參考基因組的數(shù)據(jù)理論上皆為意蜂工蜂中腸本身的數(shù)據(jù)。

A：Pathway富集分析enrichment analysis；B：Hippo信號(hào)通路概貌General picture of Hippo signaling pathway

A：上調(diào)lncRNA的lncRNA-miRNA-mRNA網(wǎng)絡(luò)lncRNA-miRNA-mRNA networks of up-regulated lncRNAs；B：下調(diào)lncRNA的lncRNA-miRNA-mRNA網(wǎng)絡(luò) lncRNA-miRNA-mRNA networks of down-regulated lncRNAs

A: XR_409562.2; B: TCONS_00011956; C: TCONS_00012589; D: TCONS_00038012; E: TCONS_00029692

LncRNA的功能與其編碼基因坐標(biāo)臨近的蛋白編碼基因相關(guān)，位于上游的lncRNA可能與啟動(dòng)子或共表達(dá)基因的其他順式作用元件有交集，從而在轉(zhuǎn)錄或轉(zhuǎn)錄后水平對(duì)基因的表達(dá)進(jìn)行調(diào)控[39]。本研究中，DElncRNA的上下游基因達(dá)到1 697個(gè)，GO富集分析結(jié)果顯示分別有299個(gè)和4個(gè)上下游基因涉及代謝進(jìn)程和發(fā)育進(jìn)程，有259個(gè)上下游基因涉及氧化磷酸化，表明DElncRNA參與意蜂工蜂中腸的物質(zhì)和能量代謝過(guò)程調(diào)控。本研究發(fā)現(xiàn)，315個(gè)上下游基因涉及細(xì)胞進(jìn)程，表明DElncRNA廣泛參與意蜂工蜂中腸的細(xì)胞生命活動(dòng)。此外，還發(fā)現(xiàn)分別有92個(gè)和3個(gè)上下游基因涉及應(yīng)激反應(yīng)和免疫系統(tǒng)進(jìn)程，暗示相應(yīng)的DElncRNA在意蜂工蜂中腸發(fā)育中發(fā)揮調(diào)控作用。王漿主蛋白（MRJP）為幼蟲(chóng)發(fā)育提供必要的可利用氮元素，在蜜蜂行為中有潛在的調(diào)控作用，還參與決定蜂群中雌性蜂的級(jí)型分化等，其中MRJP1在蜂王漿主蛋白中含量最豐富，所占比例達(dá)31%，占蜂王漿中水溶性蛋白的48%[40]。本研究發(fā)現(xiàn)，TCONS_00020918可調(diào)控西方蜜蜂MRJP1編碼基因，表明該lncRNA可能在意蜂的營(yíng)養(yǎng)吸收、級(jí)型分化中發(fā)揮重要的調(diào)控功能。此外，多個(gè)DElncRNA可能調(diào)控未知蛋白編碼基因，具體調(diào)控關(guān)系的明確有賴于Nr數(shù)據(jù)庫(kù)蛋白功能注釋信息的完善及進(jìn)一步的實(shí)驗(yàn)驗(yàn)證。

蜜蜂腸道是物質(zhì)和能量合成與代謝的主要場(chǎng)所。本研究中，DElncRNA上下游基因的KEGG富集分析結(jié)果顯示，共有221個(gè)上下游基因富集在多達(dá)65個(gè)物質(zhì)代謝通路上，包括甘油酯代謝（5 genes）和脂肪酸的生物合成（2 genes）等10條脂代謝相關(guān)通路，檸檬酸鹽循環(huán)（6 genes）和丙酮酸鹽代謝（2 genes）等13條碳水化合物代謝相關(guān)通路，賴氨酸降解（3 genes）和精氨酸生物合成（2 genes）等13條氨基酸代謝相關(guān)通路，以及嘌呤代謝（10 genes）和嘧啶代謝（1 gene），表明DElncRNA廣泛參與意蜂工蜂中腸發(fā)育過(guò)程中的物質(zhì)代謝調(diào)控。物質(zhì)合成與代謝往往伴隨著旺盛的能量代謝，本研究發(fā)現(xiàn)共有15個(gè)上下游基因富集在甲烷代謝（5 genes）、氧化磷酸化（5 genes）和硫代謝（2 genes）等能量代謝通路，表明相應(yīng)的DElncRNA在意蜂工蜂中腸的能量代謝方面具有重要的調(diào)控功能。蜜蜂腸道是物質(zhì)消化吸收的主要位置。本研究發(fā)現(xiàn)28個(gè)上下游基因富集在7個(gè)消化系統(tǒng)通路上，包括胰腺分泌（8 genes）、蛋白消化和吸收（2 genes）、胃酸分泌（7 genes）、脂肪體消化和吸收（4 genes）、唾液分泌（4 genes）、維生素消化和吸收（1 genes）和膽汁分泌（2 genes），表明DElncRNA參與調(diào)控意蜂工蜂中腸對(duì)物質(zhì)的消化和吸收。Hippo信號(hào)通路可以調(diào)節(jié)器官的大小[41-42]，也能與其他信號(hào)通路相互作用共同調(diào)節(jié)中腸組織的穩(wěn)態(tài)[43]。Notch信號(hào)影響細(xì)胞正常形態(tài)發(fā)生的多個(gè)過(guò)程，包括多能細(xì)胞分化、細(xì)胞凋亡、細(xì)胞增殖及細(xì)胞邊界的形成[44]。本研究發(fā)現(xiàn)，分別有36、19和8個(gè)上下游基因分別富集在Hippo、Wnt和Notch信號(hào)通路。此外，還發(fā)現(xiàn)7、3、3和7個(gè)上下游基因富集在TGF-beta、mTOR、Hedgehog和胰島素信號(hào)通路。結(jié)果表明相應(yīng)的DElncRNA可通過(guò)參與調(diào)控上述信號(hào)通路對(duì)意蜂工蜂中腸發(fā)育進(jìn)行調(diào)控。

角質(zhì)層和圍食膜是昆蟲(chóng)免疫防御的第一道防線[45]。如果病原微生物突破第一道防線，包括蜜蜂在內(nèi)的昆蟲(chóng)將會(huì)啟動(dòng)細(xì)胞和體液免疫抵抗病原的入侵。蜜蜂的細(xì)胞免疫包括內(nèi)吞作用、黑化作用、吞噬作用、蛋白的酶促水解等，體液免疫主要為抗菌肽的合成與釋放[46]。本研究中，DElncRNA上下游基因的KEGG pathway富集分析結(jié)果顯示分別有11、9、7、5和4個(gè)上下游基因富集在溶酶體、內(nèi)吞作用、泛素蛋白水解、吞噬體和黑化作用等細(xì)胞免疫通路，此外，分別有13、3、1和1個(gè)上下游基因富集在MAPK、Jak-STAT、NF-Kappa B和Toll-like受體等體液免疫通路。上述結(jié)果表明相應(yīng)的DElncRNA同時(shí)參與了宿主的細(xì)胞和體液免疫的調(diào)控過(guò)程，推測(cè)其在宿主的免疫防御過(guò)程發(fā)揮關(guān)鍵作用。

Salmena等提出內(nèi)源性競(jìng)爭(zhēng)RNA（ceRNA）假說(shuō)[38]，該假說(shuō)認(rèn)為lncRNA和mRNA等RNA可以通過(guò)miRNA應(yīng)答原件（MRE）與miRNA競(jìng)爭(zhēng)性結(jié)合，從而對(duì)基因表達(dá)進(jìn)行調(diào)控。此后，ceRNA假說(shuō)已被越來(lái)越多的研究[47]所證明。為進(jìn)一步揭示DElncRNA的作用，本研究通過(guò)靶向關(guān)系預(yù)測(cè)DElncRNA結(jié)合的miRNA以及miRNA結(jié)合的mRNA，并構(gòu)建三者間的調(diào)控網(wǎng)絡(luò)，發(fā)現(xiàn)部分上調(diào)與下調(diào)lncRNA位于調(diào)控網(wǎng)絡(luò)中心位置且結(jié)合較多的miRNA。其中，部分lncRNA可結(jié)合多個(gè)miRNA，如TCONS_00004891、XR_001704571.1和TCONS_00024939分別結(jié)合36、28和27個(gè)miRNA，表明這些lncRNA可能在意蜂工蜂中腸發(fā)育過(guò)程中發(fā)揮重要的調(diào)控功能。此外，也有部分lncRNA共同靶向結(jié)合一個(gè)miRNA，如有多達(dá)50個(gè)lncRNA能夠靶向結(jié)合ame-miR-6001-3p。COLLINS等在研究真社會(huì)昆蟲(chóng)等級(jí)分化過(guò)程中發(fā)現(xiàn)，bte-miR-6001-3p在歐洲熊蜂（）等級(jí)分化的幼蟲(chóng)中差異表達(dá)顯著，ame-miR-6001-3p在西方蜜蜂等級(jí)分化的幼蟲(chóng)中差異表達(dá)較弱，說(shuō)明miR-6001-3p影響蜜蜂卵巢發(fā)育[48]。上述結(jié)果表明DElncRNA可作為ceRNA吸附結(jié)合miRNA，從而對(duì)意蜂工蜂中腸發(fā)育過(guò)程進(jìn)行調(diào)控。目前，lncRNA的功能研究多集中在人類、哺乳動(dòng)物和重要疾病領(lǐng)域[17]，在果蠅中也有少量研究[49]，但包括蜜蜂[18-21]在內(nèi)的絕大多數(shù)昆蟲(chóng)的lncRNA研究仍處于初級(jí)階段，相關(guān)的lncRNA信息十分有限。本研究預(yù)測(cè)出的多數(shù)lncRNA，其具體的調(diào)控關(guān)系尚不明確，需要進(jìn)一步克隆出lncRNA全長(zhǎng)序列并深入探究其功能，這是未來(lái)的研究方向。

本研究?jī)H對(duì)意蜂7和10日齡工蜂中腸進(jìn)行了測(cè)序和lncRNA相關(guān)分析，若要全面解析意蜂工蜂中腸的發(fā)育機(jī)理及調(diào)控機(jī)制，需要對(duì)更多日齡的工蜂中腸進(jìn)行測(cè)序，進(jìn)而在全局水平進(jìn)行更加深入的分析，例如對(duì)所有DElncRNA進(jìn)行基因權(quán)重共表達(dá)分析（WGCNA）。

4 結(jié)論

通過(guò)對(duì)意蜂工蜂中腸發(fā)育過(guò)程中的lncRNA及其功能進(jìn)行深入分析，提供了意蜂工蜂中腸發(fā)育過(guò)程中的DElncRNA信息，揭示了DElncRNA廣泛參與意蜂工蜂中腸的新陳代謝、細(xì)胞活動(dòng)和免疫調(diào)控并作為競(jìng)爭(zhēng)性內(nèi)源RNA（ceRNA）發(fā)揮作用，為關(guān)鍵lncRNA的篩選和功能研究提供了必要的數(shù)據(jù)支持。

[1] PARK D, JUNG J W, CHIO B S, JAYAKODI M, LEE J, LIM J, YU Y, CHOI Y S, LEE M L, PARK Y, CHOI I Y, YANG T Y, EDWARDS O R, NAH G, KWON H W. Uncovering the novel characteristics of Asian honey bee,by whole genome sequencing., 2015, 16(1): 1.

[2] National Research Council, Division on Earth and Life Studies, Board on Agriculture and Natural Resources, Board on Life Sciences. Committee on the Status of Pollinators in North America.. Washington, D.C: National Academies Press, 2007.

[3] 周冰峰. 蜜蜂飼養(yǎng)管理學(xué). 廈門: 廈門大學(xué)出版社, 2002.

ZHOU B F.. Xiamen: Xiamen University Publishing Company, 2002. (in Chinese)

[4] DJEBALI S, DAVIS C A, MERKEL A, DOBIN A, LASSMANN T, MORTAZAVI A, TANZER A, LAGARDE J, LIN W, SCHLESINGER F, XUE C, MARINOV G K, KHATUN J, WILLIAMS B A, ZALESKI C, ROZOWSKY J, RODER M, KOKOCINSKI F, ABDELHAMID R F, ALIOTO T, ANTOSHECHKIN I, BAER M T, BAR N S, BATUT P, BELL K, BELL I, CHAKRABORTTY S, CHEN X, CHRAST J, CURADO J, DERRIEN T, DRENKOW J, DUMAIS E, DUMAIS J, DUTTAGUPTA R, FALCONNET E, FASTUCAM, FEJES-TOTH K, FERREIRA P, FOISSAC S, FULLWOOD M J, GAO H, GONZALEZ D, GORDON A, GUNAWARDENA H, HOWALD C, JHA S, JOHNSON R, KAPRANOV P, KING B, KINGSWOOD C, LUO O J, PARK E, PERSAUD K, PREALL J B, RIBECA P, RISK B, ROBYR D, SAMMETH M, SCHAFFER L, SEE L H, SHAHAB A, SKANCKE J, SUZUKI A M, TAKAHASHI H, TILGNER H, TROUT D, WALTERS N, WANG H, WROBEL J, YU Y, RUAN X, HAYASHIZAKI Y, HARROW J, GERSTEIN M, HUBBARD T, REYMOND A, ANTONARAKIS S E, HANNON G, GIDDINGS M C, RUAN Y, WOLD B, CARNINCI P, GUIGO R, GINGERAS T R.Landscape of transcription in human cells., 2012, 489(7414): 101-108.

[5] GORODKIN J, HOFACKER I L. From structure prediction to genomic screens for novel non-coding RNAs., 2011, 7(8): e1002100.

[6] KANDURI C. Kcnq1ot1: a chromatin regulatory RNA., 2011, 22(4): 343-350.

[7] TRIPATHI V, SHEN Z, CHAKRABORTY A, GIRI S, FREIER S M, WU X, ZHANG Y, GOROSPE M, PRASANTH S G, LAL A, PRASANTH K V. Long noncoding RNA MALAT1 controls cell cycle progression by regulating the expression of oncogenic transcription factor B-MYB., 2013, 9(3): e1003368.

[8] LI M, SUN X, CAI H, SUN Y, PLATH M, LI C, LAN X, LEI C, LIN F, BAI Y, CHEN H. Long non-coding RNA ADNCR suppresses adipogenic differentiation by targeting miR-204., 2016, 1859(7): 871-882.

[9] ULITSKY I, BARTEL D P. LincRNAs: genomics, evolution, and mechanisms., 2013, 154(1): 26-46.

[10] HUNG T, WANG Y, LIN M F, KOEGEL A K, KOTAKE Y, GRANT G D, HORLINGS H M, SHAH N, UMBRICHT C, WANG P, WANG Y, KONG B, LANEROD A, BORRESEN-DALE A L, KIM S K, VAN D V M, SUKUMAR S, WHITFIELD M L, KELLIS M, XIONG Y, WONG D J, CHANG H Y. Extensive and coordinated transcription of noncoding RNAs within cell cycle promoters., 2011, 43(7): 621-629.

[11] YOON J H, ABDELMOHSEN K, SRIKANTAN S, YANG X, MARTINDALE J L, DE S, HUARTE M, BECKER K G, GOROSPE M. LincRNA-p21 suppresses target mRNA translation., 2012, 47(4): 648-655.

[12] CHOONIEDASS-KOTHARI S, EMBERLEY E, HAMEDANI M K, TROUP S, WANG X, CZOSNEK A, HUBE F, MUTAWE M, WATSON P H, LEYGUE E. The steroid receptor RNA activator is the first functional RNA encoding a protein., 2004, 566(1/3): 43-47.

[13] OGAWA Y, SUN B K, LEE J T. Intersection of the RNA interference and X-inactivation pathways., 2008, 320(5881): 1336-1341.

[14] KENIRY A, OXLEY D, MONNIER P, KYBA M, DANDOLO L, SMITS G, REIK W. The H19 lincRNA is a developmental reservoir of miR-675 which suppresses growth and Igf1r., 2012, 14(7): 659-665.

[15] MUDGE J M, HARROW J. Creating reference gene annotation for the mouse C57BL6/J genome assembly., 2015, 26(9/10): 366-378.

[16] HON C C, RAMILOWSKI J A, HARSHBARGER J, BERTIN N, GOUGH J, DENISENKO E, SCHMEIER S, POULSEN T M, SEVERIN J, LIZIO M, KAWAJI H, KASUKAWA T, ITOH M, BURROUGHS A M, NOMA S, DJEBALI S, ALAM T, MEDVEDEVA Y A, TESTA A C, LIPOVICH L, YIP C W, ABUGESSAISA I, MENDEZ M, HASEGAWA A, TANG D, LASSMANN T, HEUTINK P, BABINA M, WELLS C A, KOJIMA S, NAKAMURA Y, SUZUKI H, DAUB C O, DE-HOON M J, ARNER E, HAYASHIZAKI Y, CARNINCI P, FORREST A R. An atlas of human long non-coding RNAs with accurate 5′ ends., 2017, 543(7644): 199-204.

[17] 朱斌, 梁沛, 高希武. 長(zhǎng)鏈非編碼RNA (lncRNA)及其在昆蟲(chóng)學(xué)研究中的進(jìn)展. 昆蟲(chóng)學(xué)報(bào), 2016, 59(11): 1272-1281.

ZHU B, LIANG P, GAO X W. Long noncoding RNAs (lncRNAs) and their research advances in entomology., 2016, 59(11): 1272-1281. (in Chinese)

[18] 郭昱, 蘇松坤, 陳盛祿, 張少吾, 陳潤(rùn)生. LncRNA在蜜蜂級(jí)型分化中的功能研究. 生物化學(xué)與生物物理進(jìn)展, 2015, 42(8): 750-757.

GUO Y, SU S K, CHEN C L, ZHANG S W, CHEN R S. The function of lncRNAs in the caste determination of the honeybee., 2015, 42(8): 750-757. (in Chinese)

[19] HUMANN F C, TIBERIO G J, HARTFELDER K. Sequence and expression characteristics of long noncoding RNAs in honey bee caste development–potential novel regulators for transgressive ovary size., 2013, 8(10): e78915.

[20] CHEN X, MA C, CHEN C, LU Q, SHI W, LIU Z G, WANG H H, GUO H K. Integration of lncRNA-miRNA-mRNA reveals novel insights into oviposition regulation in honey bees., 2017, 5: e3881.

[21] JAYAKODI M, JUNG J W, PARK D, AHN Y J, LEE S C, SHIN S Y, CHIN C,YANG T J, KWON H W. Genome-wide characterization of long intergenic non-coding RNAs (lincRNAs) provides new insight into viral diseases in honey beesand., 2015, 16(1): 680.

[22] BABENDREIER D, JOLLER D, ROMEIS J, BIGLER F, WIDMER F. Bacterial community structures in honeybee intestines and their response to two insecticidal proteins., 2007, 59(3): 600-610.

[23] KOCH H, ABROL D P, LI J, SCHMID-HEMPEL P. Diversity and evolutionary patterns of bacterial gut associates of corbiculate bees., 2013, 22(7): 2028-2044.

[24] ELLEGAARD K M, TAMARIT D, JAVELIND E, OLOFSSON T C, ANDERSSON S G, VASQUEZ A. Extensive intra-phylotype diversity in lactobacilli and bifidobacteria from the honeybee gut., 2015, 16(1): 284.

[25] MOHR K I, TEBBE C C. Diversity and phylotype consistency of bacteria in the guts of three bee species (Apoidea) at an oilseed rape field., 2006, 8(2): 258-272.

[26] ENGEL P, MARTINSON V G, MORAN N A. Functional diversity within the simple gut microbiota of the honey bee., 2012, 109(27): 11002-11007.

[27] GILLIAM M, VALENTINE D K. Enterobacteriaceae isolated from foraging worker honey bees,., 1974, 23(1): 38-41.

[28] LANGMEND B, TRAPNELL C, POP M, SALZBERG S L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome., 2009, 10(3): R25.

[29] KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R, SALZBERG S L. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions., 2013, 14(4): R36.

[30] Honeybee Genome Sequencing Consortium.Insights into social insects from the genome of the honeybee., 2006, 443(7114): 931-949.

[31] KONG L, ZHANG Y, YE Z Q, LIU X Q, ZHAO S Q, WEI L, GAO G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine., 2007, 35(Web Server issue): 345-349.

[32] SUN L, LUO H, BU D, ZHAO G, YU K, ZHANG C, LIU Y, CHEN R, ZHAO Y.Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts., 2013, 41(17): e166.

[33] ROBINSON M D, MCCARTHY D J, SMYTH G K. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

[34] REHMSMEIER M, STEFFEN P, HOCHSMANN M, GIEGERICH R. Fast and effective prediction of microRNA/target duplexes., 2004, 10(10): 1507-1517.

[35] BETEL D, WILSON M, GABOW A, MARKS D S, SANDER C. The microRNA.org resource: targets and expression., 2008, 36(Database issue): 149-153.

[36] ALLEN E, XIE Z, GUSTAFSON A M, CARRINGTON J C. MicroRNA-directed phasing during trans-acting siRNA biogenesis in plants., 2005, 121(2): 207-221.

[37] SMOOT M E, ONO K, RUSCHEINSKI J, WANG P L, IDEKER T. Cytoscape 2.8: new features for data integration and network visualization., 2011, 27(3): 431-432.

[38] SALMENA L, POLISENO L, TAY Y KATS L, PANDOLFI P P. Ahypothesis: The Rosetta stone of a hidden RNA language?, 2011, 146(3): 353-358.

[39] FLIRIAN K, JOSHUA T M. Functional classi?cation and experimental dissection of long noncoding RNAs., 2018, 172(3): 393-407.

[40] 趙亞周, 田文禮, 胡熠凡, 彭文君. 蜜蜂蜂王漿主蛋白 (MRJPs)的研究進(jìn)展. 應(yīng)用昆蟲(chóng)學(xué)報(bào), 2012, 49(5): 1345-1353.

ZHAO Y Z, TIAN W L, HU Y F, PENG W J. Research advances in major royal jelly proteins of honeybee., 2012, 49(5): 1345-1353. (in Chinese)

[41] HALDER G, JOHNSON R L. Hippo signaling: growth control and beyond., 2011, 138(1): 9-22.

[42] PAN D. The hippo signaling pathway in development and cancer., 2010, 19(4): 491-505.

[43] CAMARGO F D, GOKHALE S, JOHNNIDIS J B, FU D, BELL G W, JAENISCH R, BRUMMELKAMP T R. YAP1 increases organ size and expands undifferentiated progenitor cells., 2007, 17(23): 2054-2060.

[44] BITEAU B, HOCHMUTH C E, JASPER H.JNK activity in somatic stem cells causes loss of tissue homeostasis in the aginggut., 2008, 3(4): 442-455.

[45] ORIHEL T C. The peritrophic membrane: its role as a barrier to infection of the arthropod host//. Academic Press, 1975: 65-73.

[46] ARONSTEIN K A, MURRAY K D. Chalkbrood disease in honey bees., 2010, 103(Suppl. 1): 20-29.

[47] KARRENTH F A, TAY Y, PERNA D, ALA U, TAN S M, RUST A G, DENICOLA G, WEBSTER K A, WEISS D, PEREZ-MANCERA P A, KRAUTHAMMER M, HALABAN R, PROVERO P, ADAMS D J, TUVESON D A, PANDOLFI P P.identification of tumor suppressive PTEN ceRNAs in an oncogenic BRAF-induced mouse model of melanoma., 2011, 147(2): 382-395.

[48] COLLINS D H, MOHORIANU I, BECKERS M, MOULTON V, DALMAY T, BOURKE A F. MicroRNAs associated with caste determination and differentiation in a primitively eusocial insect., 2017, 7: 45674.

[49] Li E H, Zhao X, Zhang C, LIU W. Fragile X mental retardation protein participates in non-coding RNA pathways., 2018, 40(2): 87-94.

（責(zé)任編輯 岳梅）

Differential expression analysis of long non-coding RNAs during the developmental process ofworker’s Midgut

GUO Rui1, GENG Sihai1, XIONG Cuiling1, ZHENG Yanzhen1, FU Zhongmin1, WANG Haipeng1, DU Yu1, TONG Xinyu1, ZHAO Hongxia2, CHEN Dafu1

(1College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou 350002;2Guangdong Institute of Applied Biological Resources, Guangzhou 510260)

【Objective】Long non-coding RNA (lncRNA) plays an important role in regulation of gene expression, epigenetics and cell cycle in eukaryotes. The objective of this study is to investigate the expression profile and role of lncRNAs in the developmental process ofworker’s midgut. 【Method】In this study, 7- and 10-day-old worker’s midguts of(Am7, Am10) were sequenced using RNA-seq technology and strand-specific library construction method. using Perl script, raw reads were filtered to obtain clean reads with high-quality. Bowtie tool was used to compare clean reads to the ribosome database, and TopHat2 software was employed to compare unmapped clean reads to the reference genome. CPC and CNCI softwares were utilized to predict coding capacity of the transcripts. RT-PCR was performed to identify partial lncRNAs. Investigation of differentially expressed lncRNAs (DElncRNAs) was carried out with edgeR, followed by prediction of upstream and downstream genes, for which GO and KEGG pathway enrichment analyses were performed. RNAhybrid, Miranda and TargetScan softwares were utilized together to predict target miRNAs of DElncRNAs and target genes of miRNAs, and DElncRNAs-miRNAs-mRNAs regulation networks were visualized via Cytoscape. Finally, RT-qPCR was conducted to verify reliability of the sequencing data.【Result】134 802 058 and 147 051 470 raw reads were gained from deep sequencing of Am7 and Am10, respectively, and after stringent filtration, 134 166 157 and 146 293 288 were obtained. In total, 6 353 lncRNAs were predicted, and 3 890 DElncRNAs were obtained based on expression calculation, including 2 005 up-regulated lncRNAs and 1 885 down-regulated lncRNAs. The result of RT-PCR suggested the expected signal bands could be amplified from 8 lncRNAs, implying their true existence. There were 1 793 upstream and downstream genes of DElncRNAs, which were involved in 42 GO terms, including metabolic processes, developmental processes, cellular processes, stress responses, immune system processes and so forth. They were also associated with 251 KEGG pathways, including material metabolism pathways such as carbon metabolism, purine metabolism and fatty acid biosynthesis; energy metabolism pathways such as sulfur metabolism, methane metabolism and oxidative phosphorylation; signaling pathways such as Hippo, Wnt and Notch signaling pathways; cellular immune pathways such as lysosome, endocytosis and ubiquitin mediated proteolysis; humoral immune pathways such as MAPK, Jak-STAT and NF-kappa B pathways, these results demonstrated the DElncRNAs were involved in the material and energy metabolism, cell life activity and immunity regulation in the developmental process ofworker’s midgut. Further analysis showed TCONS_00020918 might play a regulatory part in the nutrient absorption and caste differentiation in the worker’s midgut.Analysis of regulation networks demonstrated that complex networks existed between DElncRNAs and target miRNAs and mRNAs, partial DElncRNAs lie in the central of the networks and link many miRNAs, and partial miRNAs could be bound by many DElncRNAs, which indicated that these DElncRNAs might play an important role during the developmental process of the worker’s midgut. Finally, 5 DElncRNAs were randomly selected for RT-qPCR assay, and the result proved the reliability of sequencing data in this study.【Conclusion】DElncRNA is widely involved in the metabolism, cellular activity and immune regulation ofworker’s midgut, and plays a role as a competitive endogenous RNA (ceRNA). The results provide the necessary data support for the screening and functional study of key lncRNA.

; midgut; development; long non-coding RNA; upstream and downstream genes

10.3864/j.issn.0578-1752.2018.18.016

2018-03-17；

2018-05-08

國(guó)家自然科學(xué)基金（31702190）、國(guó)家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)專項(xiàng)資金（CARS-44-KXJ7）、福建省教育廳中青年教師教育科研項(xiàng)目（JAT170158）、福建農(nóng)林大學(xué)科技創(chuàng)新專項(xiàng)基金（CXZX2017343）、福建省大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃（201610389053）

郭睿，E-mail：ruiguo@fafu.edu.cn。耿四海，E-mail：15737313592@163.com。郭睿和耿四海為同等貢獻(xiàn)作者。通信作者陳大福，E-mail：dfchen826@fafu.edu.cn