杜宇，周丁丁，万洁琦，卢家轩，范小雪，范元婵，陈恒，熊翠玲，郑燕珍，付中民，徐国钧，陈大福，郭睿

意大利蜜蜂工蜂中肠发育过程中的差异基因表达谱及调控网络

杜宇，周丁丁，万洁琦，卢家轩，范小雪，范元婵，陈恒，熊翠玲，郑燕珍，付中民，徐国钧，陈大福，郭睿

（福建农林大学动物科学学院（蜂学学院），福州 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

（责任编辑 岳梅）