意大利蜜蜂幼蟲腸道發育過程中的差異表達 microRNA及其調控網絡

郭睿，杜宇，熊翠玲，鄭燕珍，付中民，徐國鈞，王海朋，陳華枝，耿四海，周丁丁，石彩云，趙紅霞，陳大福

?

意大利蜜蜂幼蟲腸道發育過程中的差異表達 microRNA及其調控網絡

郭睿1，杜宇1，熊翠玲1，鄭燕珍1，付中民1，徐國鈞1，王海朋1，陳華枝1，耿四海1，周丁丁1，石彩云1，趙紅霞2，陳大福1

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

【目的】微小RNA（microRNA，miRNA）是一類在轉錄后水平對mRNA進行負調控的關鍵調控因子。本研究旨在通過分析意大利蜜蜂（，簡稱意蜂）幼蟲腸道發育過程中差異表達miRNA（differentially expressed miRNA，DEmiRNA）及其調控網絡，提供miRNA的表達譜和差異表達信息，揭示DEmiRNA在幼蟲腸道發育中的作用。【方法】利用small RNA-seq（sRNA-seq）技術對意蜂4、5和6日齡幼蟲腸道樣品（Am4、Am5和Am6）進行測序，將質控后的數據與西方蜜蜂（）參考基因組進行比對，然后將比對上的序列標簽（tags）注釋到miRBase數據庫，并利用TPM（tags per million）算法歸一化處理所得miRNA的表達量，再通過相關生物信息學軟件對miRNA進行表達量聚類、前體二級結構預測及差異表達分析。利用TargetFinder軟件預測DEmiRNA的靶基因，使用Blast軟件對靶基因進行GO和KEGG數據庫注釋，進而通過Cytoscape軟件構建miRNA-mRNA的調控網絡。采用莖環實時熒光定量PCR（Stem-loop RT-qPCR）驗證測序數據的可靠性?！窘Y果】意蜂幼蟲腸道樣品的測序分別得到10 841 644、12 037 678和9 230 496條有效序列標簽；Am4 vs Am5比較組包含16個上調和10個下調miRNA，Am5 vs Am6比較組包含5個上調和7個下調miRNA。其中，novel-m0031-3p為兩個比較組所共有，并結合5個與蛻皮激素誘導蛋白相關的靶基因;二者的特有DEmiRNA數分別為25和11個。Am4 vs Am5的26個DEmiRNA結合5 742個靶基因，其中2 725個靶基因可注釋到GO數據庫中的46個GO term，并主要富集在結合、細胞進程、代謝進程和單組織進程等；Am5 vs Am6的12個DEmiRNA結合3 733個靶基因，且其中2 725個靶基因富集在結合、細胞進程、單組織進程和代謝進程等41個GO term；此外，兩個比較組中的DEmiRNA分別有1 046和676個靶基因可注釋到116和92條KEGG代謝通路，且Am4 vs Am5比較組的DEmiRNA的靶基因富集在Wnt信號通路、Hippo信號通路、嘌呤代謝和內吞作用等通路上的數量均高于Am5 vs Am6比較組。進一步分析結果顯示Am4 vs Am5中的上調和下調miRNA可分別結合611和85個靶基因，其中ame-miR-6052結合的靶基因數最多，可通過結合5個靶基因,參與對細胞色素P450的調控；miR-281-x結合49個靶基因,并間接調控組氨酸代謝、TGF-信號通路以及Hippo信號通路；Am5 vs Am6中的上調和下調miRNA可分別結合43和431個靶基因，其中miR-iab-4-x結合靶基因數量最多，并廣泛參與調控背腹軸的形成、Hippo信號通路、Wnt信號通路、FoxO信號通路、Notch信號通路以及mTOR信號通路等與生長發育相關的通路。調控網絡分析結果表明，DEmiRNA與靶基因間形成較為復雜的調控網絡，DEmiRNA居于調控網絡的中心位置，而mRNA位于調控網絡的外周。最后，通過莖環實時熒光定量PCR對隨機選取的3個DEmiRNA進行驗證，結果證實了測序數據的可靠性?！窘Y論】在全基因組水平對意蜂幼蟲腸道的DEmiRNA及其靶基因進行預測和分析，并對DEmiRNA的調控網絡進行構建及分析，發現意蜂幼蟲通過調節包括ame-miR-6052、miR-iab-4-x、miR-281-x和novel-m0031-3p在內的多個miRNA的表達水平對腸道生長和發育進行調控，研究結果不僅提供了意蜂腸道發育過程的miRNA表達譜及差異表達信息，也為闡明意蜂幼蟲腸道的發育機理打下基礎。

意大利蜜蜂；幼蟲腸道；發育；微小RNA；靶基因；調控網絡

0 引言

