意大利蜜蜂工蜂中腸發育過程中的差異基因表達譜及調控網絡

杜宇，周丁丁，萬潔琦，盧家軒，范小雪，范元嬋，陳恒，熊翠玲，鄭燕珍，付中民，徐國鈞，陳大福，郭睿

意大利蜜蜂工蜂中腸發育過程中的差異基因表達譜及調控網絡

杜宇，周丁丁，萬潔琦，盧家軒，范小雪，范元嬋，陳恒，熊翠玲，鄭燕珍，付中民，徐國鈞，陳大福，郭睿

（福建農林大學動物科學學院（蜂學學院），福州 350002）

【】前期已對意大利蜜蜂（，簡稱意蜂）7日齡工蜂中腸（Am7）、10日齡工蜂中腸（Am10）進行全轉錄組測序，本研究基于高質量的組學數據探究中腸發育過程中的差異基因表達譜及其調控網絡，以期解析意蜂工蜂中腸發育的分子機理。根據FPKM（fragments per kilobase of transcript per million mapped reads）算法計算基因表達量，并以|log2fold change|≥1且≤0.05作為標準篩選得到差異表達基因（differentially expressed gene，DEG）。利用TargetFinder軟件預測ame-miR-6001-3p的靶mRNA。利用相關生物信息學軟件，對全部DEG進行GO和KEGG數據庫注釋。篩選出與AMPK、P13K-Akt、Wnt、cAMP、FoxO、Hippo、mTOR、Jak-STAT、Toll-like受體、TGF-beta、Notch、MAPK和NF-κB 13條信號通路存在富集關系的DEG，以及與ame-miR-6001-3p存在靶向結合關系的DEG，通過Cytoscape軟件構建調控網絡將其富集關系與調控關系可視化。利用莖環反轉錄PCR（Stem loop RT-PCR）和實時熒光定量PCR（RT-qPCR）驗證ame-miR-6001-3p以及DEG在Am7和Am10中的差異表達情況。Am7 vs Am10比較組中共有1 038個DEG，包括515個上調基因和523個下調基因。這些DEG涉及細胞進程、代謝進程和催化活性等功能條目，并顯著富集在氧化磷酸化、氨基糖與核苷酸糖代謝、脂肪酸代謝和嘌呤代謝等能量和物質代謝通路，表明工蜂中腸內存在旺盛細胞生命與新陳代謝活動。表達量聚類分析發現分別有20、18、15和14個DEG富集在AMPK信號通路、P13K-Akt信號通路、內吞作用和Hippo信號通路。57個DEG與P13K-Akt、Wnt、Jak-STAT等上述13條與生長發育和免疫防御相關的信號通路存在富集關系，且1個DEG可與多條信號通路存在富集關系。調控網絡分析結果顯示，分別有54個上調基因和44個下調基因可被ame-miR-6001-3p靶向結合；上調基因富集在磷酸肌醇代謝、胰島素信號通路、Hippo信號通路和谷胱甘肽代謝等43條代謝通路，而下調基因富集在Hippo信號通路、新陳代謝途徑、谷胱甘肽代謝和花生四烯酸代謝等20條代謝通路。RT-qPCR結果顯示隨機挑選的6個DEG的表達量變化趨勢與測序數據一致，證實了本研究中基因差異表達真實可靠。此外ame-miR-6001-3p在Am7和Am10內均真實表達，并在Am10中的相對表達量顯著較低。對意蜂工蜂中腸發育過程的DEG表達譜和DEG與ame-miR-6001-3p之間的調控關系，以及DEG的潛在作用進行深入分析和探討，發現DEG可參與TGF-beta、Wnt、Hippo、Notch、PI3K-Akt、mTOR、AMPK和NF-κB等各類信號通路進而影響中腸的生長發育和免疫防御，DEG可通過與顯著下調表達的ame-miR-6001-3p形成復雜的調控網絡參與意蜂工蜂中腸發育過程中胰島素信號通路等多條代謝途徑。

意大利蜜蜂；中腸；發育；差異表達基因；競爭性內源RNA；調控網絡

0 引言

