东方蜜蜂微孢子虫胁迫下意大利蜜蜂工蜂的piRNA差异表达谱及潜在功能

许雅静， 孙明会， 刘佳美， 郭意龙， 胡 颖， 张佳欣，赵 萧， 张凯遥， 熊翠玲， 陈大福,2， 郭 睿,2

(1.福建农林大学动物科学学院 蜂学学院， 福州 350002； 2.福建农林大学蜂疗研究所， 福州 350002)

1 引 言

意大利蜜蜂(Apismelliferaligustica)作为西方蜜蜂(Apismellifera)的亚种之一，在我国和其他国家的养蜂生产中广泛使用，具有重要的生态和经济价值[1].蜜蜂是一种真社会性昆虫，因群居性而易遭受细菌、真菌和病毒等病原微生物的侵袭[1].东方蜜蜂微孢子虫(Nosemaceranae)是一种专性侵染成年蜜蜂中肠上皮细胞的真菌病原，对蜂王、工蜂、雄蜂及蜜蜂幼虫均具有侵染性，导致的蜜蜂微孢子虫病严重危害蜜蜂健康并影响养蜂业的可持续发展[2].东方蜜蜂微孢子虫侵染可对蜜蜂宿主造成诸多不良影响，包括能量胁迫、免疫抑制、细胞凋亡抑制、中肠上皮细胞结构破坏及寿命缩短等[3].

非编码RNA(non-coding RNA, ncRNA)起初被认为不具有蛋白编码能力和生物学功能，是基因转录过程的“噪音”和“副产物”，但随着技术手段的不断进步和相关研究的持续深入，ncRNA被证实在许多生物学过程中起到关键作用[4].根据长度，ncRNA可分为长链非编码RNA(long non-coding RNA, lncRNA)和小ncRNA，后者又包括微小RNA(microRNA, miRNA)、小干扰RNA(Small interfering RNA, siRNA)和PIWI蛋白互作RNA(PIWI-interacting RNA, piRNA)等.成熟piRNA的长度一般介于26～30 nt，与miRNA(20～24 nt)和siRNA(21～25 nt)的长度相近[5].在真核生物中，piRNA能参与控制转座子的活性及维持基因组的稳定性[6].此外，有研究表明piRNA通过与靶mRNA的碱基互补配对发挥基因表达调控作用[7].PIWI蛋白最早在果蝇(Drosophilamelanogaster)体内被发现，并被证实参与了果蝇种系干细胞的维持和自我更新[8].随着深度测序技术的兴起与发展，人们在哺乳动物[9]、昆虫[10]和鱼类[11]等物种中均陆续鉴定出piRNA.与人和果蝇等少数模式生物相比，昆虫的piRNA研究较为滞后，迄今仅在果蝇[12]、白纹伊蚊(Aedesalbopictus)[13]和蜜蜂[5]等少数昆虫中有过相关报道.Jones等人[14]研究表明piRNA途径可能在维持果蝇基因组完整性、组织功能、代谢稳态和免疫防御方面起到一定作用.Guo等通过过表达和敲减piRNA-3312发现piRNA-3312可靶向结合肠道酯酶1基因进而负调控淡色库蚊Culexpipienspallens对拟除虫菊酯的抗性[15].

目前，关于piRNA参与蜜蜂宿主与病原互作的研究未见报道；东方蜜蜂微孢子虫胁迫下意蜂工蜂的piRNA应答研究仍然缺失.本研究拟基于已获得的small RNA(sRNA)组学数据对东方蜜蜂微孢子虫胁迫下意蜂工蜂中肠的DEpiRNA进行筛选、分析和验证，进而对DEpiRNA进行靶向预测及分析，以期解析意蜂工蜂中肠响应东方蜜蜂微孢子虫胁迫的piRNA差异表达谱及调控网络，并揭示DEpiRNA调控宿主应答的潜在作用.

2 材料与方法

2.1 材 料

2.1.1 供试蜜蜂及微孢子虫 意蜂工蜂取自福建农林大学动物科学学院(蜂学学院)教学蜂场的实验蜂群.东方蜜蜂微孢子虫孢子由福建农林大学动物科学学院(蜂学学院)蜜蜂保护实验室纯化和保存[16].