【研究意義】微小RNA（microRNA，miRNA）是一類長度約為18—25個核苷酸的非編碼RNA（non-coding RNA，ncRNA），通過對mRNA的負調控作用而廣泛參與動物[1-2]、植物[3-4]和微生物[5-6]的各類生物學過程。意大利蜜蜂（，簡稱意蜂）是我國養蜂生產中的主要蜂種，在經濟創收和生態維持等方面具有重要價值。利用small RNA-seq（sRNA-seq）技術對意蜂4、5和6日齡幼蟲腸道進行測序，通過生物信息學分析方法對腸道發育過程中的差異表達miRNA（differentially expressed miRNA，DEmiRNA）進行全面分析，并分別構建上調和下調miRNA的調控網絡，對揭示意蜂幼蟲腸道發育過程中miRNA的差異表達規律、腸道發育相關關鍵miRNA的篩選和功能研究以及闡明腸道發育的機理具有重要意義。【前人研究進展】西方蜜蜂（）基因組[7]的公布為其組學和分子生物學研究奠定了基礎。此后，隨著高通量測序技術、生物信息學分析方法和軟件的革新、融合及應用，蜜蜂的轉錄組學研究取得了較大進展[8-9]。對蜜蜂的miRNA開展了一些研究，如Ashby等[10]研究發現，miR-184、miR-bantam和miR315參與了蜜蜂的細胞分化、組織結構重塑以及級型分化；陳曉[11]對蜂王卵巢激活、產卵抑制和產卵恢復過程中的DEmiRNA進行靶基因預測和分析，發現多個DEmiRNA與卵巢激活和產卵密切相關；石元元[12]研究發現，ame-bantam、ame-let-7和ame-miR-8參與了東方蜜蜂（）4日齡蜂王幼蟲和工蜂幼蟲的Wnt信號通路，并與雌性蜜蜂級型分化的分子機理密切相關。此外，蜜蜂腸道不僅是食物消化、營養吸收和利用的主要場所，也是抵御病原入侵的重要免疫器官。前人對蜜蜂腸道的研究主要集中在腸道菌群的群落結構和功能預測方面[13-14]，但有關蜜蜂及其幼蟲腸道發育的研究極為滯后，腸道發育的分子機理仍不明確。筆者所在課題組前期已在mRNA組學水平對意蜂幼蟲腸道進行了全面的轉錄組學研究，通過差異表達基因（DEG）和趨勢分析揭示了意蜂幼蟲腸道發育過程中的基因表達譜和差異表達規律[15]，意蜂幼蟲腸道響應球囊菌（）脅迫的免疫應答[16]，以及意蜂幼蟲與球囊菌之間的互作[17]。【本研究切入點】目前為止，意蜂幼蟲腸道發育過程中的miRNA表達譜仍然缺失，相關miRNA在幼蟲腸道發育過程中的作用還不明確?！緮M解決的關鍵問題】結合sRNA-seq技術對意蜂4、5和6日齡幼蟲腸道樣品進行測序和分析，通過DEmiRNA及其靶基因的預測和分析、調控網絡的構建和分析，提供意蜂幼蟲腸道發育過程中的miRNA表達譜和差異表達信息，揭示DEmiRNA在腸道發育中的作用。

1 材料與方法

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

1.1 生物材料

供試意蜂幼蟲取自福建農林大學蜂學學院教學蜂場。

1.2 意蜂幼蟲人工飼養及樣品測序

按照王倩等[18]的方法配制幼蟲飼料，從健康意蜂蜂群中提取巢脾，將2日齡幼蟲移至已預置50 μL飼料的24孔細胞培養板中，在35℃，相對濕度（RH）為90%的培養箱中飼養，每隔24 h更換飼料。分別剖取4、5和6日齡意蜂幼蟲腸道（分別記為Am4、Am5和Am6），裝入EP管后放入液氮速凍，隨后保存在-80℃的超低溫冰箱。進行3次生物學重復。Am4的3個生物學重復分別為Am4-1、Am4-2和Am4-3，Am5的3個生物學重復分別為Am5-1、Am5-2和Am5-3，Am6的3個生物學重復分別為Am6-1、Am6-2和Am6-3。上述9個腸道樣品委托廣州基迪奧生物科技有限公司對進行單端測序，測序平臺為Illumina MiSeq。原始數據已上傳NCBI SRA數據庫，BioProject號：PRJNA408312。

1.3 測序數據的質控及定位

對原始下機數據，利用Perl腳本剔除銜接子（adaptor）、未知堿基N和低質量的reads，從而得到高質量的序列（clean reads），保證后續數據分析的準確性。通過Bowite軟件將各樣品的tags序列比對到GenBank及Rfam（11.0）數據庫，過濾核糖體RNA（rRNA）、胞質小RNA（scRNA）以及核內小RNA（snoRNA）等ncRNA和重復序列，得到miRNA的非注釋標簽序列（unannotated tags）。繼而利用Bowit 軟件將非注釋標簽序列與西方蜜蜂的參考基因組（assembly Amel 4.5）的序列進行比對，得到相關tags在參考基因組上的位置信息，即為mapped tags。

1.4 DEmiRNA的預測

利用miRDeep2軟件[19]將mapped tags與miRBase數據庫中已知miRNA前體序列進行比對，從而鑒定已知miRNA的表達，同時得到可能的前體序列。對各樣本中miRNA進行表達量的統計，并通過TPM（tags per million）算法公式（TPM=T×106/N，T表示miRNA的tags，N表示總miRNA的tags）對全部miRNA的表達量進行歸一化處理。利用R軟件計算各樣品之間的相關性系數。DEmiRNA（Am4 vs Am5、Am5 vs Am6和Am4 vs Am6）的篩選標準為-value≤0.05且|log2Fold change|≥1。

1.5 DEmiRNA靶基因預測及分析

