意大利蜜蜂工蜂中肠piRNA的鉴定与分析

范小雪, 隆 琦, 孙明会, 郭意龙, 赵浩东, 宋岳梅,康育欣, 顾小雨, 陈大福,2,*, 郭 睿, 2,*

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

非编码RNA(non-coding RNA, ncRNA)起初被认为是转录过程的“噪音”和副产物，不具有任何生物学功能(Kaikkonen and Adelman, 2018)。但随着组学和分子生物学技术的持续进步和相关研究的不断深入，ncRNA被证实在动植物和微生物的生长发育和新陈代谢等过程中发挥重要的调控作用(Quinn and Chang, 2016)。ncRNA根据长度可分为长度较长的长链非编码RNA(long non-coding RNA, lncRNA)和长度较短的小ncRNA，后者可进一步分为小干扰RNA(small interfering RNA, siRNA)、微小RNA(microRNA, miRNA)和Piwi蛋白互作RNA(PIWI-interacting RNA, piRNA)等。笔者所在团队前期联用基于链特异性cDNA建库的RNA-seq与small RNA-seq(sRNA-seq)技术对意大利蜜蜂Apismelliferaligustica7日龄和10日龄工蜂中肠进行测序，并系统解析了意大利蜜蜂工蜂中肠发育过程的ncRNA差异表达谱和调控网络以及差异表达ncRNA的调控作用(郭睿等, 2018a, 2018b; 熊翠玲等, 2018; 杜宇等, 2020a, 2020b)。

piRNA较miRNA(长度约为22～24 nt)和siRNA(长度约为20～25 nt)稍长，且3′末端含有2′-O-甲基修饰位点，可与Argonaute蛋白的亚族蛋白之一的PIWI蛋白家族相互作用形成piRNA介导的沉默复合体(piRNA-induced silencing complex, piRISC)。较多的研究表明piRNA在生殖干细胞分化、胚胎发育、DNA完整性维持、表观遗传学调控及物种的性别决定等方面发挥重要作用(Daietal., 2020)。近十多年来，piRNA在黑腹果蝇Drosophilamelanogaster(Aravinetal., 2001)、人Homosapiens(Girardetal., 2006)、小鼠Musmusculus(Aravinetal., 2006)等昆虫、哺乳动物中被陆续鉴定出来。虽然piRNA的首次发现是在黑腹果蝇的生殖腺中(Aravinetal., 2001)，但与脊椎动物相比，昆虫的piRNA研究十分滞后。piRNA能以一种类似于miRNA的方式通过与靶mRNA的碱基互补配对调控靶基因表达，Kotelnikov等(2009)研究发现piRNA能与Aub结合形成piRNA-Aub，进而在转录后水平特异性降解黑腹果蝇Stellate的mRNA，从而实现对Stellate的沉默。

近期的研究结果表明昆虫piRNA与抗病毒免疫过程密切相关(Kolliopoulouetal., 2019)，蜜蜂的内源性小RNA(small RNA, sRNA)研究迄今仍比较有限且主要集中在miRNA(Cristinoetal., 2014)，piRNA相关研究尤为匮乏。目前，piRNA调控意大利蜜蜂工蜂中肠发育的研究缺失，中肠发育的分子调控机制不明。本研究基于已获得的sRNA组学数据对piRNA进行全转录组范围的鉴定和分析，进而筛选出中肠发育过程的差异表达piRNA(differentially expressed piRNA, DEpiRNA)并进行靶向预测及分析，以期丰富西方蜜蜂的piRNA信息，并为探明piRNA调控意大利蜜蜂工蜂中肠发育的分子机理提供基础。

1 材料与方法

1.1 试虫来源

实验用意大利蜜蜂工蜂取自福建农林大学动物科学学院(蜂学学院)教学蜂场。

1.2 sRNA组学数据来源

