三种微小RNA在意大利蜜蜂工蜂蛹期发育过程中的表达谱及潜在功能

祝智威, 付中民, 隆 琦, 杜 宇, 张文德, 胡 颖,赵 萧, 史小玉, 徐细建, 陈大福,2,*, 郭 睿,2,*

(1.福建农林大学动物科学学院(蜂学学院), 福州 350002; 2.福建农林大学蜂疗研究所, 福州 350002;3.江西省养蜂研究所, 南昌 330000)

意大利蜜蜂Apismelliferaligustica是我国养蜂生产中的主要蜂种，其工蜂在12 d的封盖子蛹期历经组织重构、器官发育及体色和眼色加深等关键的变态发育过程(Soaresetal., 2007)。研究表明温度等外界因素可显著影响封盖子的正常发育、成活率、羽化后的短期记忆行为及社会性行为(Tautzetal., 2003; Jonesetal., 2005; 姚丹等, 2020)。

微小RNA(microRNA, miRNA)是一类长度约为19～25个核苷酸的内源性非编码RNA(non-coding RNA, ncRNA)，广泛参与基因的转录后调控，进而影响细胞增殖、新陈代谢、组织发育、器官形成和免疫应答等生物学过程(Bartel, 2004; He and Hannon, 2004; Friedmanetal., 2009)。2006年完成西方蜜蜂Apismellifera的基因组测序并公布(Honey Bee Genome Sequencing Consortium, 2006)，随后miRNA逐渐成为蜜蜂分子生物学研究领域的热点之一(Liuetal., 2012; Ashbyetal., 2016; Zhuetal., 2017)。此前，ame-miR-13b, ame-miR-100和ame-miR-bantam在蜜蜂中的研究已见诸报道(Herranzetal., 2012; Liuetal., 2012; Lozanoetal., 2015; Shietal., 2015; Ashbyetal., 2016; Macedoetal., 2016)，但总体仍比较有限。Liu等(2012)曾对西方蜜蜂哺育蜂和采集蜂头部的miRNA进行检测，发现ame-miR-13b在采集蜂中上调表达，推测ame-miR-13b可能参与调控西方蜜蜂的职能转变。Macedo等(2016)通过实时荧光定量PCR检测发现，相比于有活性的意大利蜜蜂工蜂卵巢和交配后的蜂王卵巢，ame-miR-13b在无活性的西方蜜蜂工蜂卵巢和处女王中表达水平更高，作者推测ame-miR-13b与西方蜜蜂的级型分化有关。ame-miR-100被证实与蜕皮激素(ecdysteroid, Ec)和保幼激素(juvenile hormone, JH)的效应具有相关性，可参与调节蜜蜂的级型分化和蜂王的繁殖过程(Herranzetal., 2012; Ashbyetal., 2016; Macedoetal., 2016)。Ashby等(2016)研究发现miR-bantam可参与蜜蜂的细胞增殖与分化以及组织结构重塑等过程。此外，miR-13b, miR-100和miR-bantam能够影响德国小蠊Blattellagermanica和黑腹果蝇Drosophilamelanogaster的变态发育(Bashirullahetal., 2003; Lozanoetal., 2015)以及果蝇的翅发育(Becametal., 2011)。蜜蜂是一种完全变态昆虫，一生历经幼虫、蛹和成虫3个发育阶段，其中蛹期内部组织和器官发生剧烈地消解和重组(梁勤和陈大福, 2009)。然而，蜜蜂蛹期发育相关的miRNA研究缺失，ame-miR-13b, ame-miR-100和ame-miR-bantam在蜜蜂蛹期发育中的调控作用迄今未明。

本研究利用分子生物学手段对ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂蛹期不同日龄封盖子中的真实表达进行验证，并检测三者在工蜂蛹期发育过程的表达谱，然后通过生物信息学方法预测上述3个miRNA的靶mRNA并分析miRNA-mRNA调控网络，以期为探明意大利蜜蜂工蜂蛹期发育的分子调控机理提供理论和实验依据。

1 材料与方法

1.1 供试生物材料

意大利蜜蜂封盖子取自福建农林大学动物科学学院(蜂学学院)教学蜂场蜂群。选取群势较强的3群意大利蜜蜂，将蜂王限制放在空脾上产卵8 h，随后将蜂王取出，不让蜂王再次进入该巢脾产卵。待蜂子发育至封盖后第3天(白眼蛹)，将巢脾提取至实验室，置于恒温恒湿培养箱内(35±0.5℃，RH 70%)。用镊子小心将蛹取出，并立即投入液氮速冻，每3头放入1个RNA-Free的EP管，随后放入-80℃进行保存。白眼蛹记为P1，此后每24 h连续取样至羽化前1 d，P2, P3, P4, P5, P6, P7和P8分别代表每24 h取得的样品。每个样品设置3个生物学重复。