根據DEmiRNA與對應物種的基因序列信息，利用TargetFinder軟件[20]進行靶基因預測，并利用Blast軟件將預測靶基因序列與GO（Gene Ontology）、KEGG數據庫比對，獲得靶基因的注釋信息。利用RNAhybrid、TargetScan和Miranda軟件預測DEmiRNA結合的靶基因，根據上述靶向結合關系構建miRNA-mRNA的調控網絡并通過Cytoscape軟件[21]將其可視化。

1.6 DEmiRNA的莖環實時熒光定量PCR（Stem-loop RT-qPCR）驗證

為了驗證sRNA-seq數據的可靠性，隨機選取3個DEmiRNA（miR-7964-y、miR-8516-x和miR-3747-x）進行Stem-loop RT-qPCR驗證。根據所選DEmiRNA的序列，參照Chen等[22]的方法，利用DNAMAN軟件（Lynnon Biosoft公司，美國）設計特異性的Stem-loop引物、上游引物和下游引物，委托上海生工生物工程有限公司進行引物合成。選擇snRNA U6作為內參。利用RNA抽提試劑盒（Axygen公司，美國）分別提取Am4、Am5和Am6的總RNA，利用Stem-loop引物進行反轉錄得到相應的cDNA，作為模板進行qPCR。反應體系（20 μL）中含有SYBR Green Dye 10 μL，上下游引物各1 μL，cDNA模板DNA 1 μL，Rox 0.44 μL，DEPC水補至20 μL。在ABI 7500熒光定量PCR儀（ABI公司，美國）中進行反應，反應條件：95℃預變性1 min，95℃變性15 s，48℃延伸30 s，共40個循環，最后72℃延伸45 s。所選miRNA的相對表達量采用2-△△Ct法計算。每個反應進行3個生物學重復和3次平行重復。

2 結果

2.1 數據質控與評估

本研究中，3個幼蟲腸道樣品的sRNA-seq分別產生13 186 921、14 255 967和10 921 897條clean reads，經嚴格過濾后得到的clean tags數分別為10 841 644（82.22%）、12 037 678（84.44%）和9 230 496（84.51%）條（表1）。Am4、Am5和Am6組內各生物學重復之間的Pearson相關系數均在0.9734以上，說明各樣品的重復性較好（圖1）。上述結果說明本研究的測序數據質量良好，可用于進一步分析。

表1 sRNA-seq數據統計

2.2 意蜂幼蟲腸道的DEmiRNA分析

Am4 vs Am5比較組中有26個miRNA差異表達，包括16個上調miRNA和10個下調miRNA。Am5 vs Am6 比較組共有12個DEmiRNA，包括5個上調miRNA和7個下調miRNA。Am 4 vs Am6比較組中有41個miRNA差異表達，包括15個上調miRNA和26個下調miRNA。分別對Am4 vs Am5和Am5 vs Am6的DEmiRNA進行表達量聚類分析，結果顯示2個比較組中DEmiRNA的變化倍數差異明顯，多數DEmiRNA的差異表達幅度較?。▓D2-A、2-B）。進一步分析發現，Am4 vs Am5和Am5 vs Am6中DEmiRNA分別包含3個和2個novel miRNA，其中novel-m0031-3p為2個比較組所共有。利用miRDeep2軟件對novel miRNA的前體二級結構進行預測，結果顯示它們均有標志性的莖環結構（圖2-C—2-F）。

2.3 意蜂幼蟲腸道的DEmiRNA的靶基因預測及功能注釋

利用TargetFinder軟件對意蜂幼蟲腸道的DEmiRNA進行靶基因預測，Am4 vs Am5、Am5 vs Am6的DEmiRNA可分別預測出5 742和3 733個靶基因。對上述靶基因進行GO數據庫注釋，結果顯示Am4 vs Am5的DEmiRNA的2 725個靶基因共涉及46個GO term，富集基因數最多的是結合（1 612 gene）、細胞進程（1 537 gene）、代謝進程（1 232 gene）、單組織進程（1 197 gene）、催化活性（1 020 gene）、細胞膜（761 gene）、細胞（757 gene）、細胞組件（757 gene）、細胞膜組件（745 gene）、細胞器（578 gene）等（圖3）；Am5 vs Am6的DEmiRNA的1 785個靶基因可注釋到41個GO term，富集基因數最多的是結合（1 068 gene）、細胞進程（1 026 gene）、單組織進程（847 gene）、代謝進程（799 gene）、催化活性（591 gene）、細胞膜（533 gene）、細胞膜組件（529 gene）、細胞（448 gene）、細胞組件（448 gene）、生物學調控（348 gene）等（圖3）。說明隨著發育時間的延長，意蜂幼蟲腸道的DEmiRNA數、靶基因數及其涉及的GO term數逐漸減少；DEmiRNA廣泛參與意蜂幼蟲腸道的新陳代謝、細胞生命活動及免疫防御。

圖1 各意蜂幼蟲腸道樣品的不同生物學重復間的Pearson相關性

進一步對DEmiRNA的靶基因進行KEGG代謝通路（pathway）富集分析，結果顯示Am4 vs Am5中的1 046個DEmiRNA的靶基因可注釋到116條pathway，其中富集基因數最多的是Wnt信號通路（140 gene）、嘌呤代謝（87 gene）、Hippo信號通路（80 gene）、內吞作用（64 gene）、光傳導（60 gene）、神經活性配體-受體相互作用（57 gene）、FoxO信號通路（43 gene）、內質網中蛋白質的加工（41 gene）、磷脂酰肌醇信號系統（38 gene）、泛素介導的蛋白水解（36 gene）等（圖4-A），說明相應的miRNA參與到意蜂4和5日齡幼蟲腸道發育過程中的新陳代謝、蛋白質合成、免疫防御以及相關信號通路的調控。