【研究意義】意大利蜜蜂（，簡稱意蜂）是經濟價值和生態價值優越的社會性模式昆蟲，廣泛用于國內外的養蜂生產。蜜蜂中腸作為營養消化吸收、能量物質代謝以及免疫應答防御的主要場所，其發育機理至今仍不明了。因此，通過分析差異表達基因（differentially expressed gene，DEG）在腸道發育過程中的表達譜、DEG的潛在作用及調控網絡，可為明確意蜂中腸發育的分子機理提供理論依據。【前人研究進展】昆蟲腸道的內環境穩態可影響個體的發育進程。COON等研究發現，埃及伊蚊（）和岡比亞按蚊（）腸道內的部分菌群是維持個體正常發育的關鍵因子[1]。抗生素處理可導致斯氏按蚊（）幼蟲發育緩慢[2]，而植物乳桿菌（）可以通過影響激素信號受體促進黑腹果蠅（）的個體增長[3]。目前有關蜜蜂腸道的研究主要涉及腸道菌群的結構和功能預測、外源營養素與腸道發育的關系等方面。ZHENG等[4]通過比較無菌西方蜜蜂（）與正常西方蜜蜂的生長速率和腸道質量，發現腸道菌群會促進西方蜜蜂腸道的發育并增加個體質量；馮倩倩[5]利用添加維生素A的代用花粉飼喂意蜂，通過檢測中腸蛋白質濃度、堿性磷酸酶活性和圍食膜厚度等指標發現適量的維生素A有益于意蜂腸道的發育。近年來以RNA-seq為代表的第二代轉錄組測序技術迅猛發展，為蜜蜂中腸發育機理的深入研究提供了有力工具。筆者團隊前期利用二代測序技術和趨勢分析方法解析了意蜂幼蟲腸道發育過程中的差異基因表達譜[6]，并探究了微小RNA（microRNA，miRNA）在幼蟲腸道發育過程中的調控作用[7]。前人研究發現ame-miR-6001-3p作為關鍵的基因表達調控因子，可參與調控蜜蜂蛻皮激素的分泌[8]和級型分化過程[9]。結合鏈特異性建庫的長鏈非編碼RNA（long non-coding RNA，lncRNA）測序技術、small RNA-seq（sRNA-seq）技術，筆者團隊前期已對意蜂7和10日齡工蜂中腸（Am7和Am10）進行了深度測序，獲得了高質量的全轉錄組數據，并通過深入的生物信息學分析發現意蜂工蜂中腸內的lncRNA和環狀RNA（circular RNA，circRNA）可作為競爭性內源RNA（ceRNA）吸附結合ame-miR-6001-3p，從而間接調控mRNA的表達，影響中腸的細胞分化、新陳代謝、信號傳導和免疫調控[10-11]。【本研究切入點】然而，前期研究均是在單一ncRNA組學層面探究意蜂工蜂中腸的發育機理，缺乏mRNA組學層面的分析和闡釋，也沒有綜合多個層面組學數據對編碼RNA與ncRNA的相互調控關系進行深入探討。【擬解決的關鍵問題】基于已獲得的高質量全轉錄組數據，通過生物信息學方法對意蜂工蜂中腸發育過程中的DEG及其潛在作用進行深入分析，并結合前期在lncRNA組學和circRNA組學層面的分析篩選出的ame-miR-6001-3p，構建和分析mRNA與miRNA之間的調控網絡，以期更深入地解析意蜂工蜂中腸的發育機理。

1 材料與方法

試驗于2017—2019年在福建農林大學動物科學學院（蜂學學院）蜜蜂保護實驗室進行。

1.1 供試蜜蜂

選取福建農林大學教學蜂場內的意蜂蜂群3群，提取9張老熟封蓋子脾至實驗室的培養箱內，箱內溫度（34±0.5）℃，相對濕度70%，收集剛羽化出房工蜂100只。

1.2 意蜂工蜂中腸全轉錄組數據來源

中腸樣品的制備步驟簡述如下：每10只剛羽化（記為0日齡）的工蜂放入1個四周打孔通風的干凈塑料盒內，盒子上方裝有50%（w/v）無菌糖水的飼喂器，置于培養箱中培養（條件同1.1）。在超凈工作臺內分別剖取飼養至7日齡工蜂中腸（Am7）和10日齡工蜂中腸（Am10），每3只中腸置于1個RNA-Free的EP管中，液氮速凍后迅速放入-80℃超低溫冰箱保存。Am7和Am10均設置3個生物學重復。鑒于東方蜜蜂微孢子蟲（）完成一個生活史需要6 d[12]以及10 d仍處于孢子不斷增殖時期[13]，同時為滿足探究意蜂工蜂中腸發育機理及意蜂對主要真菌病原的脅迫應答，本研究選擇上述兩個取樣時間點。

上述6個中腸樣品的cDNA文庫構建、lncRNA- seq實驗技術流程參數參照郭睿等[10]方法，委托廣州基迪奧生物技術有限公司進行，測序平臺分別為Illumina HiseqTM4000。相關原始數據已上傳至NCBI SRA數據庫，BioProject號：PRJNA406998。

1.3 測序數據的質控

利用Perl腳本剔除測序數據原始讀段中未知核苷酸含量＞5%以及含量＞50%的質量低的讀段（reads），獲得有效讀段（clean reads）。過濾可比對上核糖體RNA數據庫的reads，進而通過TopHat2將unmapped reads比對到西方蜜蜂參考基因組（Amel_4.5）（https://www.ncbi.nlm.nih.gov/assembly/GCF_000002195.4）。校準參數如下：（1）最大讀取不匹配為2；（2）配對讀取目標區域為50 bp；（3）配對讀取目標區域誤差為±80 bp。利用Cufflinks將對比到轉錄本上的reads進行組裝。根據FPKM（fragments per kilobase of transcript per million mapped reads）算法計算基因表達量，并利用R軟件計算Pearson相關性系數來評估同組樣品的生物學重復性。

1.4 DEG分析

將|log2fold change|≥1且值≤0.05的基因定義為DEG。利用Blast軟件將DEG注釋到GO數據庫和KEGG數據庫，繪制生物學進程、細胞組分和分子功能三大類的GO富集分類柱狀統計圖，并通過在線分析軟件OmicsShare（http://www.omicshare.com/ tools/Home/Index/index.html）繪制表達量聚類熱圖。ame-miR-6001-3p是基于意蜂工蜂中腸的全轉錄組數據及前期分析結果篩選得到的候選靶標，該靶標可靶向結合多個意蜂lncRNA和circRNA，且位于調控網絡的中心[10-11]。利用TargetFinder軟件對ame-miR- 6001-3p的靶mRNA進行預測，并與本研究的DEG進行比較。最后，根據DEG與ame-miR-6001-3p以及KEGG代謝通路（pathway）的關系構建調控網絡，并利用Cytoscape軟件進行網絡可視化[14]。