1.2 miRNA的stem-loop RT-PCR和Sanger测序验证

参照前期已建立的方法(Chenetal., 2005; 杜宇等, 2019)，通过stem-loop RT-PCR验证ame-miR-13b, ame-miR-100和ame-miR-bantam以及内参基因actin-1在意大利蜜蜂工蜂蛹期的表达情况。利用Primer Premier 6软件分别设计miRNA的stem-loop引物、特异性上游引物和通用下游引物以及actin-1基因的上下游引物(表1)，委托生工生物工程(上海)股份有限公司合成引物。利用RNA抽提试剂盒(Axygen公司，美国)提取意大利蜜蜂工蜂蛹期(P1-P8)的总RNA，利用miRNA 1st Strand cDNA Synthesis Kit(by Stem-loop)(诺唯赞，南京)和stem-loop引物进行反转录合成cDNA模板，用于3种miRNA的PCR反应，以无菌水为模板作为阴性对照，RT-PCR反应体系(20 μL): PCR Mix 10 μL, DEPC水 7 μL, cDNA模板1 μL, 上下游引物(10 μmol/L)各1 μL。反应程序: 94℃ 5 min; 94℃ 30 s, 56℃ 30 s, 72℃ 1 min, 共33个循环；最后72℃5 min。PCR产物经1.8%琼脂糖凝胶电泳检测。将目的片段(大小约为70 bp)切胶回收，利用零背景TOPO-TA克隆试剂盒(YEASEN，上海)按照说明书连接pESI-T载体；然后将连接产物转化大肠杆菌EscherichiacoliDH5α感受态细胞(唯地，上海)，37℃ 180 r/min振荡培养30 min。4 000 r/min离心1 min后弃上清，保留80 μL上清重悬沉淀，均匀涂布于含氨苄青霉素(100 μg/mL)的LB固体培养基上，37℃培养8 h，挑斑摇菌后进行菌液PCR，将阳性结果对应的菌液送生工生物工程(上海)股份有限公司单端测序。

1.3 miRNA的靶mRNA的预测与分析

参照笔者所在实验室前期建立的技术流程(熊翠玲等, 2018b; 郭睿等, 2019)，联用RNAhybrid, miRanda和TargetScan软件预测ame-miR-13b, ame-miR-100和ame-miR-bantam的靶mRNA，采用软件默认参数。然后取3种软件预测结果的交集作为可靠的靶mRNA集合。进而采用BLAST软件将上述靶mRNA比对GO数据库(http:∥geneontology.org/)和KEGG数据库(https:∥www.kegg.jp/)，从而得到相应的功能和通路注释。

1.4 miRNA和靶mRNA的RT-qPCR检测

利用RNA抽提试剂盒(Promega公司，美国)分别提取意大利蜜蜂工蜂蛹期(P1-P8)不同日龄封盖子的总RNA，利用miRNA 1st Strand cDNA Synthesis Kit(by Stem-loop)(诺唯赞，南京)进行反转录合成相应的cDNA模板，用于3种miRNA及内参actin-1的RT-qPCR；利用反转录试剂盒(YEASEN，上海)反转录合成P1-P8的cDNA模板，用于靶mRNA及内参actin-2的RT-qPCR检测。采用Ct值比较法进行基因相对定量，利用ABI QuantStudio 3荧光定量PCR系统(ABI 公司，美国)对miRNA和1.3节预测的靶mRNA XM_016916453.1进行RT-qPCR检测。反应体系(20 μL): SYBR Green Dye 10 μL, 1.2节设计的上下游引物(10 μmol/L)各1 μL, cDNA模板 1 μL, DEPC水补至20 μL。miRNA的反应条件: 95℃ 5 min; 95℃ 30 s, 56℃ 30 s, 共40个循环。mRNA的反应条件: 95℃ 5 min; 95℃变性15 s, 56℃ 30 s, 72℃ 20 s, 共40个循环，熔解曲线程序均默认系统设置。将上述3个miRNA及靶mRNA XM_016916453.1在P1时期的表达量作为对照组，在其他日龄的miRNA及靶mRNA的相对表达量采用2-ΔΔCt法进行计算。每个反应进行3次平行重复和3次技术重复。