Am5 vs Am6中的676個DEmiRNA的靶基注釋到92條pathway，并主要富集在Wnt信號通路（109 gene）、光傳導（69 gene）、Hippo信號通路（66 gene）、神經活性配體-受體相互作用（57 gene）、嘌呤代謝（42 gene）、內吞作用（38 gene）、背腹軸形成（30 gene）、mRNA監視（28 gene）、Hedgehog信號通路、晝夜節律（24 gene）等（圖4-B），說明相應的miRNA同樣廣泛參與了意蜂幼蟲5和6日齡腸道的生長發育過程中各類新陳代謝以及信號通路的調控過程。

2.4 意蜂幼蟲腸道DEmiRNA的調控網絡分析

利用軟件預測DEmiRNA結合的靶基因并通過Cytoscape軟件進行可視化，結果顯示Am4 vs Am5中上調miRNA可結合611個靶基因，下調miRNA可結合85個靶基因，上調（或下調）miRNA與靶基因形成較為復雜的調控網絡，DEmiRNA居于調控網絡的中心位置，而靶基因處于調控網絡的外周；其中，所有DEmiRNA均可連接2個靶基因以上，ame-miR- 6052結合的靶基因數多達204個（圖5-A、5-B）。Am5 vs Am6中上調和下調miRNA可分別結合43和431個靶基因，上調（或下調）miRNA同樣與靶基因形成復雜的調控網絡，所有DEmiRNA均可連接3個以上的靶基因，其中miR-iab-4-x結合的靶基因數最多，達到125個（圖5-C、5-D）。

A：Am4 vs Am5中DEmiRNA的表達量聚類Expression clustering of DEmiRNA in Am4 vs Am5；B：Am5 vs Am6中DEmiRNA的表達量聚類Expression clustering of DEmiRNA in Am5 vs Am6；C：novel-m0004-3p前體的二級結構Secondary structure of precursor of novel-m0004-3p；D：novel-m0017-3p前體的二級結構Secondary structure of precursor of novel-m0017-3p；E：novel-m0031-3p前體的二級結構Secondary structure of precursor of novel-m0031-3p；F：novel-m0037-5p前體的二級結構Secondary structure of precursor of novel-m0037-5p。黃色區域為novel miRNA的成熟序列Yellow regions indicate mature sequences of novel miRNA

圖3 意蜂幼蟲腸道DEmiRNA的靶基因的GO數據庫注釋

2.5 意蜂幼蟲腸道DEmiRNA的RT-qPCR驗證

隨機挑取3個DEmiRNA（miR-7964-y、miR-8516-x和miR-3747-x）進行RT-qPCR驗證，結果顯示它們的表達水平的變化趨勢和測序數據中相應DEmiRNA的表達水平的變化趨勢一致（圖6），說明本研究中的測序數據真實可靠。

A：Am4 vs Am5中的DEmiRNA的靶基因Target genes of DEmiRNA in Am4 vs Am5；B：Am5 vs Am6中的DEmiRNA的靶基因Target genes of DEmiRNA in Am5 vs Am6

A：Am4 vs Am5中上調miRNA的miRNA-mRNA網絡miRNA-mRNA networks of up-regulated miRNA in Am4 vs Am5；B：Am4 vs Am5中下調miRNA的miRNA-mRNA網絡miRNA-mRNA networks of down-regulated miRNA in Am4 vs Am5；C：Am5 vs Am6中上調miRNA的miRNA-mRNA網絡miRNA-mRNA networks of up-regulated miRNA in Am5 vs Am6；D：Am5 vs Am6中下調miRNA的miRNA-mRNA網絡miRNA-mRNA networks of down-regulated miRNA in Am5 vs Am6

A：miR-7964-y；B：miR-8516-x；C：miR-3747-x

3 討論

MiRNA作為一種重要的基因表達調控因子，在昆蟲的各種生物學進程發揮關鍵作用[23]。家蠶（）和小菜蛾（）不同發育階段和不同組織中的miRNA表達譜研究揭示了miRNA在其生長發育[24]、變態[25]及行為反應[26]等過程中的重要調控功能。相比于果蠅、家蠶等模式昆蟲，蜜蜂的miRNA研究相對滯后。陳璇等[27-28]曾對三型蜂不同發育階段的混合小RNA文庫進行了測序，發現了267個novel miRNA，并進一步對蜂王和工蜂4和5日齡幼蟲的DEmiRNA及其靶基因進行了分析，為蜜蜂發育和級型分化關鍵時期的miRNA調控網絡研究打下了基礎。目前，有關蜜蜂幼蟲腸道發育過程miRNA表達譜的研究十分滯后，miRNA在腸道發育中的調控機理仍不明確。本研究利用sRNA-seq技術對意蜂4、5和6日齡幼蟲腸道進行測序，通過生物信息學方法對鄰近日齡幼蟲腸道的miRNA進行差異表達分析，分別預測出26和12個DEmiRNA，說明DEmiRNA的數量在意蜂幼蟲腸道發育過程中有逐漸減少的趨勢；進一步分析發現Am4 vs Am5和Am5 vs Am6中的特有DEmiRNA數分別為25（miR-342-y、ame-miR-3759和miR-4331-y等）和11個（miR-4955-x、miR-46-y和ame-miR-3478等），二者包含1個共有DEmiRNA（novel-m0031-3p），且novel-m0031-3p的表達水平隨日齡的增加呈逐漸下調趨勢，推測其通過下調表達量減少對靶基因的抑制作用，從而在意蜂幼蟲腸道的不同發育階段均發揮基礎性的調控作用，而特有DEmiRNA在幼蟲腸道的不同發育階段發揮特殊的調控功能。