1.5 DEG及ame-miR-6001-3p的驗證

為驗證測序數據的可靠性，隨機選取3個上調基因（XM_016912284.1、TCONS_00023605和XM_ 016912389.1）和3個下調基因（TCONS_00001615、XM_001123053.4和XM_016915022.1）進行RT-qPCR驗證。首先利用RNA抽提試劑盒（TaKaRa公司，中國）提取Am7和Am10樣品的總RNA，反轉錄得到cDNA作為模板進行RT-qPCR。按照SYBR Green Dye試劑盒（Vazyme公司，中國）操作說明書進行。反應體系20 μL包含：SYBR Green Dye 10 μL，正、反向引物（10.0 μmol·L-1）各1 μL，cDNA模板1 μL，DEPC水7 μL。RT-qPCR反應在ABI 7500熒光定量PCR儀（ABI公司，美國）上進行，反應條件如下：95℃預變性1 min；95℃變性15 s，60℃延伸30 s，共40個循環；最后72℃延伸45 s。本試驗進行3次技術性重復。以作為內參基因，所選基因的相對表達量采用2-??Ct法計算。

利用總RNA和ame-miR-6001-3p的特異性Stem-loop引物進行反轉錄得到cDNA模板，結合上下游引物，進行常規PCR反應，PCR產物經1.5%瓊脂糖凝膠電泳檢測。以snRNA U6為內參，參照筆者團隊前期已建立的RT-qPCR流程[7,15]對ame-miR- 6001-3p在Am7和Am10中的相對表達量進行檢測。

通過GraphPad Prism軟件繪制上述DEG及ame- miR-6001-3p的相對變化量的柱狀圖，利用SPSS軟件計算組間顯著性。本研究所用相關引物信息詳見表1。

2 結果

2.1 測序數據質控及DEG分析

前期研究中已對Am7和Am10的lncRNA-seq數據進行了嚴格質控，結果表明意蜂工蜂中腸樣品的生物學重復性且全轉錄組數據質量良好[10]，可滿足本研究的分析需要。差異分析結果顯示Am7 vs Am10比較組包含1 038個DEG，其中上調和下調基因的數量分別為515和523個（圖1）。

表1 本研究所用引物

圖1 意蜂工蜂中腸發育過程中DEG的火山圖

2.2 DEG的GO分類

GO分類結果顯示，DEG共注釋到46個功能條目，其中上調和下調基因可分別注釋到36和44個功能條目。進一步分析發現，細胞組分包含基因數比較多的分別是細胞組件（158個DEG）、細胞（158個DEG）、細胞膜（118個DEG）、細胞膜組件（114個DEG）和細胞器（111個DEG）；分子功能包含基因數比較多的分別是結合（350個DEG）、催化活性（258個DEG）、轉運子活性（59個DEG）、分子傳感器活性（29個DEG）和分子功能調節器（28個DEG）；生物學進程包含基因數比較多的分別是細胞進程（285個DEG）、代謝進程（273個DEG）、單組織進程（240個DEG）、生物學調控（112個DEG）和定位（95個DEG）（圖2）。

圖2 意蜂工蜂中腸發育過程中DEG的GO分類

2.3 DEG的KEGG代謝通路富集分析

KEGG代謝通路富集分析結果顯示，DEG共富集在260條代謝通路。其中，AMPK信號通路（20個DEG）、PI3K-Akt信號通路（18個DEG）、內吞作用（15個DEG）、Hippo信號通路（14個DEG）、泛素介導的蛋白質水解（14個DEG）、胰島素信號通路（12個DEG）、磷脂酰肌醇信號轉導系統（12個DEG）、Wnt信號通路（12個DEG）、磷脂酶D信號通路（11個DEG）、鞘脂信號通路（11個DEG）等通路富集的DEG數量較多；此外，部分DEG富集于生長發育、新陳代謝、免疫防御等相關通路，包括溶酶體（11個DEG）、FoxO信號通路（9個DEG）、吞噬體（9個DEG）、cAMP信號通路（9個DEG）、磷酸肌醇代謝（8個DEG）、Hedgehog信號通路（8個DEG）、碳代謝（8個DEG）、mTOR信號通路（7個DEG）、氨基糖與核苷糖代謝（6個DEG）、脂肪酸代謝（6個DEG）、MAPK信號通路（6個DEG）、嘌呤代謝（5個DEG）、Toll-like受體信號通路（4個DEG）。

進一步對富集在AMPK信號通路、PI3K-Akt信號通路、Hippo信號通路和內吞作用上的DEG進行表達量聚類分析，結果顯示各通路富集的上調基因數分別為12、8、13和9個，均高于下調基因數，表明上述4條代謝通路被不同程度的激活（圖3）。