1.5 miRNA及靶mRNA的调控网络构建与分析

参照范小雪等(2021)的方法，从1.3节中的靶向预测结果中进一步筛选得到结合自由能小于-20 kcal/mol的靶向结合关系，并据此构建miRNA-mRNA调控网络，然后利用Cytoscape软件实现调控网络的可视化。

1.6 数据分析

通过GraphPad Prism 8软件绘制miRNA和靶mRNA表达趋势的线形图。采用SPSS软件对3个miRNA在P1-P8蛹期期间的表达量进行单因素方差分析(ANOVA)，以P<0.05为显著性阈值，利用Tukey氏检验法分析实验数据的差异显著性。

2 结果

2.1 意大利蜜蜂工蜂蛹期的体色和眼色变化

意大利蜜蜂工蜂蛹期的体色和眼色逐渐发生变化：P1眼色与体色均为白色；P2眼色为粉红色，体色为白色；P3眼色明显加深变为红色，体色变为淡黄色；P4眼色继续加深变为深红色，体色变为黄色；P5-P8眼色进一步加深至黑色，体色趋于深褐色。

2.2 ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂工蜂蛹期的表达和序列验证

RT-PCR产物的琼脂糖凝胶电泳结果显示，ame-miR-13b, ame-miR-100和ame-miR-bantam在P1-P8中均能扩增出符合预期大小(约70 bp)的目的片段，证实了ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂工蜂蛹期的真实表达(图2: A-C)。进一步的Sanger测序结果显示扩增片段与相应的测序序列一致(图2: E-G)，表明上述3个miRNA真实存在。

图2 意大利蜜蜂工蜂蛹期3种miRNA的stem-loop RT-PCR及Sanger测序验证

2.3 3个miRNA的靶mRNA预测及其功能注释

靶向预测结果显示，ame-miR-13b, ame-miR-100和ame-miR-bantam分别靶向结合850, 136和506条靶mRNA。GO数据库注释结果显示，上述3个miRNA的靶mRNA均涉及生物学进程、细胞组分和分子功能三大类的多个条目，其中ame-miR-13b靶向的mRNA可注释到32个条目，在生物学进程大类注释靶mRNA数量最多的条目是细胞进程(281条)、单一组织进程(229条)、代谢进程(180条)、应激反应(86条)和生物学调控(82条)；在细胞组分大类注释靶mRNA数量最多的条目是细胞(110条)，细胞部分(110条)，细胞膜(107条)，细胞膜部分(107条)，细胞器(98条)；在分子功能大类注释靶mRNA数量最多的条目是结合(294条)，催化活性(153条)，分子传感器活性(54条)，信号传感器活性(47条)和转运子活性(45条)(图3: A)。ame-miR-100靶向的mRNA可注释到24个条目，在生物学进程大类中富集靶mRNA数量最多条目的是细胞进程(36条)，代谢进程(30条)，单一组织进程(27条)，多细胞生物进程(16条)，信号(7条)；在细胞组分大类中富集靶mRNA数量最多的条目是细胞(28条)，细胞部分(28条)，细胞器(21条)，细胞膜(5条)，细胞膜部分(5条)；在分子功能大类中富集靶mRNA数量最多的条目是结合(35条)，核酸转录因子活性(16条)，催化活性(6条)，转运子活性(6条)，信号传感器活性(4条)(图3: B)。ame-miR-bantam靶向的mRNA可注释到33个条目，在生物学进程大类中富集靶mRNA数量最多的条目是细胞进程(116条)，代谢进程(113条)，单组织进程(95)，应激反应(42条)，信号(41条)；在细胞组分大类中富集靶mRNA数量最多的条目是细胞(69条)，细胞部分(69条)，细胞器(62条)，细胞膜(44条)，细胞膜部分(43条)；在分子功能大类中富集靶mRNA数量最多的条目是结合(154条)，催化活性(55条)，核酸结合转录因子活性(24条)，分子传感器活性(23条)，信号传感器活性(21条)(图3: C)。

图3 意大利蜜蜂工蜂蛹期ame-miR-13b(A), ame-miR-100(B)和ame-miR-bantam(C)靶mRNA的GO数据库注释