蜜蜂幼蟲腸道內存在簡單的共生菌，為盡量減小意蜂幼蟲腸道內共生菌對測序數據的影響，本研究一方面在人工飼喂過程中盡量減少環境微生物的帶入，例如對所用器具進行高壓蒸汽滅菌，用75%酒精對人工氣候箱進行擦洗消毒，更換飼料時在酒精燈旁操作等；另一方面通過對原始數據進行嚴格的過濾和質控、比對西方蜜蜂基因組以消除腸道共生菌數據的影響，昆蟲和微生物的物種親緣關系較遠，二者的基因保守性很低，因而比對西方蜜蜂基因組的數據理論上應為意蜂幼蟲腸道本身的數據。

蜜蜂腸道不僅是食物消化、營養吸收和利用的主要場所，也是抵御病原入侵的重要免疫器官。本研究發現，Am5 vs Am6中的12個DEmiRNA的靶基因富集在了各類能量和物質代謝通路，如氮代謝（1 gene）、氨基酸代謝（32 gene）、碳水化合物代謝（69 gene）和脂質代謝（74 gene）等，表明相應的DEmiRNA參與了意蜂幼蟲腸道能量和物質代謝的調控；還發現部分靶基因注釋到與生長、發育相關的代謝通路和信號通路，包括背腹軸形成（30 gene）、Hippo信號通路（66 gene）、Wnt信號通路（109 gene）等，表明相應的DEmiRNA在生長和發育過程中發揮重要的調控作用。此外，分別有38和13個靶基因注釋到內吞作用和泛素介導的蛋白質水解等細胞免疫通路，分別有6和4個靶基因注釋到MAPK和Jak-STAT等體液免疫通路，表明相應的DEmiRNA參與了意蜂幼蟲腸道的免疫防御過程。本研究中，Am4 vs Am5中的26個DEmiRNA的靶基因同樣也廣泛參與生長發育、新陳代謝以及免疫防御相關的代謝通路。Am4 vs Am5中的DEmiRNA及其靶基因數均高于Am5 vs Am6，意蜂5日齡幼蟲體積近乎4日齡幼蟲的2倍，但6日齡幼蟲體積稍大于5日齡幼蟲，推測意蜂4—5日齡幼蟲的發育時期需要更多的DEmiRNA參與到生長發育、新陳代謝等方面的調控。

蛻皮激素（ecdysteroid，Ec）是昆蟲體內的重要激素之一，參與昆蟲生長、變態和生殖的整個生命活動，其滴度的階段性增加是昆蟲生長和發育的必要條件[29]。novel-m0031-3p在Am4 vs Am5和Am5 vs Am6中均有表達，且其結合的5個靶基因（XM_006564318.2、XM_006564319.2、XM_006564320.2、XM_006564321.2和XM_016914663.1）均涉及Ec誘導蛋白的調控，有研究表明蛻皮觸發激素基因可在家蠶幼蟲蛻皮和變態發育過程中起到關鍵調控作用[30]，因此推測novel-m0031-3p通過調控Ec滴度使其達到動態平衡，從而保證意蜂幼蟲腸道的正常發育。miR-281是一種高度保守的miRNA，在家蠶幼蟲階段中呈高量表達，與神經發育、組織生長密切相關[25]。周艷河[31]研究發現，受到高劑量登革病毒侵染的白紋伊蚊（）中腸內miR-281會靶向調控自身基因從而影響病毒的復制水平；熊慧萍[32]對位于可變內含子區的黑腹果蠅（）的miR-281-1/2基因轉錄和啟動子進行分析，發現miR-281只在3齡幼蟲、蛹和成蟲期表達，且在不同發育階段miR-281轉錄起始位點的轉錄效率存在差異。本研究中，miR-281-x是Am4 vs Am5中的特有DEmiRNA，其結合的49個靶基因涉及TGF-信號通路（8 gene）、Hippo信號通路（7 gene）、背腹軸的形成（3 gene）和組氨酸代謝（3 gene）等代謝通路，推測其參與對意蜂幼蟲腸道早期生長和發育的調控，但miR-281-x的時空表達譜需進一步研究。miR-iab-4可通過調控Hox編碼蛋白的基因表達間接影響黑腹果蠅的翅形成過程，miR-iab-4的突變會影響果蠅幼蟲的自身調節能力[33]。本研究中，miR-iab-4-x為Am5 vs Am6中的特有DEmiRNA，其結合的125個靶基因涉及Hippo信號通路（14 gene）、Wnt信號通路（9 gene）、FoxO信號通路（8 gene）、Notch信號通路（7 gene）、mTOR信號通路（6 gene）和背腹軸的形成（5 gene）等與生長發育相關的代謝通路，推測其參與意蜂幼蟲腸道發育后期的調控，miR-iab-4-x的下調表達可能對意蜂幼蟲腸道發育過程發揮重要作用。