57個DEG在13條信號通路（AMPK、P13K-Akt、Wnt、cAMP、FoxO、Hippo、mTOR、Jak-STAT、Toll-like受體、TGF-beta、Notch、MAPK和NF-B）的富集網絡關系，利用Cytoscape軟件進行可視化，結果顯示一個基因可同時富集在多條信號通路（圖4），例如磷脂酰肌醇3-激酶調節亞基-類編碼基因（XM_016917733.1和XM_392119.6）富集在7條信號通路（AMPK、PI3-Akt、cAMP、FoxO、mTOR、Jak-STAT和Toll-like受體）；Ras相關蛋白rac1亞型編碼基因（TCONS_00021690和TCONS_00038549）富集在5條信號通路（P13K-Akt、cAMP（NF-B、AMPK）、Toll-like受體、Wnt、MAPK）（圖4）；Smad（mothers against decapentaplegic homolog）3亞型編碼基因（TCONS_00001946）與4條信號通路（Wnt、FoxO、Hippo和TGF-beta）存在富集關系；類胰島素樣受體編碼基因（XM_016918001.1）、mTOR的調節相關蛋白編碼基因（TCONS_00038733）、磷酸烯醇丙酮酸羧激酶編碼基因（XM_006560464.2）等8個DEG均與3條信號通路存在富集關系。

利用TargetFinder軟件對ame-miR-6001-3p與本研究中的DEG進行靶向關系預測，發現ame-miR-6001-3p可靶向結合98個DEG，其中上調基因54個，下調基因44個（圖5）。進一步分析發現，ame-miR-6001-3p靶向結合的上調基因富集在新陳代謝途徑（8個上調基因）、磷酸肌醇代謝（2個上調基因）、胰島素信號通路（2個上調基因）、Hippo信號通路（2個上調基因）、谷胱甘肽代謝（1個上調基因）、氧化磷酸化（1個上調基因）、鈣離子信號通路（1個上調基因）、cAMP信號通路（1個上調基因）、FoxO信號通路（1個上調基因）、AMPK信號通路（1個上調基因）和PI3K-Akt信號通路（1個上調基因）等43條代謝通路；靶向結合的下調基因富集在Hippo信號通路（3個下調基因）、新陳代謝途徑（3個下調基因）、谷胱甘肽代謝（1個下調基因）和花生四烯酸代謝（1個下調基因）等20條代謝通路。

2.4 DEG及ame-miR-6001-3p的驗證

利用RT-qPCR對隨機選擇的3個上調基因和3個下調基因進行表達量檢測，結果顯示上述6個DEG與轉錄組測序中相應的基因表達量變化趨勢一致（圖6-A—F），證明了本研究中基因差異表達情況真實可靠。

此外，瓊脂糖凝膠電泳結果顯示ame-miR-6001-3p在Am7和Am10內均有表達（圖6-G）；RT-qPCR檢測結果顯示ame-miR-6001-3p在Am10中的相對表達量較其在Am7中的相對表達量顯著下調（= 0.0003）（圖6-H）。

3 討論

腸道是蜜蜂吸收消化以及免疫防御的主要場所[16-17]，腸道的發育情況及內環境穩態共同影響著不同蟲態蜜蜂的正常生命活動。前期研究中，筆者團隊對意蜂幼蟲腸道發育及響應蜜蜂球囊菌（）脅迫的mRNA和miRNA應答進行了全面細致的組學研究[7,15,18-19]；也分別在lncRNA和circRNA組學層面探究了意蜂工蜂中腸的發育機理[10-11]。基于前期已獲得的意蜂工蜂中腸的全轉錄組數據，本研究對意蜂工蜂中腸發育過程的DEG及其潛在作用，以及融合DEG與ame-miR-6001-3p的調控網絡進行深入分析，發現Am7 vs Am10比較組的上調基因和下調基因數分別為515和523個，其中分別有285、273和258個DEG涉及細胞進程、代謝進程和催化活性等功能條目，分別有13、8、8、6、6和5個DEG富集在氧化磷酸化、磷酸肌醇代謝、碳代謝、氨基糖與核苷酸糖代謝、脂肪酸代謝和嘌呤代謝等能量和物質代謝通路，說明意蜂工蜂中腸發育過程伴隨著較為旺盛細胞生命與新陳代謝活動。

3.1 DEG通過參與復雜的信號通路cross-talk影響意蜂工蜂中腸的發育和免疫