前期研究中，笔者所在团队已利用sRNA-seq技术对意大利蜜蜂7日龄(Am7)和10日龄(Am10)工蜂中肠进行测序，获得了高质量的sRNA组学数据，质控结果显示经sRNA-seq后Am7和Am10分别测得13 511 614条和13 448 175条原始读段(raw reads)，过滤后分别得到13 119 002条和13 010 751条有效读段(clean reads)，两组内生物学重复间的Pearson相关性分别达到0.9870和0.9988以上(杜宇等, 2020a)。原始数据已上传NCBI SRA数据库，获得BioProject号为PRJNA408312。高质量的sRNA组学数据可为本研究中piRNA的鉴定和分析提供数据基础。

1.3 piRNA的生物信息学预测及分析

首先将上述clean reads比对西方蜜蜂参考基因组(Assembly Amel_4.5)，获得能比对上参考基因组的clean reads即mapped reads，然后按照以下步骤过滤其他类型的sRNA并预测piRNA：(1)通过比对数据库滤除rRNA, scRNA, snoRNA, snRNA和tRNA等类型的非编码小RNA；(2)进一步滤除剩余clean reads中的miRNA；(3)根据piRNA的长度特征，保留长度介于24～33 nt的sRNA；考虑到部分sRNA比对基因组时会出现同时比对到多个不同位置的情况，因此进一步滤除比对到多个位置的sRNA，仅保留比对到唯一位置的sRNA作为候选piRNA。 根据TPM(tags per million)法(TPM=T×106/N，T表示piRNA的clean reads，N表示总sRNA的clean reads)对Am7和Am10组中piRNA的表达量进行归一化处理。采用GraphPad Prism 7软件分析各组中piRNA的结构特征并绘图。

1.4 DEpiRNA的筛选及靶向分析

按照|log2fold change|≥1且P≤0.05的标准，筛选Am7vsAm10比较组的DEpiRNA。使用TargetFinder软件(Allenetal., 2005)预测DEpiRNA的靶mRNA，采用软件默认参数。利用BLAST工具将DEpiRNA的靶mRNA比对到GO数据库(https:∥www.geneontology.org)和KEGG数据库(https:∥www.genome.jp/kegg/)以获得功能和通路注释，根据上述DEpiRNA与mRNA的靶向结合关系，通过Cytoscape软件(Smootetal., 2011)进行调控网络的可视化。

1.5 piRNA的Stem-loop RT-PCR验证

利用FastPure®Cell/Tissue Total RNA Isolation Kit V2(Vazyme公司，南京)分别提取意大利蜜蜂7日龄和10日龄工蜂中肠样品的总RNA，3只中肠作为一个生物学重复，3次生物学重复，3次技术重复。利用超微量分光光度仪(Thermo Fisher，美国)测定RNA质量。随机选取Am7和Am10两组共有的6个piRNA(piR-ame-1084826, piR-ame-11093, piR-ame-14476, piR-ame-24995, piR-ame-39500和piR-ame-774987)进行验证。根据上述piRNA的核酸序列设计相应的特异性Stem-loop引物和上游引物(F)以及通用下游引物(R)，委托上海生物工程有限公司合成引物(表1)。参照HiScript®1st Strand cDNA Synthesis Kit试剂盒说明书(Vazyme公司，南京)利用Stem-loop引物反转录合成cDNA作为piRNA的PCR模板，利用随机引物和oligo (dT)引物进行常规反转录，合成的cDNA作为内参基因U6 snRNA(GenBank登录号: 725641)的PCR模板。PCR反应体系(20 μL): PCR Mix(Vazyme公司，南京)10 μL, 上游引物(2.5 pmol/μL) (F) 1 μL,通用下游引物(2.5 pmol/μL) (R) 1 μL, cDNA模板1 μL, DEPC处理水7 μL。PCR反应程序: 95℃ 5 min; 95℃ 10 s, 55℃ 30 s, 40个循环。产物经1.8%的琼脂糖凝胶电泳检测。采用凝胶图像成像分析系统(上海培清)进行凝胶观察和拍照。

1.6 DEpiRNA的RT-qPCR验证