一個miRNA可以同時靶向調控多個mRNA，反之亦然。為進一步揭示DEmiRNA的作用，本研究通過序列匹配關系預測DEmiRNA結合的靶基因，并構建了二者的調控網絡（圖5），發現部分上調和下調的miRNA位于調控網絡的中心位置且結合較多的mRNA，如ame-miR-6052與miR-iab-4-x可分別結合204和125個靶基因，具有很高的連通性，表明二者可能在意蜂幼蟲腸道的生長和發育過程中發揮關鍵的調控功能。細胞色素P450作為一種重要的解毒酶，廣泛存在于昆蟲的脂肪體、馬氏管和中腸內，且在中腸的含量最高；其也可參與內源性物質代謝，在生物體內發揮重要的作用[34-35]。本研究發現Am4 vs Am5中特有的ame-miR-6052表達量上調，其結合的5個靶基因（XM_006559340.2、XM_006559341.1、XM_016912202.1、XM_016916903.1和XM_623618.5）與細胞色素P450的調節有關，推測ame-miR-6052可參與意蜂幼蟲腸道對異源物質的降解、新陳代謝過程的調控，下一步可通過合成miRNA mimic、miRNA inhibitor對其進行過表達和敲減，從而深入探究其功能。

4 結論

結合高通量測序技術和生物信息學分析方法對意蜂幼蟲腸道的DEmiRNA進行了全基因組水平的預測和分析，并對DEmiRNA-mRNA調控網絡進行構建及分析，發現意蜂幼蟲通過調節包括ame-miR-6052、miR-iab-4-x、miR-281-x和novel-m0031-3p在內的多個miRNA的表達水平對腸道生長和發育進行調控，研究結果不僅提供了意蜂腸道發育過程的miRNA表達譜及差異表達信息，也為闡明意蜂幼蟲腸道的發育機理打下了基礎。

[1] Azzam G, Smibert P, Lai E C, LIU J L.Argonaute 1 and its miRNA biogenesis partners are required for oocyte formation and germline cell division., 2012, 365(2): 384-394.

[2] Asgari S. MicroRNA functions in insects., 2013, 43(4): 388-397.

[3] Xu L, Hu Y G, Cao Y, LI J R, MA L G, LI Y, QI Y J. An expression atlas of miRNAs in., 2018, 61(2): 178-189.

[4] Liu M, Yu H, Zhao G, HUANG Q, LU Y, OUYANG B. Profiling of drought-responsive microRNA and mRNA in tomato using high-throughput sequencing., 2017, 18(1): 481.

[5] Yu Z Q, Gao X L, Liu C N, LV X P, ZHENG S M. Analysis of microRNA expression profile in specific pathogen-free chickens in response to reticuloendotheliosis virus infection., 2017, 101(7): 2767-2777.

[6] Ojha C R, Rodriguez M, Dever S M, MUKHOPADHYAY R, EL-HAGE N. Mammalian microRNA: an important modulator of host-pathogen interactions in human viral infections., 2016, 23(1): 74.

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

[8] 劉芳. 意蜂哺育蜂與采集蜂頭部mRNAs與miRNAs表達譜Solexa測序比較分析及其調控網絡研究[D]. 杭州: 浙江大學, 2012.

LIU F. Integrating of Solexa high-abundance mRNAs and miRNAs in: comparison between nurses and foragers to identify regulatory network[D]. Hangzhou: Zhejiang University, 2012. (in Chinese)

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

GUO Y, SU S K, CHEN S 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)

[10] Ashby R, FORêT S, Searle I, MALESZKA R. MicroRNAs in honey bee caste determination., 2016, 6: 18794.

[11] 陳曉. 蜜蜂卵巢激活和產卵過程差異表達的編碼RNA與非編碼RNA的篩選和鑒定[D]. 北京: 中國農業科學院, 2017.

CHEN X. Identification of differentially expressed coding and noncoding RNAs during ovary activation and oviposition in honey bees[D]. Beijing: Chinese Academy of Agricultural Sciences, 2017. (in Chinese)

[12] 石元元. 東方蜜蜂遺傳圖譜構建以及雌性蜜蜂發育分子機理[D]. 南昌: 江西農業大學, 2014.

SHI Y Y. Construction of genetic linkage map inand molecular mechanism of development in females honey bee[D]. Nanchang: Jiangxi Agricultural University, 2014. (in Chinese)

[13] 常偉. 蜜蜂消化道共生細菌及其多態性研究初探[D]. 福州: 福建農林大學, 2010.

CHANG W. Diversity of symbiotic bacteria in honeybee alimentary tract[D]. Fuzhou: Fujian Agriculture and Forestry University, 2010. (in Chinese)

[14] 賈慧茹. 亞致死劑量吡蟲啉對意大利蜜蜂中腸菌群的影響[D]. 北京: 中國農業科學院, 2015.

JIA H R. Effect of the sublethal doses of imidacloprid on the bacterial diversity in the midgut of[D]. Beijing: Chinese Academy of Agricultural Sciences, 2015. (in Chinese)

[15] 郭睿, 解彥玲, 熊翠玲, 尹偉軒, 鄭燕珍, 付中民, 陳大福. 意大利蜜蜂4、5和6日齡幼蟲腸道發育過程中差異表達基因的趨勢分析. 上海交通大學學報(農業科學版), 2018, 36(4): 14-21, 29.