昆蟲體內具有高度復雜的信號通路，同一基因可能涉及多條信號通路（圖4）。有研究表明，TGF-beta信號通路主要參與調節細胞增殖、組織分化和免疫防御等生物學過程，并能與Wnt、FoxO、Hippo、MAPK、Notch等信號通路共同作用，影響昆蟲色素沉淀和附肢發育[20]，蜜蜂幼蟲的腸道發育[7]，以及調控蜂王的卵巢激活和產卵活動[21]等。本研究中，分別有14、12、9、6、3、2個DEG富集在Hippo、Wnt、FoxO、MAPK、TGF-beta、Notch信號通路，其中SMAD3亞型編碼基因（TCONS_00001946）顯著下調表達（=0.0489，log2fold change=-7.1）并同時富集在TGF-beta、Wnt、Hippo和FoxO 4條信號通路。寧越等[22]研究發現，SMAD1可正向調控秦川牛（）原代成肌細胞的分化。SMAD轉導因子位于TGF-beta通路的下游，在果蠅幼蟲發育過程中Activin/TGF-beta通路可通過SMAD2促進細胞的生長[23]。因此，推測SMAD3亞型編碼基因通過TGF-beta信號通路影響細胞的增殖和分化，進而調控意蜂工蜂中腸的發育過程。AMPK信號通路作為一種高度保守的細胞能量調節器，在維持細胞內能量平衡方面發揮重要的調節作用，并可通過介導NF-B信號通路進而抑制具有調節昆蟲蛻皮激素分泌、海藻糖代謝及組織發育等功能的胰島素信號通路的活躍程度[24-25]。本研究中，絲裂原活化蛋白激酶亞型編碼基因（TCONS_00038549）顯著下調表達（=0.0002，log2fold change=-7.3）且同時富集在AMPK和NF-B信號通路；類胰島素受體樣編碼基因（XM_016918001.1和XM_016918002.1）顯著上調（=0.0004，log2fold change=9.2；=0.0031，log2fold change=4.2），并富集在胰島素耐受性信號通路；此外，還發現有3個與海藻糖轉運蛋白相關的編碼基因（XM_016916206.1、NM_001113739.1和XM_006565793.2）的表達水平顯著上調（=0.0292，log2fold change=10.0；=3.46E-39，log2fold change= 9.5；=0.0233，log2fold change=1.5）。推測絲裂原活化蛋白激酶亞型編碼基因通過下調其表達水平，減輕對胰島素信號通路的抑制作用間接增強類胰島素受體樣編碼基因的表達，從而影響促進意蜂工蜂中腸對海藻糖等能量物質的吸收和代謝作用，進而影響中腸的發育過程。

A：AMPK信號通路AMPK signaling pathway；B：內吞作用Endocytosis；C：P13K-Akt信號通路PI3K-Akt signaling pathway；D：Hippo信號通路Hippo signaling pathway

六邊形代表信號通路hexagons indicate signaling pathways；紅色圓形代表上調基因red circles indicate up-regulated genes；綠色圓形代表下調基因green circles indicate down-regulated genes；六邊形和圓形的大小代表連接DEG或信號通路的數量the size of hexagons and circles indicates the number of DEGs or signaling pathways connected

圖5 意蜂工蜂中腸的DEG與ame-miR-6001-3p的調控網絡

圖6 DEG與ame-miR-6001-3p的驗證

PI3K-Akt信號通路作為細胞內重要信號轉導通路之一，通過影響下游多種效應分子的活化狀態，調節昆蟲細胞增殖、組織分化以及個體發育等多個生物學進程[26]。此外，PI3K-Akt和mTOR信號通路可共同作用，通過加強干擾素刺激基因的翻譯過程激活抗病原免疫應激反應[27]。本研究發現，磷脂酰肌醇3-激酶調節亞基-類編碼基因（XM_016917733.1和XM_392119.6）、mTOR的調節相關蛋白編碼基因（TCONS_00038733）和結核蛋白亞型編碼基因（XM_016917423.1）4個DEG同時富集在PI3K-Akt和mTOR信號通路，此外，干擾素相關發育調節因子編碼基因（TCONS_00019888）呈顯著上調表達（=0.0033，log2fold change=12.4）。推測在意蜂工蜂腸道發育過程中PI3K-Akt和mTOR通路共同作用刺激干擾素相關基因的表達，從而激活中腸的抗病原免疫應激反應，這仍有待于進一步的實驗驗證。

昆蟲腸道發育過程復雜，伴隨著非編碼RNA和編碼RNA的相互作用[28]。為探究意蜂幼蟲腸道的發育機理，筆者團隊利用二代測序和趨勢分析初步解析了意蜂幼蟲腸道發育過程中差異基因表達譜，發現呈顯著上調趨勢的DEG廣泛富集在核苷酸代謝、脂代謝、氨基酸代謝等各類新陳代謝通路，以及FoxO、Hippo、mTOR、MAPK、Notch、Wnt、TGF-beta、Hedgehog、Jak-STAT等信號通路上，表明隨著幼蟲腸道發育時間的延長，與發育和免疫相關通路上富集基因的表達水平逐漸增強[6]。此外，筆者團隊還發現在意蜂幼蟲腸道發育過程中，共有560個miRNA與16 479個靶mRNA存在靶向結合關系，宿主miRNA通過調控相關通路及富集基因參與對腸道發育和免疫等生命活動的調控[29]。本研究發現隨著意蜂工蜂中腸發育時間的延長，富集在AMPK、PI3K-Akt、Hippo等免疫和發育相關信號通路上的上調基因數均高于下調基因數，表明上述通路存在不同程度的激活。綜上，推測意蜂幼蟲腸道和工蜂中腸隨著發育歷期的延長，其內部與發育和免疫相關的信號通路及富集基因激活表達，以適應腸道生長和發育所需的物質能量代謝，并維持腸道內環境穩態以應答外界復雜環境。

3.2 意蜂工蜂中腸發育過程中ceRNA調控網絡