参照笔者所在团队前期已建立的方法(杜宇等, 2020a)，随机挑选6个DEpiRNA(piR-ame-1084826, piR-ame-11093, piR-ame-14476, piR-ame-24995, piR-ame-39500和piR-ame-774987)进行RT-qPCR验证，U6作为内参基因(GenBank登录号: 725641)。利用反转录试剂盒(Vazyme公司，南京)对1.5节中7日龄和10日龄工蜂中肠样品的总RNA进行Stem-loop反转录得到cDNA，作为模板进行RT-qPCR反应。反应按照SYBR Green Dye试剂盒(Vazyme公司，南京)操作说明书进行。反应体系(20 μL): SYBR Green Dye 10 μL, cDNA模板1.3 μL, 正反向引物(2.5 pmol/μL)各1 μL, DEPC水6.7 μL。反应在ABI QuantStudio 3荧光定量PCR仪(ABI，美国)上进行，反应程序: 95℃ 5 min; 95℃ 10 s, 60℃ 30 s, 共40个循环。3只中肠作为一个生物学重复，3次生物学重复，每个反应均进行3个平行重复和3次技术重复，引物信息见表1。采用2-ΔΔCt法计算piRNA相对表达量。最后利用Graph Prism 7软件对qPCR数据进行Student氏t检验、计算组间差异显著性及绘图。

表1 本研究所用引物Table 1 Primers used in this study

2 结果

2.1 意大利蜜蜂工蜂中肠piRNA的数量与结构特征

意大利蜜蜂7日龄(Am7)和10日龄(Am10)工蜂中肠的sRNA组学数据中滤除rRNA, scRNA, snoRNA, snRNA, tRNA和miRNA后的piRNA长度分布介于24～33 nt，Am7组中分布在24～30 nt范围的piRNA数量高于Am10组的，但分布在31～33 nt的piRNA数量低于Am10组的(图1)；其中在Am7组和Am10组中共鉴定到596个piRNA。

Am7中长度介于24～29 nt的piRNA的首位碱基主要偏向U，长度介于30～32 nt的piRNA首位碱基主要偏向于G，而长度为33 nt的piRNA首位碱基主要偏向A(图2: A)；Am10中长度介于24～25 nt的piRNA的首位碱基主要偏向A，长度介于26～29 nt的piRNA首位碱基主要偏向C，而长度介于30～32 nt的piRNA的首位碱基主要偏向G(图2: B)。

图2 意大利蜜蜂Am7(A)和Am10(B)中肠中piRNA的首位碱基偏向性Fig. 2 Bias of the first base of piRNAs in the midgut of Am7 (A) and Am10 (B) of Apis mellifera ligustica

2.2 意大利蜜蜂工蜂中肠发育过程DEpiRNA的靶mRNA

在Am7vsAm10比较组中共筛选出41个DEpiRNA，包括1个上调piRNA(piR-ame-1233597, log2fold change=1.767,P=4.38E-03)和40个下调piRNA，其中下调幅度最大的piRNA为piR-ame-1008436(log2fold change=-15.234,P=1.57E-05)，其次为piR-ame-24995(log2fold change=-14.99,P=8.57E-06)和piR-ame-638412(log2fold change=-14.86,P=1.65E-04)。

靶向预测结果显示仅piR-ame-11093(log2fold change=-13.91,P=9.23E-06), piR-ame-1111451(log2fold change=-13.15,P=9.80E-07), piR-ame-190949(log2fold change=-1.63,P=0.0069)和piR-ame-932156(log2fold change=-1.54,P=0.0089)可分别靶向1 195, 1 018, 4 040和1 063条mRNA，其余piRNA未预测到靶mRNA。

GO数据库注释结果显示上述4个piRNA的靶mRNA可注释到细胞进程(1 440条)、代谢进程(1 235条)、单一组织进程(1 106条)和生物学调控(421条)等18个生物学过程相关功能条目，细胞(631条)、细胞组件(631条)和细胞膜(626条)等17个细胞组分相关条目，结合(1 482条)、催化活性(1 069条)和转运活性(184条)等10个分子功能相关条目，共计45个功能条目(图3)。