GUO R, XIE Y L, XIONG C L, YI W X, ZHENG Y Z, FU Z M, CHEN D F. Trend analysis for differentially expressed genes in developmental process of 4-, 5- and 6-day-old larval guts of, 2018, 36(4): 14-21, 29. (in Chinese)

[16] Chen D, Guo R, Xu X, XIONG C l, LIANG Q, ZHENG Y z, LUO Q, ZHANG Z, HUANG Z J, KUMAR D, XI W J, ZOU X, LIU M. Uncovering the immune responses of, larval gut to, infection utilizing transcriptome sequencing., 2017, 621: 40-50.

[17] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 黃枳腱, 張曌楠, 張璐, 李汶東, 童新宇, 席偉軍. 脅迫意大利蜜蜂幼蟲腸道的球囊菌的轉錄組分析. 昆蟲學報, 2017, 60(4): 401-411.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, HUANG Z J, ZHANG Z N, ZHANG L, LI W D, TONG X Y, XI W J. Transcriptomic analysis ofstressing larval gut of(Hyemenoptera: Apidae)., 2017, 60(4): 401-411. (in Chinese)

[18] 王倩, 孫亮先, 肖培新, 劉鋒, 康明江, 胥保華. 室內人工培育中華蜜蜂幼蟲技術研究. 山東農業科學, 2009(11): 113-116.

WANG Q, SUN L X, XIAO P X, LIU F, KANG M J, XU B H. Study on technology for indoor artificial feeding oflarvae., 2009(11): 113-116. (in Chinese)

[19] Friedl?nder M R, Mackowiak S D, Li N, CHEN W, RAJEWSKY N. MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades., 2012, 40(1): 37-52.

[20] 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.

[21] 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.

[22] Chen C, Ridzon D A, Broomer A J, ZHOU Z, LEE D H, NGUYEN J T, BARBISIN M, XU N L, MAHUVAKAR V R, ANDERSEN M R, LAO K Q, LIVAK K J, GUEGLER K J. Real-time quantification of microRNAs by stem-loop RT-PCR., 2005, 33(20): e179.

[23] Lucas K J, Zhao B, Liu S, RAIKHEL A S. Regulation of physiological processes by microRNAs in insects., 2015, 11: 1-7.

[24] Zhang Y, Zhou X, Ge X, LI M, JIA S, YANG X, KAN Y, MIAO X, ZHAO G, LI F, HUANG Y. Insect-specific microRNA involved in the development of the silkworm., 2009, 4(3): e4677.

[25] Liu S, Gao S, Zhang D, YIN J, XIANG Z, XIA Q. MicroRNAs show diverse and dynamic expression patterns in multiple tissues of., 2010, 11: 85.

[26] Liang P, Feng B, Zhou X G, GAO X W. Identification and developmental profiling of microRNAs in diamondback moth,(L.)., 2013, 8(11): e78787.

[27] 陳璇, 俞曉敏, 鄭火青, 蔡亦梅, 胡福良. 西方蜜蜂(L.) sRNA的富集與文庫檢測. 中國農業科學, 2009, 42(8): 2943-2948.

CHEN X, YU X M, ZHENG H Q, CAI Y M, HU F L. Separation and enrichment of sRNAs from honeybee (L.) and its quality detection by library construction., 2009, 42(8): 2943-2948. (in Chinese)

[28] 陳璇. 蜜蜂() microRNA的全基因組挖掘及在雌性蜜蜂級型分化關鍵時期轉錄組水平調控作用[D]. 杭州: 浙江大學, 2012.

CHEN X. Genome-wide identification of microRNAs and their regulation of transcriptome on female caste determination of honey bee ()[D]. Hangzhou: Zhejiang University, 2012. (in Chinese)

[29] 汝玉濤, 王勇, 周敬林, 王德意, 馬月月, 姜義仁, 高清, 秦利. 蛻皮激素受體和超氣門蛋白基因在柞蠶發育過程及激素誘導后的表達模式. 蠶業科學, 2017(4): 594-602.

RU Y T, WANG Y, ZHOU J L, WANG D Y, MA Y Y, JIANG Y R, GAO Q, QIN L. The expression patterns of ecdysone receptor and ultraspiracle genes induring development and hormone-induced process., 2017(4): 594-602. (in Chinese)

[30] 舒旭. 家蠶蛻皮觸發激素及其受體基因的克隆和表達分析[D]. 重慶: 西南大學, 2009.

SHU X. Cloning and expression analysis of genes encoding ecdysis triggering hormone and its receptor in the silkworm,[D]. Chongqing: Southwest University, 2009. (in Chinese)

[31] 周艷河. 白紋伊蚊中腸特異性高表達miRNA-miR-281對登革病毒復制的調節作用[D]. 廣州: 南方醫科大學, 2014.

ZHOU Y H. Dengue virus replication is regulated by miR-281: an abundant midgut-specific miRNA of vector mosquito[D]. Guangzhou: Southern Medical University, 2014. (in Chinese)

[32] 熊慧萍. 位于可變內含子區的果蠅microRNA-281-1/2基因轉錄和啟動子分析[D]. 南京: 南京農業大學, 2008.