LncRNA和circRNA可作為ceRNA，通過miRNA應答元件競爭性結合miRNA，以減輕miRNA對mRNA的負調控作用，間接影響基因的表達，最終調節各類生物進程[30]。郭睿等[10-11]曾對意蜂7和10日齡工蜂中腸發育過程的lncRNA和circRNA的差異表達譜及ceRNA調控網絡進行解析，研究結果顯示這兩種ncRNA均能通過靶向結合ame-miR-6001-3p、let-7-z和miR-136-x等對靶mRNA的表達水平進行間接調控，經比較分析發現，本研究中的283個DEG能被112個DElncRNA靶向的111個miRNA所調控，120個DEG可被141個DEcircRNA靶向的107個miRNA所調控。功能注釋結果顯示上述DEG參與意蜂工蜂中腸發育過程中的各類信號轉導及細胞生命活動。前期分析結果顯示lncRNA和circRNA也可分別通過調節上下游基因和來源基因的表達水平，參與意蜂工蜂中腸發育過程中的TGF-beta、Wnt、Hippo和Notch等信號通路，以及各類細胞生命和新陳代謝活動。上述結果表明意蜂工蜂中腸發育過程伴隨著極其復雜的分子事件，mRNA、lncRNA和circRNA可通過競爭性結合miRNA形成復雜的ceRNA調控網絡。前人研究發現ame-miR-6001-3p在蜜蜂蛻皮激素分泌[8]和級型分化[9]過程中具有調控作用。意蜂工蜂中腸發育過程中分別有29個顯著性上調lncRNA，以及14個顯著性上調circRNA參與競爭性結合ame-miR-6001- 3p[10-11]。本研究以顯著下調表達的ame-miR-6001-3p為介導，構建DEG與ame-miR-6001-3p的調控網絡，發現ame-miR-6001-3p靶向結合54個上調基因和44個下調基因。根據ceRNA假說，推測lncRNA和circRNA通過改變表達水平影響對ame-miR-6001-3p的吸附結合，從而調節ame-miR-6001-3p的表達水平，進而部分解除ame-miR-6001-3p對靶mRNA的調控。進一步對ame-miR-6001-3p靶向結合的上調基因進行功能及代謝通路注釋，發現它們主要涉及胰島素信號通路（2個上調基因）、Hippo信號通路（2個上調基因）、PI3K-Akt信號通路（1個上調基因）等多條信號通路，以及磷酸肌醇代謝（2個上調基因）和谷胱甘肽代謝（1個上調基因）等各類代謝途徑。綜上，推測ceRNA網絡在意蜂工蜂中腸發育過程中發揮廣泛的調控作用。下一步的工作重點是結合意蜂工蜂中腸的miRNA組學數據構建mRNA-miRNA-lncRNA- circRNA調控網絡，進一步篩選出關鍵代謝通路、關鍵基因和關鍵ncRNA作為候選靶標，并通過過表達和敲減進行功能研究。

4 結論

意蜂工蜂中腸發育過程中DEG通過參與TGF- beta、Wnt、Hippo、Notch、PI3K-Akt、mTOR、AMPK和NF-B等信號通路影響中腸的生長發育和免疫防御；ceRNA網絡在中腸發育中發揮廣泛的調控作用；顯著下調表達的ame-miR-6001-3p靶向結合的54個上調基因可參與胰島素信號通路等多條代謝途徑。

[1] COON K L, VOGEL K J, BROWN M R, STRAND M R. Mosquitoes rely on their gut microbiota for development., 2014, 23(11): 2727-2739.

[2] CHOUAIA B, ROSSI P, MONTAGNA M, RICCI I, CROTTI E, DAMIANI C, EPIS S, FAYE I, SAGNON N, ALMA A, FAVIA G, DAFFONCHIO D, BANDI C. Molecular evidence for multiple infections as revealed by typing ofbacterial symbionts of four mosquito species., 2010, 76(22): 7444-7450.

[3] STORELLI G, DEFAYE A, ERKOSAR B, HOLS P, ROYET J, LEULIER F.promotessystemic growth by modulating hormonal signals through TOR-dependent nutrient sensing., 2011, 14(3): 403-414.

[4] ZHENG H, POWELL J E, STEELE M I, DIETRICH C, MORAN N A. Honeybee gut microbiota promotes host weight gain via bacterial metabolism and hormonal signaling., 2017, 114(18): 4775-4780.

[5] 馮倩倩. 意大利蜜蜂春繁階段代用花粉中VA、VE適宜水平的研究[D]. 泰安: 山東農業大學, 2012.

FENG Q Q. Optimal level of VA and VE in pollen substitutes ofduring the period of spring multiplication[D]. Taian: Shandong Agricultural University, 2012. (in Chinese)

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

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

[7] 郭睿, 杜宇, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 陳華枝, 耿四海, 周丁丁, 石彩云, 趙紅霞, 陳大福. 意大利蜜蜂幼蟲腸道發育過程中的差異表達microRNA及其調控網絡. 中國農業科學, 2018, 51(21): 4197-4209.

GUO R, DU Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, CHEN H Z, GENG S H, ZHOU D D, SHI C Y, ZHAO H X, CHEN D F. Differentially expressed microRNA and their regulation networks during the developmental process oflarval gut., 2018, 51(21): 4197-4209. (in Chinese)

[8] MELLO T R P, ALEIXO A C, PINHEIRO D G, NUNES F M F, BITONDI M M G, HARTFELDER K, BARCHUK A R, SIM?ES Z L P. Developmental regulation of ecdysone receptor (EcR) and EcR-controlled gene expression during pharate-adult development of honeybees ()., 2014, 5: 445.

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

[10] 郭睿, 耿四海, 熊翠玲, 鄭燕珍, 付中民, 王海朋, 杜宇, 童新宇, 趙紅霞, 陳大福. 意大利蜜蜂工蜂中腸發育過程中長鏈非編碼RNA的差異表達分析. 中國農業科學, 2018, 51(18): 3600-3613.