2.1.2 sRNA组学数据来源 前期研究中，笔者所在团队已利用small RNA-seq(sRNA-seq)技术对东方蜜蜂微孢子虫胁迫7 d和10 d的意蜂工蜂中肠(AmT1和AmT2)及未受胁迫的工蜂中肠(AmCK1和AmCK2)进行测序，共测得16 589 574条原始读段(raw reads)，过滤后得到132 028 990条有效读段(clean reads)，各组中clean reads占比均达到77.04%及以上；AmCK1、AmT1、AmCK2和AmT2组内3个不同生物学重复之间的Pearson相关性系数分别达到99.93%、87.92%、99.67%和91.14%[17].高质量的sRNA-seq组学数据可为本研究中piRNA的差异表达及调控网络分析提供可靠的数据基础.

2.2 方 法

2.2.1 DEpiRNA的筛选与分析 根据edgeR软件(http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)对AmCK1 vs AmT1和AmCK2 vs AmT2比较组中的DEpiRNA进行筛选，筛选标准为|log2(Fold change)|≥1(|log2FC|≥1)且P≤0.05，其余采用软件默认参数.利用基迪奥在线云平台(www.omicshare.com)的相关工具对上述两个比较组中的DEpiRNA进行Venn分析，采用默认参数.

2.2.2 DEpiRNA的靶mRNA预测与分析 参照杜宇等[18]的方法，利用TargetFinder软件[19]对DEpiRNA进行靶向预测.通过BLAST工具将靶mRNA比对到GO数据库(https://www.geneontology.org)和KEGG数据库(https://www.genome.jp/kegg/)以获得相应的功能和通路注释.最后利用Cytoscape软件[20]分别对上述两个比较组中的DEpiRNA及其靶mRNA、共有DEpiRNA及其靶mRNA的调控网络进行可视化.均采用默认参数.

2.2.3 DEpiRNA的Stem-loop RT-PCR与qPCR验证 利用FastPure®Cell/Tissue Total RNA Isolation Kit V2(Vazyme公司，中国)分别提取意蜂7日龄和10日龄工蜂中肠样品的总RNA.利用超微量分光光度仪(Thermo Fisher，美国)测定RNA质量.从上述两个比较组中随机选取4个DEpiRNA(piR-ame-967173、piR-ame-459976、piR-ame-779990、piR-ame-1044429)进行验证.根据上述piRNA的核酸序列设计相应的特异性Stem-loop引物和上游引物(F)以及通用下游引物(R)，委托上海生物工程有限公司合成引物(附表1).参照HiScript®1st Strand cDNA Synthesis Kit试剂盒说明书利用Stem-loop引物反转录合成cDNA作为piRNA的PCR模板；利用随机引物和oligo (dT)引物进行常规反转录，合成的cDNA作为内参基因snRNA U6的PCR模板.PCR体系(20 μL)包括：10 μL PCR mix(南京诺唯赞，中国)，1 μL上游引物(F)、1 μL通用下游引物(R)、1 μL cDNA模板以及7 μL DEPC处理水.PCR反应条件为：95 ℃ 5 min，95 ℃ 10 s，55 ℃ 30 s，40个循环.产物经1.8%的琼脂糖凝胶电泳检测.采用凝胶图像成像分析系统(上海培清，中国)进行凝胶观察和拍照.

参照笔者所在团队前期已建立的方法[18]对上述4个DEpiRNA进行qPCR验证.利用Stem-loop引物进行反转录得到相应的cDNA，作为模板进行qPCR反应.反应按照SYBR Green Dye试剂盒(Vazyme公司，中国)操作说明书进行.反应体系(20 μL)包含：SYBR Green Dye 10 μL，cDNA模板1.3 μL，正、反向引物各1 μL，DEPC水6.7 μL.反应在ABI QuantStudio 3荧光定量PCR仪(ABI，美国)上进行，程序设置为：95 ℃预变性5 min；95 ℃变性10 s，60 ℃退火延伸30 s，共40个循环.每个反应进行3次技术重复.采用2-ΔΔCt法计算DEpiRNA相对表达量.最后利用Graph Prism 7软件对qPCR数据进行student’st检验，并计算组间差异显著性及绘图，采用默认参数.

3 结果与分析

3.1 东方蜜蜂微孢子虫胁迫下意蜂工蜂中肠的piRNA差异表达谱

差异表达分析结果显示，AmCK1 vs AmT1比较组包含50个DEpiRNA，上调piRNA和下调piRNA的数量分别为46和4个；其中上调倍数最高的是piR-ame-1207641(log2FC=16.81,P=6.37E-9)，其次是piR-ame-500944(log2FC=13.13,P=4.45E-6)和piR-ame-748814(log2FC=12.49,P=5.14E-7)；下调倍数最高的是piR-ame-1221425(log2FC=-19.97,P=0.0001)，其次是piR-ame-1182394(log2FC=-19.65,P=0.0001)和piR-ame-1219424(log2FC=-17.26,P=0.0002).AmCK2 vs AmT2比较组包含207个DEpiRNA，上调piRNA和下调piRNA的数量分别为56和151个；其中上调倍数最高的是piR-ame-1008436(log2FC=15.01,P=3.11E-7)，其次是piR-ame-771408(log2FC=15.01,P=3.11E-7)和piR-ame-716466(log2FC=14.95,P=5.36E-7)；下调倍数最高的是piR-ame-706008(log2FC=-14.46,P=7.49E-6)，其次是piR-ame-634882(log2FC=-14.09,P=9.88E-7)和piR-ame-1253391(log2FC=-13.96,P=3.08E-5).Venn分析结果显示，上述两个比较组的共有DEpiRNA为10个(表1)，特有DEpiRNA分别为40和197个.

表1 东方蜜蜂微孢子虫胁迫7 d和10 d的意蜂工蜂中肠共有DEpiRNA的详细信息

3.2 东方蜜蜂微孢子虫胁迫下意蜂工蜂中肠DEpiRNA的靶mRNA预测及注释

靶向预测结果显示，AmCK1 vs AmT1比较组中piR-ame-748814、piR-ame-762269、piR-ame-1128833、piR-ame-175077可分别靶向2035、927、3021和1273条mRNA；AmCK2 vs AmT2比较组中piR-ame-1048676、piR-ame-1055898、piR-ame-11093、piR-ame-1128833、piR-ame-165683、piR-ame-190949、piR-ame-202265、piR-ame-471566、piR-ame-503081、piR-ame-630674、piR-ame-706008、piR-ame-742536、piR-ame-797627、piR-ame-932156可分别靶向1290、631、1195、3021、3549、4040、1152、2583、512、4935、1297、3191、593、1163条mRNA.

GO数据库注释结果显示，AmCK1 vs AmT1比较组中DEpiRNA的靶mRNA可注释到细胞进程(1381)和代谢进程(1214)等18个生物学进程相关条目，细胞器(503)和细胞(687)等16个细胞组分相关条目，结合(1509)和转运活性(161)等11个分子功能相关条目(图1a)；AmCK2 vs AmT2比较组中DEpiRNA的靶mRNA可注释到细胞进程(3744)和生物学调控(1087)等19个生物学进程相关条目，细胞膜组件(1529)和细胞膜(1581)等17个细胞组分相关条目，信号转导活性(417)和结构分子活性(76)等11个分子功能相关条目(图1b).括号内的数字代表注释在该条目的靶mRNA数目.

KEGG数据库注释结果显示，AmCK1 vs AmT1比较组中DEpiRNA的靶mRNA共注释到117条代谢通路，包括Jak-STAT(6)、磷脂酰肌醇信号系统(46)、Wnt信号通路(42)、mTOR信号通路(9)、氨基酸代谢(10)、碳水化合物代谢(14)及脂质代谢(10)等(图2a)；AmCK2 vs AmT2比较组中DEpiRNA的靶mRNA共注释到140条代谢通路，包括Wnt信号通路(117)、FoxO信号通路(82)、MAPK信号通路(26)、氨基酸代谢(13)、碳水化合物代谢(14)及辅助因子和维生素代谢(12)等(图2b).括号内的数字代表注释在该通路的靶mRNA数目.

3.3 东方蜜蜂微孢子虫胁迫下意蜂工蜂中肠DEpiRNA及其靶mRNA的调控网络

调控网络分析结果显示，AmCK1 vs AmT1比较组中同一个DEpiRNA可靶向多个mRNA，如piR-ame-748814可靶向多达2035个mRNA；此外，部分mRNA可被1～2个DEpiRNA靶向结合，如NM_001014992.1可被piR-ame-1128833和piR-ame-175077共同靶向结合，NM_001040233.1可被piR-ame-762269和piR-ame-175077共同靶向结合(图3a).类似地，AmCK2 vs AmT2比较组中也存在同一个DEpiRNA靶向多个mRNA、部分mRNA被多个DEpiRNA靶向结合的现象(图3b).

图1 AmCK1 vs AmT1 (a)和AmCK2 vs AmT2 (b)比较组DEpiRNA的靶mRNA的GO数据库注释Fig.1 GO database annotation of mRNAs targeted by DEpiRNAs in AmCK1 vs AmT1 (a) and AmCK2 vs AmT2 (b) comparison groups

图2 AmCK1 vs AmT1 (a)和AmCK2 vs AmT2 (b)比较组DEpiRNA的靶mRNA的KEGG数据库注释

图3 东方蜜蜂微孢子虫胁迫7 d和10 d的意蜂工蜂中肠中DEpiRNA及其靶mRNA的调控网络Fig.3 Regulation network of DEpiRNAs and their target mRNAs in A. m. ligustica workers’ midguts at 7 d and 10 d post N.ceranae stress

进一步分析发现，上述两个比较组共有的DEpiRNA piR-ame-1128833靶向的mRNA中有30个涉及Hippo、Notch、FoxO、Wnt和Hedgehog等发育相关的信号通路，有54个涉及内吞作用、MAPK信号通路、Jak-STAT信号通路及泛素介导的蛋白水解等免疫途径.

3.4 DEpiRNA的Stem-loop RT-PCR与RT-qPCR验证

利用Stem-loop RT-PCR对随机选取4个DEpiRNA进行扩增，产物的电泳结果显示均能扩增出符合预期大小(约80 bp)的目的片段(图4).进一步的qPCR结果显示上述4个DEpiRNA的表达趋势与测序数据中的趋势一致.以上结果证实了本研究中piRNA表达趋势的真实性和测序数据的可靠性.

4 讨 论

本研究选择东方蜜蜂微孢子虫胁迫7 d和10 d两个时间点的意蜂工蜂中肠主要基于三点依据：(1) 东方蜜蜂微孢子虫孢子接种西方蜜蜂工蜂后，病原孢子载量在宿主细胞可持续增长到 20 d[21]；(2) 东方蜜蜂微孢子虫孢子接种意蜂工蜂后，接种组工蜂的累积死亡率在1～10 dpi阶段均高于未接种组工蜂，但两组累计死亡率只在7 dpi和10 dpi两个时间点具有显著性差异[22]；(3) 笔者团队前期

图4 DEpiRNA的Stem-loop RT-PCR(a)和qPCR验证(b～e)Fig.4 Stem-loop RT-PCR (a) and qPCR(b～e) validation of DEpiRNAs

已基于东方蜜蜂微孢子虫接种7 d和10 d的意蜂工蜂中肠及未接种工蜂中肠转录组数据解析了宿主响应胁迫的免疫基因应答[23]，分析和探讨了miRNA和lncRNA在宿主的胁迫应答中的调控作用[17,22]，本研究是在前期基础上进一步探讨东方蜜蜂微孢子虫胁迫下意蜂工蜂的piRNA差异表达谱及潜在功能，以期全面深入地阐释意蜂工蜂响应东方蜜蜂微孢子虫胁迫的作用机制.相比于未受胁迫的意蜂工蜂中肠，我们在东方蜜蜂微孢子虫胁迫7 d和10 d的工蜂中肠中分别鉴定到50和207个DEpiRNA，说明这些DEpiRNA参与了宿主的胁迫应答.AmCK1 vs AmT1和AmCK2 vs AmT2比较组共有的DEpiRNA为10个，特有的DEpiRNA分别为40和197个，推测上述共有DEpiRNA在宿主应答过程的2个时间点均发挥作用，而多数DEpiRNA在宿主应答的不同时间点发挥特定作用.

piRNA参与了许多后生动物的生殖系发育和配子发生且与人类的癌症高度相关[24].piRNA也被证实能通过靶向mRNA发挥基因表达调控的作用[25].本研究中，AmCK1 vs AmT1中仅有4个DEpiRNA可靶向7256条mRNA，AmCK2 vs AmT2中有14个DEpiRNA可靶向29152条mRNA，说明少数DEpiRNA参与宿主应答过程中基因表达调控，而多数DEpiRNA可能主要发挥保持基因组稳定性的作用.AmCK1 vs AmT1中DEpiRNA的靶mRNA可注释到细胞进程等45个功能条目(图1a)；还注释到磷脂酰肌醇信号系统等117条代谢通路(图2a).类似地，AmCK2 vs AmT2中DEpiRNA的靶mRNA可注释到细胞进程等47个功能条目(图1b)以及MAPK信号通路等140条代谢通路(图2b).以上结果表明DEpiRNA通过调节诸多生物学过程潜在影响宿主对东方蜜蜂微孢子虫胁迫的应答.miRNA作为基因表达的关键调控因子，其“种子区”根据碱基互补配对原则靶向结合mRNA的3’UTR特定区域，同一个miRNA能够靶向多个mRNA，反之亦然[26].本研究也发现，同一个DEpiRNA可靶向结合多个mRNA，一些mRNA可同时被多个DEpiRNA靶向结合，说明DEpiRNA能以一种类似于miRNA的方式参与宿主应答中的基因表达调控[27].

信号转导途径在真核生物的环境适应、细胞活动、生长发育、新陈代谢和免疫应答等生物学过程中发挥举足轻重的作用[28].Hippo信号通路通过调节细胞增殖和凋亡参与器官大小控制和组织再生[29].Notch信号通路可调节果蝇中肠上皮细胞分化及细胞凋亡[30].FoxO信号通路参与调控碳水化合物代谢和能量代谢等重要代谢过程[31].同时，FoxO和Wnt信号通路在调节细胞生长、增殖和凋亡及新陈代谢等方面均起到关键作用[32-34].Hedgehog信号通路与细胞的生长和发育密切相关[35].上述两个比较组共有的DEpiRNA仅piR-ame-1128833可靶向3021条mRNA，暗示piR-ame-1128833在意蜂工蜂响应东方蜜蜂微孢子虫胁迫的过程中发挥重要作用.本研究发现，piR-ame-1128833的靶mRNA中分别有11、1、2、11和5个可注释到Hippo、Notch、FoxO、Wnt和Hedgehog等信号转导途径，表明piR-ame-1128833通过调控生长、发育和物质能量代谢通路潜在影响意蜂工蜂对东方蜜蜂微孢子虫胁迫的应答.

为抵御病原体和寄生虫的入侵，昆虫进化出了高效的先天免疫系统，包括以内吞和吞噬作用为代表的细胞免疫[36]和以Jak-STAT和MAPK等信号通路及抗菌肽释放为代表的体液免疫[37].媒介昆虫的中肠肠壁是病毒入侵的主要屏障，研究表明依赖于网格蛋白的内吞作用在促进病毒进入烟粉虱(Bemisiatabaci)中肠细胞方面起到重要作用[38].MAPK信号通路可以通过级联反应调控中肠基因表达进而引起小菜蛾(Plutellaxylostella)的抗病毒反应[39].JAK/STAT信号通路参与调节细胞生长、分化、凋亡及炎症反应，并与埃及伊蚊(Aedesaegypti)抵御真菌侵染密切相关[40].泛素介导的蛋白水解系统作为消除错误折叠或损伤细胞的主要机制，是许多生理过程如信号转导、细胞周期及免疫细胞功能的控制器[41].本研究发现，piR-ame-1128833的靶mRNA中分别有32、2、1和19个可注释到内吞作用、MAPK信号通路、Jak-STAT通路及泛素介导的蛋白水解系统，说明piR-ame-1128833可能通过调节上述细胞和体液免疫途径参与宿主对东方蜜蜂微孢子虫胁迫的应答，但仍需要进一步实验验证.目前，参照miRNA的过表达和敲减方法对真核生物组织或细胞中的piRNA进行功能研究已见诸报道[42,43].下一步我们将根据piR-ame-1128833的序列设计合成相应的模拟物和抑制物，并通过饲喂法对东方蜜蜂微孢子虫感染的意蜂工蜂内的piR-ame-1128833进行过表达和敲减，以明确该piRNA在宿主的胁迫应答中的功能.

综上，东方蜜蜂微孢子虫胁迫可引起意蜂工蜂中肠piRNA表达谱的改变；DEpiRNA可通过靶向mRNA潜在调控宿主响应胁迫中的多个生物学过程；piR-ame-1128833可通过靶向调控相关mRNA潜在调节Hippo和Wnt等信号通路及内吞作用和MAPK信号通路等免疫途径.