图3 意大利蜜蜂Am7 vs Am10比较组中肠中差异表达piRNA(DEpiRNA)的靶mRNA的GO数据库注释Fig. 3 GO database annotation of mRNAs targeted by differentially expressed piRNAs (DEpiRNAs) in the midgutof the Am7 vs Am10 comparison group of Apis mellifera ligustica

此外，KEGG数据库注释结果显示这些靶mRNA可注释到嘌呤代谢(74条)、内吞作用(55条)、吞噬体(49条)及Hippo信号通路(59条)等45条通路(图4)。

图4 意大利蜜蜂Am7 vs Am10比较组中肠中差异表达piRNA(DEpiRNA)的靶mRNA的KEGG数据库注释Fig. 4 KEGG database annotation of mRNAs targeted by differentially expressed piRNAs (DEpiRNAs) in the midgutof the Am7 vs Am10 comparison group of Apis mellifera ligustica

进一步分析发现，piR-ame-190949的25条靶mRNA涉及发育相关的Hippo, Hedgehog, FoxO, Wnt和Notch信号通路，免疫相关的泛素介导的蛋白水解(ubiquitin-mediated proteolysis)及MAPK, Jak-STAT和Toll-like受体信号通路(图5)；piR-ame-11093的13条靶mRNA涉及Hippo, Hedgehog和Wnt信号通路及泛素介导的蛋白水解、Toll和Imd信号通路(图5)；piR-ame-1111451的16条靶mRNA涉及Hippo和Hedgehog信号通路及泛素介导的蛋白水解(图5)；piR-ame-932156的5条靶mRNA涉及Hippo和Wnt信号通路及泛素介导的蛋白水解(图5)。

图5 意大利蜜蜂Am7 vs Am10比较组中肠中差异表达piRNA(DEpiRNA)及其靶mRNA调控网络Fig. 5 Regulatory networks between differentially expressed piRNAs (DEpiRNAs) in the midgut of the Am7 vs Am10comparison group of Apis mellifera ligustica and its target mRNAA: 涉及发育相关信号通路的DEpiRNA-mRNA调控网络 DEpiRNA-mRNA regulatory network in development-related signaling pathways; B: 涉及免疫相关信号通路的DEpiRNA-mRNA关系网络 DEpiRNA-mRNA regulatory network in immune-related signaling pathways. 蓝色六边形代表信号通路，绿色菱形代表DEpiRNA，红色矩形代表靶mRNA，六边形、菱形和矩形的大小分别代表连接信号通路、DEpiRNA和mRNA的数量。Blue hexagons indicate signaling pathways, green rhombuses indicate DEpiRNA, and red rectangles indicate target mRNA. The sizes of hexagons, rhombuses, and rectangles represent the numbers of signaling pathway, DEpiRNAs and mRNAs, respectively.

2.3 意大利蜜蜂工蜂中肠piRNA的Stem-loop RT-PCR验证

利用Stem-loop RT-PCR对随机选择的Am7与Am10两组共有的6个piRNA进行扩增，扩增产物的琼脂糖凝胶电泳结果显示在7日龄和10日龄意大利蜜蜂工蜂中肠中均扩增出符合预期大小(约28 bp)的目的片段(图6)，证明了piRNA表达的真实性。

图6 随机选择的意大利蜜蜂Am7和Am10中肠中共有的6个piRNA表达的Stem-loop RT-PCR验证Fig. 6 Stem-loop RT-PCR validation of expression of six randomly selected piRNAs sharedin the midgut of Am7 and Am10 of Apis mellifera ligustica

2.4 意大利蜜蜂工蜂中肠DEpiRNA的RT-qPCR验证

RT-qPCR结果显示7个DEpiRNA中的6个的表达趋势与测序数据中的变化趋势一致(图7)，证实了本研究中sRNA-seq数据的可靠性和piRNA差异表达谱趋势的真实性。

图7 意大利蜜蜂Am7 vs Am10比较组中肠中差异表达piRNA(DEpiRNA)的RT-qPCR验证Fig. 7 RT-qPCR verification of differentially expressed piRNAs (DEpiRNAs) in the midgutof the Am7 vs Am10 comparison group of Apis mellifera ligustica图中数据为平均值±标准误；柱上星号表示两组之间表达量的差异显著性 (*P<0.05; **P<0.01)(Student氏t检验)。Data in the figure are mean±SE. Asterisks above bars indicate significance of difference in the expression level between the two groups (*P<0.05; **P<0.01)(Student’s t-test).