GUO R, GENG S H, XIONG C L, ZHENG Y Z, FU Z M, WANG H P, DU Y, TONG X Y, ZHAO H X, CHEN D F. Differential expression analysis of long non-coding RNAs during the developmental process ofworker’s midgut., 2018, 51(18): 3600-3613. (in Chinese)

[11] 郭睿, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 杜宇, 王海朋, 耿四海, 周丁丁, 劉思亞, 陳大福. 意大利蜜蜂工蜂中腸發育過程中的差異表達環狀RNA及其調控網絡分析. 中國農業科學, 2018, 51(23): 4575-4590.

GUO R, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, XU G J, DU Y, WANG H P, GENG S H, ZHOU D D, LIU S Y, CHEN D F. Analysis of differentially expressed circular RNAs and their regulation networks during the developmental process ofworker’s midgut., 2018, 51(23): 4575-4590. (in Chinese)

[12] HUANG Q, CHEN Y P, WANG R W, CHENG S, EVANS J D. Host-parasite interactions and purifying selection in a microsporidian parasite of honey bees., 2016, 11(2): e0147549.

[13] HUANG W F, SOLTER L F. Comparative development and tissue tropism ofand., 2013, 113(1): 35-41.

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

[15] 郭睿, 杜宇, 童新宇, 熊翠玲, 鄭燕珍, 徐國鈞, 王海朋, 耿四海, 周丁丁, 郭意龍, 吳素珍, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌侵染前期的差異表達microRNA及其調控網絡. 中國農業科學, 2019, 52(1): 166-180.

GUO R, DU Y, TONG X Y, XIONG C L, ZHENG Y Z, XU G J, WANG H P, GENG S H, ZHOU D D, GUO Y L, WU S Z, CHEN D F. Differentially expressed microRNAs and their regulation networks inlarval gut during the early stage ofinfection., 2019, 52(1): 166-180.(in Chinese)

[16] MARTINSON V G, MOY J, MORAN N A. Establishment of characteristic gut bacteria during development of the honeybee worker., 2012, 78(8): 2830-2840.

[17] KWONG W K, MANCENIDO A L, MORAN N A. Immune system stimulation by the native gut microbiota of honey bees., 2017, 4(2): 170003.

[18] CHEN D, GUO R, XU X, XIONG C, LIANG Q, ZHENG Y, LUO Q, ZHANG Z, HUANG Z, KUMAR D, XI W, ZOU X, LIU M. Uncovering the immune responses oflarval gut toinfection utilizing transcriptome sequencing., 2017, 621: 40-50.

[19] 郭睿, 杜宇, 周倪紅, 劉思亞, 熊翠玲, 鄭燕珍, 付中民, 徐國鈞, 王海朋, 耿四海, 周丁丁, 陳大福. 意大利蜜蜂幼蟲腸道在球囊菌脅迫后期的差異表達微小RNA及其靶基因分析. 昆蟲學報, 2019, 62(1): 49-60.

GUO R, DU Y, ZHOU N H, LIU S Y, XIONG C L, ZHENG Y Z, FU Z M, XU G J, WANG H P, GENG S H, ZHOU D D, CHEN D F. Comprehensive analysis of differentially expressed microRNAs and their target genes in the larval gut ofduring the late stage ofstress, 2019, 62(1): 49-60.(in Chinese)

[20] BARRY E R, CAMARGO F D. The Hippo superhighway: signaling crossroads converging on the Hippo/Yap pathway in stem cells and development., 2013, 25(2): 247-253.

[21] 陳曉. 蜜蜂卵巢激活和產卵過程差異表達的編碼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)

[22] 寧越, 米雪, 陳星伊, 邵建航, 昝林森. SMAD1基因的沉默和過表達及對秦川牛原代成肌細胞生肌的影響. 中國農業科學, 2019, 52(10): 1818-1829.

NING Y, MI X, CHEN X Y, SHAO J H, ZAN L S. Silencing and overexpressing SMAD family member 1 () gene and its effect on myogenesis in primary myoblast of Qinchuan cattle ()., 2019, 52(10): 1818-1829.(in Chinese)

[23] BRUMMEL T, ABDOLLAH S, HAERRY T E, SHIMELL M J, MERRIAM J, RAFTERY L, WRANA J L, O’CONNOR M B. Theactivin receptor baboon signals through dSmad2 and controls cell proliferation but not patterning during larval development., 1999, 13(1): 98-111.

[24] 張振, 盧忠燕. 家蠶類胰島素信號通路及其功能的研究進展. 蠶學通訊, 2014, 34(2): 22-26.

ZHANG Z, LU Z Y. Research progress on insulin signaling pathway and its function in silkworm., 2014, 34(2): 22-26. (in Chinese)

[25] GAUTAM S, ISHRAT N, YADAV P, SINGH R, NARENDER T, SRIVASTAVA A K. 4-Hydroxyisoleucine attenuates the inflammation- mediated insulin resistance by the activation of AMPK and suppression of SOCS-3 coimmunoprecipitation with both the IR-subunit as well as IRS-1., 2016, 414(1/2): 95-104.

[26] NEUFELD T P. Body building: regulation of shape and size by PI3K/TOR signaling during development., 2003, 120(11): 1283-1296.