KEGG数据库注释结果显示，ame-miR-13b的靶mRNA可注释到204条通路，包括代谢通路(56条)、嘌呤代谢(19条)、自噬-动物(12条)、神经活性的配体-受体相互作用(11条)和背腹轴的形成(8条)；ame-miR-100的靶mRNA可注释到114条通路，包括FoxO信号通路(9条)、mTOR信号通路(7条)、Wnt信号通路(7条)、Hippo信号通路-果蝇(6条)和MAPK信号通路(5条)；ame-miR-bantam的靶mRNA可注释到171条通路，包括代谢通路(30条)、Rap1信号通路(16条)、mRNA监测通路(15条)、Wnt信号通路(12条)和cAMP信号通路(9条)。

2.4 ame-miR-13b, ame-miR-100和ame-miR-bantam及共同靶mRNA在意大利蜜蜂工蜂蛹期发育过程的表达谱

RT-qPCR结果表明，ame-miR-13b在P1-P8阶段总体表达趋势为上调；在P1-P2阶段上调表达，P2-P4阶段略有下调，P4-P7阶段持续上调，P7-P8阶段转而下调；P7时表达量达到最高(图4: A)。ame-miR-100在P1-P8阶段总体呈上升趋势；表达水平在P1-P3阶段上升，P3-P4阶段下降，P4-P8阶段持续上升；在P8时的表达水平达到峰值(图4: B)。ame-miR-bantam在P1-P8阶段总体上也呈上升趋势；在P1-P8阶段不同时间点的表达趋势与ame-miR-100基本一致(图4: C)。此外，上述3个miRNA共同的靶mRNA XM_016916453.1在P1-P8阶段总体表现为下调趋势；在P7时表达量最高，P8时表达量最低(图4: D)。

图4 意大利蜜蜂工蜂蛹期ame-miR-13b(A), ame-miR-100(B)和ame-miR-bantam(C)及靶mRNA XM_016916453.1(D)的RT-qPCR验证结果

2.5 意大利蜜蜂蛹期miRNA-mRNA调控网络

ame-miR-13b, ame-miR-100和ame-miR-bantam与靶mRNA之间形成较为复杂的调控网络(图5)。分析结果显示ame-miR-bantam靶向结合194条mRNA，数量最多；ame-miR-100靶向结合60条mRNA；ame-miR-13b靶向结合137条mRNA。进一步分析发现，ame-miR-13b和ame-miR-bantam共同靶向调控NM_001011629.1, XM_006568833.2, XM_006571679.2, XM_006571681.2, XM_016915393.1和XM_016916690.1；ame-miR-13b, ame-miR-100和ame-miR-bantam共同靶向调控XM_006563323.2, XM_006563325.2, XM_016916453.1和XM_016916455.1；此外，ame-miR-13b可同时靶向结合4个蜕皮激素相关mRNA(NM_001098215.2, NM_001159355.1, XM_006568053.2和XM_016913298.1)(图5)。

图5 意大利蜜蜂工蜂蛹期ame-miR-13b, ame-miR-100和ame-miR-bantam及相应靶mRNA间的调控网络

3 讨论

miR-13b, miR-100和miR-bantam已被研究证实参与调控蜜蜂(Ashbyetal., 2016)、家蚕Bombyxmori(Lingetal., 2015)和果蝇(Aboobakeretal., 2005)等昆虫的变态发育过程。但意大利蜜蜂工蜂蛹期变态发育过程中ame-miR-13b, ame-miR-100和ame-miR-bantam的调控作用迄今不明。本研究中，我们在实验室条件下模拟自然蜂群的内部环境对意大利蜜蜂工蜂蛹进行人工培养，发现P1-P8阶段各日龄蛹的体色和眼色变化规律与自然蜂群巢房中工蜂蛹的体色和眼色变化规律一致(图1)，证明了通过人工培养的方法将意大利蜜蜂工蜂幼虫培养至蛹期具有可行性，为将来通过饲喂模拟物(mimic)和抑制物(inhibitor)探究蛹期发育中miRNA的分子功能提供了实验依据。Stem-loop RT-PCR结果显示ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂工蜂蛹变态发育阶段(P1-P8)真实表达(图2)，暗示三者对蛹期发育具有潜在的重要调控作用。Sanger测序结果显示3个miRNA的序列与miRBase数据库中收录的miRNA序列一致，证实了ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂工蜂蛹中真实存在。此外，RT-qPCR结果显示上述3个miRNA在意大利蜜蜂工蜂蛹变态发育过程总体上都呈现上调表达趋势(图4)，进一步说明三者在蛹变态发育过程的重要性。