3 讨论

Liao等(2010)在西方蜜蜂中成功鉴定并克隆了piRNA生物合成通路的两个关键基因Am-aub和Am-ago3，首次证实西方蜜蜂中真实存在piRNA合成机制。Wang等(2017)曾通过预测发现西方蜜蜂雄蜂体内存在大量长度介于26～31 nt的piRNA，并发现雄蜂体内piRNA的表达丰度和多样性明显高于雌性的蜂王和工蜂的。但蜜蜂的piRNA研究迄今仍非常有限。本研究基于已获得的sRNA-seq数据，通过生物信息学方法从意大利蜜蜂7日龄和10日龄工蜂中肠中共鉴定到596个piRNA。鉴于目前蜜蜂的piRNA信息匮乏，本研究鉴定到的piRNA很好地丰富了西方蜜蜂的piRNA信息，为进一步开展相关功能研究提供了数据基础。推测上述piRNA在意大利蜜蜂工蜂中肠发育过程发挥调控作用。

piRNA在物种间具有高度保守性，主要体现在长度分布和碱基偏向性(Ramat and Simonelig, 2020)。piRNA的长度分布主要介于21～31 nt，例如人类胃癌细胞中鉴定到的piRNA长度介于25～31 nt之间(夏言, 2019)；黑腹果蝇雌性性腺中鉴定到的piRNA长度介于23～30 nt之间(Huangetal., 2017)；家蚕Bombyxmori中鉴定到的piRNA长度介于25～30 nt，其中分布在28～29 nt的piRNA数量最多(刘晓, 2018)。本研究中，意大利蜜蜂工蜂中肠piRNA的长度介于24～33 nt，其中分布在24～28 nt长度范围内的piRNA数量较多，这与上述物种的piRNA的长度分布规律相似，说明不同物种的piRNA的长度分布具有较高的保守性。迄今为止，在能够产生piRNA的动物中均发现单链簇通过传统的单向转录生成piRNA前体是piRNA的主要生成方式，这种在物种间高度稳定的生成方式使得piRNA序列特征高度保守，5′末端的首位碱基具有尿嘧啶(U)偏向性。这一偏向性已在小鼠(Aravinetal., 2005)、家蚕(刘晓, 2018)和果蝇(Huangetal., 2017)等动物中被证实。本研究Am7组中长度介于24～29 nt的piRNA的首位碱基主要偏向U(图2: A)，为piRNA的首位碱基偏向性提供了又一例证。

