意大利蜜蜂工蜂中腸發育過程中長鏈非編碼RNA的差異表達分析

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

?

意大利蜜蜂工蜂中腸發育過程中長鏈非編碼RNA的差異表達分析

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

（1福建農林大學蜂學學院，福州 350002；2廣東省生物資源應用研究所，廣州 510260）

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

意大利蜜蜂；中腸；發育；長鏈非編碼RNA；上下游基因

0 引言

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

1 材料與方法

試驗于2017年在福建農林大學蜂學學院蜜蜂保護實驗室進行。

1.1 生物材料

供試意大利蜜蜂工蜂取自福建農林大學蜂學學院教學蜂場。

1.2 樣品制備與Illumina測序

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

圖1 意蜂工蜂的人工飼養

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

1.3 高通量數據分析

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

lncRNA與mRNA相比在序列保守型、ORF長度等特點上具有一定差異，利用軟件CPC[31]和CNCI[32]對新轉錄本的編碼能力進行預測，選取二者的交集作為可靠的預測結果。

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

1.3.3 DElncRNA的調控網絡構建 LncRNA可以競爭結合miRNA，從而減少miRNA結合mRNA，表明lncRNA可以通過miRNA調控mRNA的表達。利用RNAhybrid[34]、Miranda[35]和TargetScan[36]軟件預測DElncRNA靶向結合的miRNA，以及miRNA結合的mRNA，根據上述靶向結合關系構建lncRNA-miRNA- mRNA的調控網絡。利用Cytoscape軟件[37]對上述調控網絡進行可視化。

1.4 lncRNA的RT-PCR驗證

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

1.5 DElncRNA的實時熒光定量PCR（RT-qPCR）驗證

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

2 結果

2.1 數據質控與評估

在實驗室條件下對意蜂工蜂進行人工飼養（圖1），Am7和Am10中腸樣品的測序得到raw reads平均分別為134 802 058和147 051 470條，過濾后得到clean reads平均分別為134 166 157和146 293 288條。各樣品的平均Q20和Q30分別為97.34%和93.86%（表1）。Am7與Am10的組內各生物學重復之間的Pearson相關系數均在0.97以上，說明各樣品的重復性較好（圖2）。上述結果說明本研究的測序數據質量良好，可用于進一步分析。

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

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

表1 RNA-seq數據概覽

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的差異表達分析結果顯示，在Am7 vs Am10比較組中共有3 890個DElncRNA，包括2 005個上調lncRNA和1885個下調lncRNA。lncRNA可通過影響其編碼基因的上下游基因而發揮作用。共預測出1 698個DElncRNA的上下游基因，其中上游基因、下游基因和重疊基因分別為485、535和593個。對DElncRNA編碼基因的上下游基因進行GO分類，結果顯示它們分布在生物學進程、細胞組分和分子功能3大類，涉及42個GO term，富集基因數最多的前10個GO term分別為結合（355 genes）、細胞進程（315 genes）、代謝進程（299 genes）、單細胞進程（262 genes）、催化活性（259 genes）、細胞（147 genes）、細胞組分（147 genes）、細胞膜（138 genes）、細胞膜組分（131 genes）和細胞器（117 genes）（圖4）。

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

2.3 DElncRNA的調控網絡分析

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

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

2.4 DElncRNA的qRT-PCR驗證

隨機挑取6個DElncRNA進行RT-qPCR驗證，結果顯示其中有5個DElncRNA表達水平的變化趨勢和轉錄組數據中相應DElncRNA表達水平的變化趨勢一致（圖7），說明本研究中的測序數據真實可靠。相關引物序列信息詳見表2。對驗證的5個DElncRNA進行上下游基因分析，發現TCONS_00029692可調控肥大細胞脫粒肽編碼基因，TCONS_00038012可調控未知蛋白LOC100576733編碼基因。

3 討論

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

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

A：上調lncRNA的lncRNA-miRNA-mRNA網絡lncRNA-miRNA-mRNA networks of up-regulated lncRNAs；B：下調lncRNA的lncRNA-miRNA-mRNA網絡 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的功能與其編碼基因坐標臨近的蛋白編碼基因相關，位于上游的lncRNA可能與啟動子或共表達基因的其他順式作用元件有交集，從而在轉錄或轉錄后水平對基因的表達進行調控[39]。本研究中，DElncRNA的上下游基因達到1 697個，GO富集分析結果顯示分別有299個和4個上下游基因涉及代謝進程和發育進程，有259個上下游基因涉及氧化磷酸化，表明DElncRNA參與意蜂工蜂中腸的物質和能量代謝過程調控。本研究發現，315個上下游基因涉及細胞進程，表明DElncRNA廣泛參與意蜂工蜂中腸的細胞生命活動。此外，還發現分別有92個和3個上下游基因涉及應激反應和免疫系統進程，暗示相應的DElncRNA在意蜂工蜂中腸發育中發揮調控作用。王漿主蛋白（MRJP）為幼蟲發育提供必要的可利用氮元素，在蜜蜂行為中有潛在的調控作用，還參與決定蜂群中雌性蜂的級型分化等，其中MRJP1在蜂王漿主蛋白中含量最豐富，所占比例達31%，占蜂王漿中水溶性蛋白的48%[40]。本研究發現，TCONS_00020918可調控西方蜜蜂MRJP1編碼基因，表明該lncRNA可能在意蜂的營養吸收、級型分化中發揮重要的調控功能。此外，多個DElncRNA可能調控未知蛋白編碼基因，具體調控關系的明確有賴于Nr數據庫蛋白功能注釋信息的完善及進一步的實驗驗證。

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

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

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

本研究僅對意蜂7和10日齡工蜂中腸進行了測序和lncRNA相關分析，若要全面解析意蜂工蜂中腸的發育機理及調控機制，需要對更多日齡的工蜂中腸進行測序，進而在全局水平進行更加深入的分析，例如對所有DElncRNA進行基因權重共表達分析（WGCNA）。

4 結論

通過對意蜂工蜂中腸發育過程中的lncRNA及其功能進行深入分析，提供了意蜂工蜂中腸發育過程中的DElncRNA信息，揭示了DElncRNA廣泛參與意蜂工蜂中腸的新陳代謝、細胞活動和免疫調控并作為競爭性內源RNA（ceRNA）發揮作用，為關鍵lncRNA的篩選和功能研究提供了必要的數據支持。

[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] 周冰峰. 蜜蜂飼養管理學. 廈門: 廈門大學出版社, 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] 朱斌, 梁沛, 高希武. 長鏈非編碼RNA (lncRNA)及其在昆蟲學研究中的進展. 昆蟲學報, 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] 郭昱, 蘇松坤, 陳盛祿, 張少吾, 陳潤生. LncRNA在蜜蜂級型分化中的功能研究. 生物化學與生物物理進展, 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)的研究進展. 應用昆蟲學報, 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.

（責任編輯 岳梅）

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

國家自然科學基金（31702190）、國家現代農業產業技術體系建設專項資金（CARS-44-KXJ7）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農林大學科技創新專項基金（CXZX2017343）、福建省大學生創新創業訓練計劃（201610389053）

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