靶向预测结果显示，ame-miR-13b, ame-miR-100和ame-miR-bantam分别靶向调控850, 136和506条靶mRNA，这些靶mRNA均可注释到细胞进程、单一组织进程、代谢进程、细胞、细胞部分和结合等GO条目(图3)，说明上述3个miRNA潜在广泛调控意大利蜜蜂工蜂蛹期的细胞生命活动和新陈代谢过程。Hippo信号通路能够调节昆虫器官的生长发育，并与影响细胞分化、细胞凋亡、细胞增殖等过程的Notch信号通路相互作用，进而影响果蝇的眼、翅及其他器官的分化过程(Halder and Johnson, 2011; Kimetal., 2011; Kalkman, 2012; Barry and Camargo, 2013)。此外，Hippo信号通路还与Wnt, FoxO, mTOR和Jak-STAT等信号通路共同参与调节昆虫的免疫应答、肠道发育及卵巢发育等过程(Herranzetal., 2012; 郭睿等, 2018; 杜宇等, 2019)。本研究中，KEGG数据库注释结果显示上述靶mRNA涉及Hippo, FoxO, Notch, Wnt和Jak-STAT等信号通路；其中ame-miR-13b的5条靶mRNA显著富集在Jak-STAT信号通路；ame-miR-100的靶mRNA显著富集在FoxO(9条)、mTOR(7条)、Wnt(7条)、Jak-STAT(4条)及Notch(3条)信号通路；ame-miR-bantam的靶mRNA显著富集在Wnt(12条)、FoxO(8条)及Notch(7条)信号通路。以上结果表明上述3个miRNA调控相关靶基因表达对上述重要的生长发育信号通路进行调节，进而影响意大利蜜蜂工蜂蛹期的变态发育过程。

蜕皮激素(Ec)和保幼激素(JH)的共同作用可激活数十个基因的表达进而参与昆虫的变态调节(Margametal., 2006)。笔者所在团队前期通过对意大利蜜蜂4, 5和6日龄幼虫肠道进行small RNA-seq(sRNA-seq)和生物信息学分析，揭示ame-miR-13b和ame-miR-2等136个差异表达miRNA潜在参与调控意大利蜜蜂幼虫的Ec滴度，使其与JH滴度保持动态平衡，从而保证幼虫的正常发育过程(熊翠玲等, 2018a)。miR-13b作为无脊椎动物特有的miR-2家族的成员之一，参与德国小蠊B.germanica的变态发育过程，并调节家蚕成虫翅发育，此外还与蜜蜂级型分化过程密切相关。本研究发现，ame-miR-13b可同时靶向结合4个与Ec相关的mRNA(NM_001098215.2, NM_001159355.1, XM_006568053.2和XM_016913298.1)；此外RT-qPCR结果显示ame-miR-13b的表达水平在P2-P8阶段呈动态上升趋势(图4: A)。推测ame-miR-13b通过靶向蜕皮激素相关基因的表达，间接影响Ec的滴度，保证意大利蜜蜂蛹期的正常变态发育过程。背后的分子机理值得进一步深入研究。

黑色素与其他色素共同作用可形成昆虫丰富的体色和斑纹(Wright, 1987)。yellow基因编码的蛋白可与多巴以及酚氧化酶的共同作用形成多巴黑色素，进而参与昆虫的体色形成过程(崔坤蓉等, 2015)。前人研究发现yellow蛋白家族与调节蜜蜂级型分化以及生长发育密切相关的王浆主蛋白家族(MRJP)之间具有较近的进化关系(李艳, 2007; 赵亚周等, 2012)。本研究发现，ame-miR-13b和ame-miR-bantam共同靶向yellow的mRNA(GenBank登录号: XM_006564945.2)，说明二者潜在参与蜜蜂蛹期的体色变化及生长发育等过程。未来可根据本研究调控网络中的miRNA-mRNA靶向关系，进一步通过过表达和敲减技术深入探究miRNA调控意大利蜜蜂蛹期变态发育及其体色变化的分子机理。

综上所述，本研究结合分子生物学和生物信息学方法对ame-miR-13b, ame-miR-100和ame-miR-bantam在意大利蜜蜂工蜂蛹期(P1-P8)发育过程的表达谱及调控网络进行分析和探讨，结果表明上述3个miRNA的表达水平在蛹期发育过程中均呈动态上升趋势，并通过调控蜕皮激素相关mRNA和yellow基因的mRNA水平表达以及Hippo, FoxO, Wnt, Jak-STAT和Notch等信号通路潜在参与调节意大利蜜蜂工蜂蛹期的变态发育过程。研究结果为进一步深入探究意大利蜜蜂蛹期变态发育的分子机理提供了理论和实验依据。