piRNA不仅在生殖细胞中起到靶向抑制转座子进而维持基因组稳定性的作用，还能够以一种类似miRNA的方式通过碱基互补配对抑制靶基因表达，进而影响诸多重要生物学过程(Daietal., 2020)。在哺乳动物中，piRNA以类似miRNA和siRNA的作用方式介导体细胞质中mRNA的降解(Zhangetal., 2018)。意大利蜜蜂工蜂中肠是物质消化和营养吸收的主要场所，发生着活跃的物质和能量代谢反应。本研究中，在Am7vsAm10比较组中共筛选到41个DEpiRNA，但其中只有piR-ame-11093(log2fold change=-13.91), piR-ame-1111451(log2fold change=-13.15), piR-ame-190949(log2fold change=-1.63)和piR-ame-932156(log2fold change=-1.54)能够潜在靶向结合mRNA，说明多数DEpiRNA的主要功能可能是维持基因组稳定性，少数DEpiRNA参与基因表达调控。前期研究中，笔者所在团队在Am7vsAm10比较组筛选到112个DEmiRNA，其中显著上调和下调的miRNA分别靶向7 434和9 559条mRNA(杜宇等, 2020a)。本研究发现，上述4个DEpiRNA共靶向结合7 316条mRNA，可注释到119个新陈代谢相关通路，包括嘌呤代谢(74条)、丙酮酸盐代谢(27条)、脂肪酸代谢(26条)、粘蛋白型O-聚糖的生物合成(12条)、不饱和脂肪酸生物合成通路(13条)及脂肪酸合成(11条)等7个物质代谢通路(图4)，说明在意大利蜜蜂工蜂中肠发育过程中上述4个DEpiRNA通过潜在调控靶基因表达以调节物质代谢。昆虫体内存在的复杂且广泛的信号通路参与调控诸多生命活动，其中Wnt, Hippo, FoxO, TGF-beta, mTOR和Notch信号通路均已被证实与昆虫的生长发育与器官形成等密切相关(Ahmad, 2017)。Jak-STAT、Toll-like受体、Imd、JNK和MAPK信号通路和泛素介导的蛋白水解可参与昆虫的免疫应答(Brutscheretal., 2015)。本研究发现，上述4个DEpiRNA的靶mRNA可注释到发育相关的Hippo, Hedgehog, Wnt, FoxO和Notch信号通路以及免疫相关的泛素介导的蛋白水解、MAPK、Jak-STAT、Toll和Imd及Toll-like受体信号通路(图5)。这说明上述4个DEpiRNA通过调节以上发育和免疫相关通路在意大利蜜蜂工蜂中肠的发育过程中发挥作用。

非编码小RNA在调控昆虫发育方面扮演着重要角色(李灏淼等, 2020; 张卫星, 2020; 祝智威等, 2021)。目前，非编码小RNA调控蜜蜂发育的研究总体有限且集中在miRNA方面(张卫星, 2020; 祝智威等, 2021)，早期的piRNA研究多集中在生殖细胞方面(Ozataetal., 2019)，随后的研究结果表明piRNA在体细胞发育方面也具有关键调控作用。本研究首次从意大利蜜蜂工蜂中肠中鉴定到piRNA，为进一步深入开展piRNA调控意大利蜜蜂和其他蜂种肠道发育的功能和机制研究提供了理论和实验依据。作为基因表达调控的又一因子，piRNA或许在蜜蜂的级型分化、职能分工、清理行为和免疫防御等重要生物学过程发挥重要作用。

越来越多的研究表明不同类型的ncRNA可通过竞争性内源RNA(competing endogenous RNA, ceRNA)调控网络而相互作用，共同调控基因表达和生物学过程(Songetal., 2017), lncRNA和circRNA均可以通过miRNA反应元件(miRNA response element, MRE)竞争性结合mRNA，从而减弱miRNA对靶基因表达的抑制(Wangetal., 2019)。前期研究中，我们在意大利蜜蜂工蜂中肠Am7vsAm10比较组中筛选出112条DElncRNA和201个DEcircRNA可分别通过ceRNA机制靶向112个DEmiRNA潜在调节意大利蜜蜂工蜂中肠的发育过程(杜宇等, 2020a)。Betting等(2021)研究发现propiR1在白纹伊蚊Aedesalbopictus的胚胎发育早期已开始表达，并通过与piRNA 3′端种子区的碱基互补配对介导lncRNA(lnc027353)的降解，进而影响白纹伊蚊的体细胞的成熟。虽然关于piRNA调控基因表达的研究还比较有限，但鉴于piRNA表现出与miRNA类似的长度分布特征及哺乳动物piRNA靶向结合mRNA的报道，lncRNA和circRNA能够通过竞争性结合piRNA调控下游靶基因的表达。lncRNA/circRNA-piRNA-mRNA轴是否参与调控意大利蜜蜂工蜂中肠发育过程以及背后的分子机制是非常有趣且值得深入探究的科学问题。

综上，本研究从意大利蜜蜂工蜂中肠中鉴定到596个piRNA，其中的piR-ame-11093, piR-ame-1111451, piR-ame-190949和piR-ame-932156 4个DEpiRNA共靶向结合7 316条mRNA，并潜在通过调控Hippo等发育相关信号通路以及Jak-STAT等免疫通路参与意大利蜜蜂工蜂中肠的发育过程。