[27] KAUR S, SASSANO A, DOLNIAK B, JOSHI S, MAJCHRZAK- KITA B, BAKER D P, HAY N, FISH E N, PLATANIAS L C. Role of the Akt pathway in mRNA translation of interferon-stimulated genes., 2008, 105(12): 4808-4813.

[28] LIU S, LUCAS K J, ROY S, HA J, RAIKHEL A S.Mosquito-specific microRNA-1174 targets serine hydroxymethyltransferase to control key functions in the gut., 2014, 111(40): 14460-14465.

[29] 熊翠玲, 杜宇, 陳大福, 鄭燕珍, 付中民, 王海朋, 耿四海, 陳華枝, 周丁丁, 吳素珍, 石彩云, 郭睿. 意大利蜜蜂幼蟲腸道的miRNAs的生物信息學預測及分析. 應用昆蟲學報, 2018, 55(6): 1023-1033.

XIONG C L, DU Y, CHEN D F, ZHENG Y Z, FU Z M, WANG H P, GENG S H, CHEN H Z, ZHOU D D, WU S Z, SHI C Y, GUO R. Bioinformatic prediction and analysis of miRNAs in thelarval gut., 2018, 55(6): 1023-1033.(in Chinese)

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

Profiling and Regulation Network of Differentially Expressed Genes During the Development Process ofWorker’s Midgut

DU Yu, ZHOU DingDing, WAN JieQi, LU JiaXuan, FAN XiaoXue, FAN YuanChan, CHEN Heng, XIONG CuiLing, ZHENG YanZhen, FU ZhongMin, XU GuoJun, CHEN DaFu, GUO Rui

(College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002)

【】The whole transcriptome sequencing of7- and 10-day-old workers’ midguts (Am7 and Am10) was previously conducted. In this study, the differential expression profile and regulation network of genes were investigated to reveal the molecular mechanism underlying the midgut development.【】Gene expressions were calculated based on FPKM (fragments per kilobase of transcript per million mapped reads) algorithm, and differentially expressed genes (DEGs) were gained following the standard |log2fold change|≥1 and≤0.05. Target mRNAs of ame-miR-6001-3p were predicted utilizing TargetFinder. Annotations of all DEGs in GO and KEGG databases were performed using related software. In addition, DEGs enriched in 13 signaling pathways including AMPK, P13K-Akt, Wnt, cAMP, FoxO, Hippo, mTOR, Jak-STAT, Toll-like receptor, TGF-beta, Notch, MAPK and NF-B, as well as DEGs targeted by ame-miR-6001-3p were screened out, followed by visualization of enrichment networks and regulation networks with Cytoscape. Finally, Stem loop RT-PCR and RT-qPCR were used to verify the expression and differential expression pattern of ame-miR-6001-3p and DEGs in Am7 and Am10.【】A total of 1 038 DEGs were identified in Am7 vs Am10 comparison group, including 515 up- and 523 down-regulated genes, respectively. These DEGs were associated with cellular process, metabolic process and catalytic activity, and significantly enriched in some material and energy metabolisms such as oxidative phosphorylation, amino sugar and nucleotide sugar metabolisms, fatty acid metabolism and purine metabolism, indicative of the active cellular and metabolic activities. Expression cluster analysis suggested that 20, 18, 15 and 14 DEGs were respectively enriched in AMPK signaling pathway, PI3K-Akt signaling pathway, endocytosis and Hippo signaling pathway. In addition, 57 DEGs were enriched in the aforementioned 13 signaling pathways associated with growth and development as well as immune defense, among them one DEG was enriched in several signaling pathways. Moreover, regulation network analysis showed that 54 up-regulated genes and 44 down-regulated genes were targets of ame-miR-6001-3p, respectively; these up-regulated genes were enriched in 43 pathways including inositol phosphate metabolism, Hippo signaling pathway, glutathione metabolism and insulin signaling pathway, while these down-regulated genes were enriched in 20 pathways including Hippo signaling pathway, metabolic pathways, glutathione metabolism and arachidonic acid metabolism. Moreover, RT-qPCR result showed that the variation trend of six randomly selected DEGs were consistent with that in sequencing data, confirming the reliability of DEGs. Finally, ame-miR-6001-3p was definitely expressed and significantly down-regulated in Am10.【】In this work, the expression pattern of DEGs and the regulation network between DEGs and ame-miR-6001-3p as well as the potential role of DEGs during the developmental process ofworker’s midgut were deeply analyzed. The results revealed that the DEGs may participate in various signaling pathways including TGF-beta, Wnt, Hippo, Notch, PI3K-Akt, mTOR, AMPK, NF-B signaling pathways, thus affecting the growth and development as well as immune defense of the midgut; DEGs were likely to regulate several signaling pathways such as insulin signaling pathway during the midgut development via formation regulation networks with significantly down-regulated ame-miR-6001-3p.

; midgut; development; differentially expressed gene (DEG); competitive endogenous RNA; regulation network

10.3864/j.issn.0578-1752.2020.01.019

2019-07-09；

2019-09-02

國家現代農業產業技術體系建設專項資金（CARS-44-KXJ7）、福建省自然科學基金（2018J05042）、福建省教育廳中青年教師教育科研項目（JAT170158）、福建農林大學杰出青年科研人才計劃（xjq201814）、福建農林大學科技創新專項基金（CXZX2017342，CXZX2017343）、福建省大學生創新創業訓練計劃（3165602032，3155006018）

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

（責任編輯 岳梅）