XIONG H P. Transcription and promoter analysis ofintronic microRNA-281-1/2 located in alternative spliced region[D]. Nanjing: Nanjing Agricultural University, 2008. (in Chinese)

[33] Ronshaugen M, Biemar F, Piel J, LEVINE M, LAI E C. ThemicroRNA iab-4 causes a dominant homeotic transformation of halteres to wings., 2005, 19(24): 2947-2952.

[34] TAMáSI V, Monostory K, Prough R A, Falus A. Role of xenobiotic metabolism in cancer: involvement of transcriptional and miRNA regulation of P450s., 2011, 68(7): 1131-1146.

[35] 黃獻彬. 小菜蛾解毒代謝相關miRNA鑒定及表達分析[D]. 福州: 福建農林大學, 2016.

HUANG X B. Identification and analysis of miRNAs involved in detoxification in the diamondback moth,[D]. Fuzhou: Fujian Agriculture and Forestry University, 2016. (in Chinese)

Differentially expressed microRNA and their regulation networks during the developmental process oflarval gut

Guo Rui1, Du Yu1, Xiong CuiLing1, Zheng YanZhen1, Fu ZhongMin1, Xu GuoJun1, Wang HaiPeng1, Chen HuaZhi1, Geng SiHai1, Zhou DingDing1, Shi CaiYun1, Zhao HongXia2, Chen DaFu1

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

【Objective】MicroRNA (miRNA) is a kind of key regulator for negative regulation of mRNA at post-transcriptional level. The objective of this study is to provide miRNA expression patterns and differential expression information, illuminate the function of differentially expressed miRNA (DEmiRNA) in the development of larval gut by comprehensively investigating the DEmiRNAs and their regulation networks during the developmental process oflarval gut.【Method】Deep sequencing of the 4-, 5- and 6-day-old larval guts ofwas conducted using small RNA-seq (sRNA-seq) technology, followed by mapping of the data after quality-control with the reference genome of,and the mapped tags were then compared to miRBase database. The miRNA expression level was normalized by TPM algorithm, and the expression clustering, prediction of secondary structure of precursor and differential expression analysis were performed using related bioinformatic softwares. TargetFinder software was used to predict target gene of DEmiRNA, which was annotated to GO and KEGG databases using Blast, furthermore, miRNA-mRNA regulation networks were constructed using Cytoscape software. Stem-loop RT-qPCR was used to verify the sequencing data in this study.【Result】High-throughput sequencing of larval gut samples produced 10 841 644, 12 037 678 and 9 230 496 clean tags, respectively. In Am4 vs Am5 comparison group, there were16 up-regulated and 10 down-regulated miRNAs, while Am5 vs Am6 comparison groupincluded 5 up-regulated and 7 down-regulated miRNAs, respectively. Among them, novel-m0031-3p was shared by both Am4 vs Am5 and Am5 vs Am6, binding 5 target genes associated with ecdysone inducible protein, 25 and 11 DEmiRNAs were specific for the above-mentioned two comparison groups.DEmiRNA in Am4 vs Am5 could bind 5 742 target genes, among them 2 725 targets could be annotated to 46 GO terms in GO database, and the largest ones were binding, cellular process, metabolic process and single-organism process. similarly, 12 DEmiRNAs in Am5 vs Am6 could link 3 733 target genes, among them 2 725 targets could be annotated to 41 GO terms, and mostly enriched terms were binding, cellular process, single-organism processand metabolic process. In addition, 1 046 and 676 target genes of two comparison groups were related to 116 and 92 KEGG pathways, and the number of DEmiRNA target genes in Am4 vs Am5 was more than that in Am5 vs Am6, which annotated to Wnt signaling pathway, Hippo signaling pathway, purine metabolism and endocytosis. further analysis demonstrated that up-regulated and down-regulated miRNAs in Am4 vs Am5 could bind 611 and 85 target genes, and ame-miR-6052 linked the most target genes and participated in regulating cytochrome P450 via 5 target genes.miR-281-x could bind 49 target genes and indirectly regulate histidine metabolism, TGF-signaling pathway and Hippo signaling pathway.InAm5 vs Am6 comparison group, up-regulated and down-regulated miRNAs could bind 43 and 431 target genes, respectively, among them miR-iab-4-x linked the most target genes, and it could participate in regulating growth and development related pathways, such as dorso-ventral axis formation, Hippo signaling pathway, Wnt signaling pathway, FoxO signaling pathway, Notch signaling pathway and mTOR signaling pathway. Regulation network analysis indicated that complex networks formed between DEmiRNAs and target genes, and DEmiRNAs lied in the center while target genes lied in the periphery. Finally, Stem-loop RT-qPCR was carried out to validate the randomly selected three DEmiRNAs, and the result confirmed the reliability of sequencing data. 【Conclusion】The DEmiRNA and corresponding target genes in thelarval gut were predicted and analyzed at genome-wide level, it was found thatare capable of regulating the expression of many miRNAs such as ame-miR-6052, miR-iab-4-x, miR-281-x and novel-m0031-3p. The results not only offer the expression pattern and differential expression information of miRNA during the developmental process oflarval gut, but also lay a foundation for clarifying the molecular mechanisms underlying the larval gut’s development.

; larval gut; development; microRNA; target gene; regulation network

10.3864/j.issn.0578-1752.2018.21.018

2018-05-22；

2018-06-28

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

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

（責任編輯 岳梅）