赫氏颗石藻(Emiliania huxleyi)响应病毒感染的microRNA转录组分析

张恩权，蔡伟聪，李桂玲，李健，刘静雯

研究报告

赫氏颗石藻()响应病毒感染的microRNA转录组分析

张恩权，蔡伟聪，李桂玲，李健，刘静雯

集美大学食品与生物工程学院，厦门 361000

海洋颗石藻病毒–宿主互作是影响海洋碳、硫生物地化循环及全球气候变化的重要环节。作为大的双链DNA病毒，颗石藻病毒进化出一种“病毒细胞代谢”模式，通过重编程宿主代谢途径以满足其代谢需求，但对这一代谢模式的调控机制尚缺乏足够的认识。MicroRNA (miRNA)作为一种基因表达调控的重要因子，能够通过调控代谢过程中的靶基因表达，从而调节相关代谢通路。本研究采用small RNA测序技术分析病毒感染颗石藻差异表达的miRNA及其靶基因功能，鉴定出26条成熟miRNA (包括2条病毒来源的miRNA)，均来自23条新的miRNA前体序列，其中5条miRNA显著差异表达，包括4条上调，1条下调。实时荧光定量PCR (quantitative real-time PCR, qRT-PCR)验证结果与miRNA-seq结果基本一致。功能富集分析显示，5个差异表达的miRNA可能参与调节糖代谢、脂代谢和氨基酸等代谢。此外，差异表达miRNA的表达水平与脂质代谢相关靶基因如、、、、、等的表达水平呈负相关，说明这些miRNA可能在病毒感染过程中对宿主的脂质代谢发挥重要的调控作用。

赫氏颗石藻；颗石藻病毒；small RNA测序；microRNA；脂代谢

海洋颗石藻(Coccolithophores)是一种全球广泛分布且具有重要生态功能的真核微型浮游植物[1～3]，其中赫氏颗石藻(, Eh)具有形成“球石粒”和高产二甲基硫化物(DMSP)的能力，且几乎每年都在大洋(尤其是高纬度海域)中形成大面积赤潮[4]，该藻赤潮的迅速大规模消亡被证实是特异性病毒(virus, EhV)感染和裂解所致[5,6]。因此，EhV-Eh的互作过程是影响海洋碳、硫生物地化循环及全球气候变化的重要环节。颗石藻病毒与宿主间侵染和抵抗的博弈形成了一种新的、稳定的协同进化代谢模式，是研究真核浮游植物宿主-病毒互作的理想模式系统[7]。

脂代谢是当今病毒-宿主互作研究的重要问题。作为大的双链DNA病毒，EhVs进化出一种全新的“病毒细胞代谢(virocell metabolism)”模式，通过重编程宿主代谢途径以满足其更高的代谢需求[8,9]。Evans等[10]2009年首次报道在赤潮消亡过程中，病毒感染诱导宿主细胞脂肪酸的组成由多不饱和向单不饱和转变，导致被病毒感染的宿主细胞减少多不饱和脂肪酸向食物链中更高级营养水平的传递，从而降低海洋生态系统的总生产力。本课题组前期研究发现，病毒感染早期宿主细胞高饱和三酰基甘油(triacylglycerols, TAGs)合成积累、并形成脂滴聚集在病毒颗粒中促进病毒外壳疏水蛋白的存储，满足病毒复制和组装的需求[11]。最新研究发现，病毒感染颗石藻能够以外泌体形式富集并分泌TAGs，以此加速病毒感染过程[12]。以色列学者Carmit团队2016年首次发现，EhV编码的丝氨酸棕榈酰转移酶(viral serine palmitoyltransferase,：鞘脂从头生物合成途径中的第一个关键限速酶)在感染细胞中能够催化合成病毒特有的新型鞘糖脂(viral-specific glycosphingolipids, vGSLs)，该物质被认为是病毒组装的必要条件[9]，这表明病毒能够利用编码的辅助代谢基因(virus-encoded auxiliary metabolic genes, vAMGs)重塑宿主鞘脂代谢途径。同时vGSLs还能作为重要的信号分子诱发ROS(H2O2)[13]和NO[14]的产生，进而启动宿主细胞凋亡程序，最终裂解细胞释放病毒粒子[15～17]。然而，目前对EhV感染重塑的脂代谢的miRNA调控机制尚缺乏足够的认识。

MicroRNA(miRNA)是一种长度分布在18～26 nt的内源性非编码小RNA。动植物miRNA来自基因组位置略有不同，动物miRNA通常位于内含子区，而植物miRNA更多地来自基因组的基因间区[18]。miRNA通过结合在靶基因的不同位点抑制或激活靶基因的表达水平，作为丰度极高的基因调控因子参与调控多种生物学过程。最近在高等动植物研究中发现，miRNA能够通过调节代谢过程中的靶基因，从而调控相关代谢通路，特别是在脂质合成、脂肪酸氧化及脂蛋白分泌等代谢网络中起重要调节作用[19]。相对于高等动植物而言，由于浮游植物种内和种间的差异性特点，且目前只有少数种类基因组信息得到较完整的诠释，因此对miRNA在海洋浮游植物生物及非生物胁迫中的重要作用了解较少。目前，主要针对氮、磷、硅及铁等营养胁迫条件下miRNA的表达情况进行了初步研究[20～22]。有关浮游植物响应病原菌感染相关miRNA的研究国内外尚未报道。

本研究以Eh-EhV模式系统为研究对象，在病毒感染早期(6 h)和晚期(45 h)收集藻细胞样品，采用高通量小RNA测序技术筛选病毒诱导的差异表达miRNA，分析miRNA在病毒重塑宿主脂代谢过程中可能的调控作用，为从表观遗传学角度深入理解EhV-Eh互作的分子机制提供新的理论认识。

1 材料与方法

1.1 赫氏颗石藻的培养及病毒感染

赫氏颗石藻BOF92及其特异性裂解病毒株系EhV99B1均由挪威卑尔根大学生物系微生物研究所Gunnar Bratbak教授馈赠并保存于本实验室，藻株和病毒株均为纯化株系。颗石藻培养采用70%海水配制的f/2-si加富培养基，培养条件为：温度16℃±2℃，光照强度为60 μmol quanta m–2s–1，光照周期为14/10 (光/暗)。在2 L指数生长期的藻液中(～106cells/mL)，以1∶50 (EhV∶Eh)体积比加入浓缩病毒裂解液(病毒初始浓度约为107viruses/mL)，另外2 L添加等量高温灭活的病毒作为对照组。基于本课题组前期mRNA转录组和脂质代谢组学分析结果，EhV感染早期(6 h)和中后期(45 h)，宿主转录组发生了显著差异表达[23]，细胞脂质代谢产物的积累也发生了明显变化[11]，因此本研究选择病毒感染的6 h和45 h作为采样时间点。分别于病毒感染6 h和45 h离心收集500 mL藻细胞样品(4℃，7000 r/min，5min)，立即置于液氮中速冻，于−80℃保存备用。每个样本设置两个生物学平行，共计8个样本：Con_6 h-1、Con_6 h-2、Exp_6 h-1、Exp_6 h-2、Con_45 h-1、Con_45 h-2、Exp_45 h-1、Exp_45 h-2。

1.2 RNA的提取、small RNA文库构建及测序

用mirVana microRNA Isolation Kit (Ambion，美国)试剂盒提取样本总RNA，NanoPhotometer (IMPLEN，德国)检测总RNA的质量和纯度，Agilent 2100 BioAnalyzer系统和RNA 6000 Nano chip (Agilent，美国)分析总RNA的完整性。利用Illumina smallRNA- seq文库构建试剂盒(KAPA Biosystems，美国)制备small RNA测序文库，在Illumina HiSeqTM 2500平台进行测序(华大基因公司，深圳)。参考基因组为近缘颗石藻株系CCMP1516 (https://www.ncbi. nlm.nih.gov/genome/2?genome_assembly_id=22489)和近缘病毒株系EhV86 (https://www.ncbi.nlm.nih. gov/genome/?term=GCA_000865825.1)。

1.3 生物信息学分析

测序获得的原始raw reads过滤筛选后得到clean reads。去除接头序列、长度小于18 bp或者大于30 bp的reads、含有N碱基的reads及低质量的reads。后续所有分析基于clean reads。筛选过程基于Cutadapt 1.7.1[24]和Fastx toolkit 0.0.14[25]软件。Clean reads经过比对Rfam 11.0数据库，过滤掉rRNA、scRNA、snoRNA、snRNA以及tRNA等非编码RNA序列，过滤后的small tags作为预测miRNA的候选序列。使用BLASTN软件将剩余的small tags比对到近缘株系CCMP1516和EhV86基因组，保留能完全匹配到参考基因组的tags。使用miRDeep2[26]软件对能匹配到参考基因组的tags进行二级结构预测。miRNA的表达水平通过TPM (transcripts per million reads)方法归一化。使用edgeR软件[27]对miRNA进行差异表达分析，差异条件设置为：|log2(foldchange)| >1以及<0.05。

1.4 miRNA保守性分析

使用序列比对工具BLAST，将本实验条件下预测的24条宿主miRNAs序列与miRBase数据库收录的38,589条miRNAs、其他微藻的2650条(数据库尚未收录的miRNAs)以及在CCMP1516中已鉴定出的18条miRNA[28]进行比对分析。鉴定保守miRNA的条件为：miRNA种子序列(第2到第8个碱基)不允许有错配，且总错配数≤2。由于颗石藻基因组具有明显的种内变异，为了进一步了解BOF92的miRNAs及其前体的功能和进化情况，分析了miRNA及其前体序列与其他13株测序颗石藻株的保守性[29]。其中，92A、Eh2和Van556为深度测序的分离株，其他10株是低覆盖度测序(覆盖度可以达到91%至 95%)。对病毒来源的2条miRNA及其前体序列也与其他13株测序病毒株进行保守性分析。

1.5 miRNA靶基因预测及功能富集分析

miRNA的靶基因通过miRanda(v3.3a)[30]软件预测，参数设置为：-sc 140 -en -20 -scale 4 -strict。靶基因的GO (Gene Ontology)和KEGG (Kyoto Ency­clopedia of Genes and Genomes)富集通过DAVID[31]在线软件分析，经过超几何分布检验后，<0.05的term为显著富集。使用R语言程序绘制柱状图及气泡图。

1.6 miRNA-脂代谢靶基因互作网络分析

选取差异表达miRNA靶基因中与脂代谢相关的基因进行miRNA-靶基因互作网络分析，包括脂肪酸合成、脂肪酸降解、甘油脂代谢、甘油磷脂代谢以及鞘脂代谢等。使用Cytoscape(v3.3.0)软件绘制互作网络图。

1.7 qRT-PCR验证miRNA及其靶基因

对6条差异表达的miRNA进行茎环荧光定量qRT-PCR验证。为了进一步分析病毒诱导的miRNA是否可能参与调节宿主的脂质代谢，从miRNA靶基因中选取了6个与脂代谢相关的关键酶基因，即乙酰辅酶A羧化酶-1 (acetyl-CoA carboxylase-1,)、丝氨酸棕榈酰转移酶(serine palmitoyltransferase,)、酰基辅酶A氧化酶(acyl-CoA oxidase,)、乙酰辅酶A酰基转移酶(acetyl-CoA acyltransferase,)、鞘氨脂碱N-棕榈转移酶(sphingoid base N-stearoyltransferase,)以及丁酰辅酶A脱氢酶(butyryl-CoA dehydrogenase,)进行常规荧光定量qRT-PCR检测，并分析它们与miRNA表达水平之间的关系。用Primer 5软件设计引物(表1)，其中通用茎环引物序列为：5′-GTCGTATCCAGTGC­AGGGTCCGAGGTATTCGCACTGGATACG-3′，均由铂尚生物技术(上海)有限公司合成。提取的RNA经DNA酶消化后反转录为cDNA，反转录试剂盒为FastQuant RT Kit with gDNase试剂盒(TIANGEN，北京)。qRT-PCR反应体系参照Power SYBR® Green PCR Master Mix (Thermo Fisher)荧光定量检测试剂盒说明书。miRNA以U6为内参，靶基因以细胞周期蛋白激酶CDKA为内参[32]。qRT-PCR实验在QuantStudio 7 Flex Real Time PCR system(Thermo Fisher)上进行，每组实验设置3个生物学重复和3个技术重复(所用样本均为测序同批次藻样)，并采用2–ΔΔCt法进行相对定量[33]。

表1 miRNA及其靶基因的引物序列

1.8 统计学分析

采用SPSS 17.0软件进行统计分析。qRT-PCR的实验数据以平均值±标准差表示，显著性差异采用检验，<0.05表示差异显著，<0.01表示差异极显著。

2 结果与分析

2.1 赫氏颗石藻miRNA的基本特征

经Illumina HiSeqTM 2500平台高通量测序后，每个small RNA文库至少获得了两千万个原始reads，其中高质量的reads占比>97%，经过筛选后每个样本中符合条件的clean tags均在70%左右，满足后续分析的需要。原始数据已上传至NCBI的SRA数据库(SRP108676)。

通过对small RNA进行过滤、质控、分类以及与参考基因组比对后，共鉴定出26条成熟miRNA (其中包括两条病毒来源的miRNA)，23个前体miRNA，且病毒的2条miRNA来自同一个前体miRNA。66.7%的miRNA来自基因间区，29.2%的miRNA来自外显子区，只有1条miRNA来自内含子区。miRNA的长度在18～30 nt之间，峰值为21 nt (图1)。这些miRNA的拷贝数从8到765,680不等，表明在病毒感染条件下miRNA表达水平差异较大。此外，赫氏颗石藻miRNA的首位碱基对C有着较强的偏好性(54.17%)。成熟miRNA的序列、测序丰度以及在基因组的定位见附表1，前体miRNA序列见附表2。

2.2 赫氏颗石藻miRNA的保守性分析

用序列比对工具将预测的26条miRNA序列与miRBase数据库收录的miRNA以及未被录入miRBase数据库的微藻miRNA进行比对分析，结果见附表3。虽然本文中的一些miRNA能够和其他miRNA的部分碱基配对(包括来源于CCMP1516的18条miRNA)，但都不满足鉴定为保守miRNA的条件。赫氏颗石藻BOF92来源的前体miRNA及成熟miRNA与其他13株赫氏颗石藻基因组的比对结果显示(附表4、5)，10条miRNA (深灰色标记)与其他所有13株藻株基因组完全匹配，另外14条的前体序列或成熟体序列的匹配较短或无法匹配，表明这10条来源于赫氏颗石藻BOF92的miRNA及其前体比其他miRNAs更为保守。对于病毒的2条miRNA，它们都不能同时比对到其他13株赫氏颗石藻病毒的基因组上，表明EhV99B1编码的miRNA在株系间的变异更为明显，但ehv-miR1-3p比ehv-miR2-5p相对保守(附表6)。

图1 miRNA的长度分布

2.3 miRNA差异表达分析

对26条miRNA的测序丰度进行差异分析，结果显示共有3条宿主miRNA和2条病毒miRNA在本实验条件下发生了差异表达(表2)。这5条miRNA茎环荧光定量qRT-PCR结果与测序结果基本一致。5条差异miRNA的二级结构如图2所示，其中3条miRNA位于前体二级结构的5′端臂，2条位于3′端臂，且病毒的2条miRNA (ehv-miR1-5p和ehv-miR2- 3p)来自同一个前体二级结构。

2.4 差异表达miRNA靶基因的预测及GO富集分析

以同批次样品转录组测序结果中的32,909条宿主mRNA作为靶基因预测的数据库(NCBI登录号为SRP189555)[29]，26条miRNA共获得了26,380个靶基因。5条差异表达的miRNA 靶向了11,408个靶基因，其中3986个靶基因发生了差异表达(在感染早期和晚期分别为1246和3043个，交集为303个基因)。差异表达的靶基因GO富集结果主要分为3个类别(图3)：生物过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)。生物过程中富集到基因数目最多的是代谢过程(metabolic process)，还包括刺激应答(response to stimulus)、信号通路(signaling)等；细胞组分中膜(membrane)、膜组分(membrane part)以及细胞(cell)等富集最显著；分子功能中，更多的靶基因与催化活性(catalytic activity)、结合(binding)以及转运活性(transporter activity)相关，此外还包括信号转导活性(signal transducer activity)和抗氧化活性(antioxidant activity)。GO富集结果表明miRNA对颗石藻各种功能基因的调控范围十分广泛。

表2 差异表达miRNA的测序差异倍数对数值及qRT-PCR的–ΔΔCt结果

*表示<0.05；**表示<0.01；a表示Con_45 h样本测序结果中未检测到丰度值，TPM值以0.01代替；b表示Con_45 h样本qRT-PCR结果中值达到最大阈值50。

图2 5条差异差异表达miRNA的前体二级结构

在每个发夹结构中，红色表示位于5′臂(5p)的成熟miRNA，绿色表示位于3′臂(3p)的成熟miRNAs。

此外，本研究也对病毒的基因进行了靶向预测分析。结果显示，赫氏颗石藻 BOF92来源的成熟miRNA靶向8个有明确功能的病毒基因(附表7)，分别是病毒衣壳蛋白(major capsid protein)、DNA拓扑异构酶(DNA topoisomerase)、核酸内切酶(endo­nuclease)、核糖核酸酶(ribonuclease)、核糖核苷酸还原酶蛋白(ribonucleoside-diphosphate reductase protein)、DNA依赖的RNA聚合酶II亚基(DNA-dependent RNA polymerase II largest subunit)、DNA指导的RNA聚合酶II亚基(DNA-directed RNA polymerase II subunit)以及丝氨酸蛋白酶(serine protease)。其中，宿主编码的上调ehx-miR2-5p和ehx-miR3-3p靶向病毒衣壳蛋白基因的得分均较高，提示ehx-miR2-5p和ehx-miR3-3p可能在宿主抗病毒防御过程中具有一定的作用。衣壳蛋白的作用是包裹病毒的遗传物质，能够协助病毒感染，且具有免疫原性。而宿主编码的miRNA可能通过靶向病毒衣壳蛋白基因以抑制其翻译，进而影响病毒的装配过程。表明ehx-miR2-5p和ehx-miR3-3p的上调表达可能是宿主的抗病毒策略。由于转录组测序数据中没有检测到病毒的转录本，因此有关宿主miRNA如何影响病毒的生理过程本文将不做进一步分析。

图3 miRNA靶基因的GO富集分析

横坐标表示不同类别的GO term，纵坐标表示每个GO term中基因的数量，每个柱子上面的数字代表miRNA靶基因占该通路中总基因数的百分比(%)；不同的颜色对应不同的显著性值。

2.5 差异表达miRNA靶基因的KEGG富集分析

差异表达miRNA的靶基因KEGG富集结果见图4。在病毒感染早期，富集最显著的两个通路分别是戊糖葡糖醛酸相互转化(pentose and glucuronate interconversions)和氨酰tRNA合成(aminoacyl-tRNA biosynthesis)此外，还包括脂代谢过程，特别是甘油脂代谢(glycerolipid metabolism)和脂肪酸代谢(fatty acid metabolism)。甘氨酸，丝氨酸和苏氨酸代谢(glycine, serine and threonine metabolism)也富集显著，丝氨酸可作为神经酰胺合成的原料参与鞘脂代谢。在病毒感染晚期，富集最显著的代谢通路包括ABC转运蛋白(ABC transporters)和氨酰tRNA合成(aminoacyl-tRNA biosynthesis)，鞘脂代谢(sphingo­lipid metabolism)、脂肪酸降解(fatty degradation)以及糖酵解/糖异生(glycolysis/gluconeogenesis)等。特别是感染早期和晚期，差异表达miRNA的靶基因均能够显著富集到宿主脂代谢过程，暗示病毒感染重构的脂代谢过程可能存在转录后水平调控。

2.6 miRNA与脂代谢靶基因的互作网络

miRNA-脂代谢靶基因的互作网络图能够帮助人们更好地了解miRNA与脂代谢相关基因之间的调控关系。从图5中可以看出，ehx-miR3-3p和ehv- miR1-5p靶基因数量最多，并且这两个miRNA的靶基因在5个脂代谢相关通路中都有分布，表明这两个miRNA可能在脂代谢调控中起到核心作用。同时，部分基因也能被多个miRNA靶向，比如(长链酰基辅酶A合成酶4)能够同时被ehv-miR2-3p和ehx-miR3-3p靶向。此外，病毒miRNA和宿主miRNA也能同时靶向某一代谢途径，表明病毒可能会利用宿主miRNA来共同调节脂代谢过程。

图4 病毒感染6 h (A)和45 h (B)差异miRNA靶基因的KEGG富集分析

富集倍数为：(某通路中差异基因数目/有注释的差异基因总数)/(某通路中的有注释的基因总数/背景基因中有注释的基因总数)；气泡的大小对应该通路中基因数量的多少，颜色表示富集的显著性值。

图5 miRNA-脂代谢靶基因互作网络图

图中红色圆形表示miRNA，三角形表示脂代谢相关靶基因，不同颜色对应不同脂代谢过程。

2.7 miRNA靶基因qRT-PCR验证

为了进一步了解miRNA与脂代谢相关靶基因的靶向关系，采用qRT-PCR方法验证了miRNA靶基因的表达水平，并分析miRNA与其靶基因表达之间的相关性。选取6个脂代谢关键靶基因(附表8)，包括参与脂肪酸合成的，参与鞘脂代谢过程的和，以及参与脂肪酸降解过程的,和酶基因。荧光定量结果显示，除了基因在病毒感染早期表达水平显著升高外，其他几个脂代谢相关酶基因的表达水平在病毒感染早期(6 h)和晚期(45 h)均不同程度下调(附图1)；而且这些酶基因的表达水平与miRNA的qRT-PCR结果均为负相关关系，符合一般的miRNA对靶基因的调控规律，表明这些靶基因的表达水平很可能受到miRNA的调控。

3 讨论

miRNA是近年来RNA生物学领域中的重大发现，其通过结合在靶基因的不同位点抑制或激活靶基因的表达水平，作为丰度极高的基因调控因子参与调控多种生物学过程。特别是新近研究发现miRNA与它作用的靶基因是响应生物及非生物环境胁迫的主要调控因子[34]。例如，miRNA作为关键因子在复杂的病毒–宿主互作网络中调控病毒的复制、免疫逃逸及宿主抗病毒等过程，并可能成为病毒感染诊断及治疗的新手段[35]。miRNA的另一个新功能是通过控制代谢过程中的靶基因，从而调控相关代谢通路，特别是在脂质代谢网络中起重要调节作用[19,36]。越来越多的研究表明，病毒感染赫氏颗石藻能够显著重塑宿主脂肪酸代谢、鞘脂代谢及甘油脂等脂代谢过程，且这些相关脂代谢与病毒的侵染、复制、组装及病毒颗粒的释放等过程密切相关[9～11]。miRNA是否参与颗石藻病毒感染重塑的脂代谢过程尚不清楚。本研究通过对EhV99B1感染的赫氏颗石藻BOF92细胞样品进行小RNA测序，鉴定出26条新的miRNA，分析其基本特征和保守性，并结合miRNA差异表达谱特征、靶基因功能富集以及qRT-PCR实验等，初步分析了病毒感染诱导的差异表达miRNA在宿主脂代谢过程中可能发挥的作用。

3.1 赫氏颗石藻miRNA的基本特征

本研究经测序分析共鉴定出了24条颗石藻miRNA，长度主要集中于18～21 nt，峰值为21 nt，符合一般植物miRNA的特征。BOF92的miRNA首位碱基对C有着较强的偏好，这一点与营养胁迫条件下CCMP1516中的miRNA相似[28]。然而，在大部分生物体中，包括其他一些浮游植物，其miRNA的5′端首位碱基通常为U[37～41]。这可能与参与RISC (RNA诱导的沉默复合体)形成过程中的2个重要蛋白即Dicer-like protein (DCL)和 Argonaute (AGO)有关。DCL能够切割初始miRNA进而形成前体miRNA，该蛋白通常包含以下几个特征结构域：DEAD、Hleicase-C、dsRBD及RNA酶Ⅲ结构域[42,43]。CCMP1516能够编码四种DCL，但是它们都缺少一个或多个特征结构域，有三种包含了DEAD，Helicase-C和DSRM结构域，但不含有RNA酶Ⅲ结构域，该结构与对于RNA的绑定和切割十分重要；另外一种DCL含有DRSM和RNA酶Ⅲ结构域，但不含DEAD和Helicase-C结构域[25]，颗石藻DCL结构上的这些差异可能会影响其切割位点的选择，进而影响miRNA首位碱基的偏好性。AGO蛋白位于RISC的中心，其最主要的功能结构域为PAZ和PIWI结构域。CCMP1516基因组中有2种AGO同源物，虽然它们都含有PAZ和PIWI结构域，但与其他微藻以及部分模式物种的AGO系统发育树的分析结果显示[28]，包括颗石藻在内的各种微藻的AGO蛋白都各自聚集成一类，表明这些微藻的AGO蛋白在进化上同源性较低。因此颗石藻可能也存在特殊的RNA沉默机制，使得miRNA的首位碱基不偏好于U。

在病毒诱导的BOF92 miRNA中，66.7%来自基因间区，29.2%来自外显子区，只有1条miRNA来自内含子区(附表1)。这与近缘株系CCMP1516的miRNA来源分布较为相似，CCMP1516的miRNA约77%来自基因间区，23%来自内含子区[28]。三角褐指藻()和团藻()的miRNA也都主要来源于基因间区，只有少部分miRNA比对到了外显子区[38,44]，这些微藻的miRNA的来源更类似于陆生高等植物[18]。莱茵衣藻()和长囊水云()的miRNA来源类似，即大部分位于内含子区，少数位于基因间区[40,39]。综上所述，不同浮游植物miRNA的来源具有多样性，表明miRNA在浮游植物的转录后水平可能具有广泛的调控功能。

3.2 miRNA保守性特点

针对本研究中鉴定的24条宿主miRNA同源性较低的问题，推测可能有以下几个原因。首先，不同微藻之间的miRNA可能是独立进化的。Zhao等[40]首先对莱茵衣藻的miRNA进行了表征，作者将鉴定出来的莱茵衣藻miRNA与miRBase数据库比对后未发现同源miRNA。随后多位学者对不同种的微藻，如假微型海链藻()[21]，三角褐指藻()[44]及赫氏颗石藻CCMP1516[28]等在不同生长条件或环境胁迫条件下的miRNA进行了研究，也均未能在miRBase数据库找到同源的序列。这暗示miRNA可能在浮游植物、动植物中分别有着独立的进化机制。其次，病毒感染诱导产生特异性miRNA的表达。miRNA的表达是动态的，具有时空特异性，不同胁迫条件如营养胁迫、氧化应激、病毒感染等都会影响miRNA的表达，某些特异的miRNA可能表达水平很低，仅在特殊的生理生态条件下表达量才增加[28]。因此，病毒感染可能使赫氏颗石藻BOF92产生了特异性miRNA。另外，赫氏颗石藻miRNA株系间的变异较为明显。根据本文24条宿主miRNA与其他13株赫氏颗石藻基因组的比对结果看，能够全部比对上其他赫氏颗石藻株系的miRNA仅有10条，占miRNA总数的42%左右。赫氏颗石藻作为一种单细胞浮游植物，基因组的种内变异更有利于其适应动态变化的海洋环境以及在与病毒的“arms race”中获得优势。本研究中鉴定的24条宿主miRNA保守性较低可能是这些因素共同作用的结果。此外，本研究也没有找到EhV99B1 miRNA的同源序列。一方面，已鉴定的病毒miRNA非常有限；另一方面，与真核生物相比，病毒基因组有着更高的变异率，这意味着病毒拥有能快速适应宿主和环境条件的进化优势。

3.3 miRNA靶向宿主脂代谢途径

miRNA靶基因功能富集结果显示，病毒感染早期(6 h)和晚期(45 h)，差异表达的miRNA均能靶向宿主脂肪酸和鞘脂代谢通路(图6)。就脂肪酸代谢而言，感染早期，宿主编码的ehx-miR1-5p下调，其对脂肪酸合成过程的限速酶靶基因的抑制作用减弱，与此同时宿主编码的ehx-miR4-3p上调，其对脂肪酸降解过程中的关键酶和靶基因的抑制作用增强，从而导致脂肪酸含量的积累(图6A)。感染晚期，则是病毒编码的ehv-miR2-5p 上调，并靶向抑制脂肪酸降解过程的相关酶基因(图6B)。因此，推测感染过程中，宿主和病毒编码的miRNA协同作用导致宿主细胞中脂肪酸含量的不断积累。脂肪酸在病毒感染过程有着十分重要的作用，一方面大型双链DNA病毒的复制和装配有更高的代谢需求，另一方面脂肪酸可作为甘油脂和鞘脂的合成原料[7]。另外，病毒感染过程中，宿主鞘脂从头合成过程中的两个关键酶基因和的表达分别受到早期宿主上调的ehx-miR3-5p (图6A)和晚期病毒上调的ehv-miR2-5p (图6B)的靶向抑制，因而宿主鞘脂从头合成途径在整个病毒感染过程中被显著抑制，并导致鞘脂类物质含量减少。上述结果可能在一定程度上补充解释了本课题组之前的脂质组学数据，即EhV99B1感染早期(6 h)和晚期(45 h)，宿主赫氏颗石藻BOF92细胞中大部分脂肪酸代谢物显著积累，而几乎所有检测到的16种鞘脂类物质的含量均显著降低[11]。众所周知，鞘脂是细胞膜上“脂筏”结构域的关键功能成分，也是病毒包膜结构的主要成分之一，对病毒的侵染、装配及释放过程十分重要，且作为信号分子，在感染晚期能够诱导宿主程序性细胞死亡。对赫氏颗石藻CCMP1516及其特异性裂解病毒EhV86的全基因测序注释结果显示，通过基因水平转移，病毒基因组“截获”了一套鞘脂从头合成途径中的关键酶基因，并通过调节相关基因的表达在一定程度上掌控了宿主的鞘脂代谢过程[45]。在EhV201感染赫氏颗石藻CCMP374过程中发现，随着感染时间的延长(感染0～36 h)，宿主基因组编码的基因表达显著降低，而病毒基因组编码的在mRNA和蛋白水平的表达均明显上升，从而大量合成并积累病毒性鞘糖脂[9]，并最终诱导宿主程序性细胞死亡[15,17]。综上所述，病毒和宿主的miRNA均可能参与宿主脂代谢的重塑过程，其中病毒的miRNA主要在感染晚期通过对脂肪酸降解过程和鞘脂合成过程中的部分基因进行靶向抑制，进而促进宿主脂肪酸的积累和鞘脂水平的降低。

图6 miRNA靶向的脂肪酸及鞘脂代谢通路

A和B分别表示病毒感染6 h和45 h差异miRNA对宿主脂肪酸和鞘脂代谢的靶向调控图。红色字体代表miRNA及靶基因高表达，蓝色填充表示miRNA及靶基因低表达；红色箭头代表脂代谢物积累，蓝色箭头表示脂代谢物减少。虚线框中的绿色字体表示不同的代谢通路

附加材料详见文章电子版www.chinagene.cn。

感谢Gunnar Bratbak教授(挪威卑尔根大学生物系)提供的BOF92藻株和virus 99B1病毒株。

[1] Tong SY, Gao KS, Hutchins DA. Adaptive evolution in the coccolithophorefollowing 1,000 generations of selection under elevated CO2., 2018, 24(7): 3055–3064.

[2] Sun J. Organic carbon pump and carbonate counter pump of living coccolithophorid., 2007, 22(12): 1231–1239.

孙军. 今生颗石藻的有机碳泵和碳酸盐反向泵. 地球科学进展，2007, 22(12): 1231–1239.

[3] Zhang J, Li JR, Yang L, Yue NN. Coccolithophores and their characteristics., 2018, 42(2): 150–158.

张健, 李佳芮, 杨璐, 岳娜娜. 球石藻及其生态功能. 海洋科学, 2018, 42(2): 150–158.

[4] Daniels CJ, Poulton AJ, Balch WM, Marañón E, Adey T, Bowler BC, Cermeño P, Charalampopoulou A, Crawford DW, Drapeau D, Feng YY, Fernández A, Fernández E, Fragoso GM, González N, Graziano LM, Heslop R, Holligan PM, Hopkins J, Huete-Ortega M, Hutchins DA, Lam PJ, Lipsen MS, López-Sandoval DC, Loucaides S, Marchetti A, Mayers KMJ, Rees AP, Sobrino C, Tynan E, Tyrrell, T. A global compilation of coccolithophore calcification rates., 2018, 10(4): 1859–1876.

[5] Ruiz E, Oosterhof M, Sandaa RA, Larsen A, Pagarete A. Emerging interaction patterns in the- EhV system., 2017, 9(3): 61.

[6] Highfield A, Evans C, Walne A, Miller PI, Schroeder DC. How many coccolithovirus genotypes does it take to terminate anbloom?, 2014, 466–467: 138–145.

[7] Rosenwasser S, Mausz MA, Schatz D, Sheyn U, Malitsky S, Aharoni A, Weinstock E, Tzfadia O, Ben-Dor S, Feldmesser E, Pohnert G, Vardi A. Rewiring host lipid metabolism by large viruses determines the fate of, a bloom-forming alga in the ocean., 2014, 26(6): 2689–2707.

[8] Rosenwasser S, Ziv C, Creveld SG, Vardi A. Virocell metabolism: metabolic innovations during host–virus interactions in the ocean., 2016, 24(10): 821–832

[9] Ziv C, Malitsky S, Othman A, Ben-Dor S, Wei Y, Zheng SN, Aharoni A, Hornemann T, Vardi A. Viral serine palmitoyltransferase induces metabolic switch in sphingolipid biosynthesis and is required for infection of a marine alga., 2016, 113(13): E1907–E1916.

[10] Evans C, Pond DW, Wilson WH. Changes infatty acid profiles during infection with EhV86: physiological and ecological implications., 2009, 55(3): 219–228.

[11] Zeng J, Liu SSY, Cai WC, Jiang HR, Lu X, Li GL, Li J, Liu JW. Emerging lipidome patterns associated with marine-virus model system., 2019, 688: 521–528.

[12] Schatz D, Rosenwasser S, Malitsky S, Wolf SG, Feldmesser E, Vardi A. Communicationextracellular vesicles enhances viral infection of a cosmopolitan alga., 2017, 2(11): 1485–1492.

[13] Sheyn U, Rosenwasser S, Ben-Dor S, Porat Z, Vardi A. Modulation of host ROS metabolism is essential for viral infection of a bloom-forming coccolithophore in the ocean., 2016, 10(7): 1742–1754.

[14] Vardi A, Haramaty L, Van Mooy BAS, Fredricks HF, Kimmance SA, Larsen A, Bidle KD. Host-virus dynamics and subcellular controls of cell fate in a natural coccolithophore population., 2012, 109(47): 19327–19332.

[15] Bidle KD, Haramaty L, Barcelos E Ramos J, Falkowski P. Viral activation and recruitment of metacaspases in the unicellular coccolithophore., 2007, 104(14): 6049–6054.

[16] Vardi A, Van Mooy BAS, Fredricks HF, Popendorf KJ, Ossolinski JE, Haramaty L, Bidle KD. Viral glycosphingolipids induce lytic infection and cell death in marine phytoplankton., 2009, 326(5954): 861–865.

[17] Liu JW, Cai WC, Fang X, Wang XT, Li GL. Virus-induced apoptosis and phosphorylation form of metacaspase in the marine coccolithophorid., 2018, 200(3): 413–422.

[18] Millar AA, Waterhouse PM. Plant and animal microRNAs: similarities and differences., 2005, 5(3): 129–135.

[19] Aryal B, Singh AK, Rotllan N, Price N, Fernández- Hernando N. MicroRNAs and lipid metabolism., 2017, 28(3): 273–280.

[20] Lopez-Gomollon S, Beckers M, Rathjen T, Moxon S, Maumus F, Mohorianu I, Moulton V, Dalmay T, Mock T. Global discovery and characterization of small non-coding RNAs in marine microalgae., 2014, 15(1): 697.

[21] Norden-Krichmar TM, Allen AE, Gaasterland T, Hildebrand M. Characterization of the small RNA transcriptome of the diatom,., 2011, 6(8): e22870.

[22] Shi XG, Lin X, Li L, Li MZ, Palenik B, Lin SJ. Transcriptomic and microRNAomic profiling reveals multi-faceted mechanisms to cope with phosphate stress in a dinoflagellate., 2017, 11(10): 2209–2218.

[23] Tian X, Cai WC, Su JJ, Wu SY, Liu JW. Transcriptome analysis of marine microalgain response to virus infection., 2019, 41(12): 103–112.

田雪, 蔡伟聪, 苏金净, 吴书燕, 刘静雯. 病毒感染海洋球石藻的转录组分析. 海洋学报, 2019, 41(12): 103–112.

[24] Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads., 2011, 17(1): 10–12.

[25] Hällman J, Avesson L, Reimegård J, Käller M, Söderbom F. Identification and verification of microRNAs by high-throughput sequencing., 2013, 983: 125–138.

[26] Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades., 2012, 40(1): 37–52.

[27] Robinson M, Mccarthy D, Smyth G. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data., 2009, 26(1): 139–140.

[28] Zhang XY, Gamarra J, Castro S, Carrasco E, Hernandez A, Mock T, Hadaegh AR, Read BA. Characterization of the small RNA transcriptome of the marine coccolithophorid,., 2016, 11(4): e0154279.

[29] Read BA, Kegel J, Klute MJ, Kuo A, Lefebvre SC, Maumus F, Mayer C, Miller J, Monier A, Salamov A, Young J, Aguilar M, Claverie JM, Frickenhaus S, Gonzalez K, Herman EK, Lin YC, Napier J, Ogata H, Sarno AF, Shmutz J, Schroeder D, de Vargas C, Verret F, von Dassow P, Valentin K, Van de Peer Y, Wheeler G,Annotation Consortium, Dacks JB, Delwiche CF, Dyhrman ST, Glöckner G, John U, Richards T, Worden AZ, Zhang X, Grigoriev IV. Pan genome of the phytoplanktonunderpins its global distribution., 2013, 499(7457): 209–213.

[30] Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in drosophila.2003, 5(1): R1.

[31] Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources., 2009, 4(1): 44–57.

[32] Zhang E, Wu S, Cai W, Zeng J, Li J, Li G, Liu J. Validation of superior reference genes for qRT-PCR and western blot analyses in marine-virus model system., 2021, 131(1): 257–271.

[33] Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCTMethod., 2001, 25(4): 402–408.

[34] Manavella1 PA, Yang SW, Palatnik J. Keep calm and carry on: miRNA biogenesis under stress., 2019, 99(5): 832–843.

[35] Bruscella P, Bottini S, Baudesson C, Pawlotsky JM, Feray C, Trabucchi M. Viruses and miRNAs: more friends than foes., 2017, 8: 724.

[36] Li C, Lu LG, Qi ZT, Zhu YQ, Su F, Zhao P, Dong H. Transcriptome and miRNome analysis provide new insight into host lipid accumulation, innate immunity, and viral persistence in hepatitis C virus infection., 2020, 11: 535673.

[37] Aravin AA, Lagos-Quintana M, Yalcin A, Zavolan M, Marks D, Snyder B, Gaasterland T, Meyer J, Tuschl T. The small RNA profile duringdevelopment., 2003, 5(2): 337–350.

[38] Li JR, Wu Y, Qi YJ. MicroRNAs in a multicellular green alga., 2014, 57(1): 36–45.

[39] Tarver JE, Cormier A, Pinzón N, Taylor RS, Carré W, Strittmatter M, Seitz H, Coelho SM, Cock JM. MicroRNAs and the evolution of complex multicellularity: identification of a large, diverse complement of microRNAs in the brown alga., 2015, 43(13): 6384–6398.

[40] Zhao T, Li GL, Mi SJ, Li S, Hannon GJ, Wang XJ, Qi YJ. A complex system of small RNAs in the unicellular green alga., 2007, 21(10): 1190–1203.

[41] Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function., 2004, 116(2): 281–297.

[42] Ha MJ, Kim VN. Regulation of microRNA biogenesis., 2014, 15(8): 509–524.

[43] Peng JJ, Yan F, Chen HR, Chen JP. Progress of studies on Dicer structure and function., 2008, 30(12): 1550–1556.

彭杰军, 燕飞, 陈海如, 陈剑平. Dicer结构和功能研究进展. 遗传, 2008, 30(12): 1550–1556.

[44] Huang AY, He LW, Wang GC. Identification and characterization of microRNAs fromby high-throughput sequencing and bioinformatics analysis., 2011, 12(1): 337.

[45] Wilson WH, Schroeder DC, Allen MJ, Holden M, Parkhill J, Barrell BG, Churcher C, Harnlin N, Mungall K, Norbertczak H, Quail MA, Price C, Rabbinowitsch E, Walker D, Craigon M, Roy D, Ghazal P. Complete genome sequence and lytic phase transcription profile of a coccolithovirus., 2005, 309(5737): 1090LP–1092.

Analysis of microRNA expression profile inin response to virus infection

Enquan Zhang, Weicong Cai, Guiling Li, Jian Li, Jingwen Liu

The interactions betweenandvirus (EhV) regulate marine carbon and sulfur biogeochemical cycle and play a prominent role in global climate change. As a large DNA virus, EhVs have developed a novel “virocell metabolism” model to meet their higher metabolic needs. However, the regulatory mechanism of this metabolic model is still largely unclear. MicroRNAs (miRNAs) can regulate biological pathways through targeting hub genes in the metabolic processes. Here, we performed high-throughput small RNA sequencing to analyse miRNA expression in EhV99B1 infectedBOF92. A total of 26 miRNAs (including 2 virus-derived miRNAs) were identified, including four up-regulated and one down-regulated miRNAs. These results were further validated through quantitative real-time PCR. Functional enrichment analysis showed that five differentially-expressed miRNAs might be involved in the regulation of carbohydrate metabolism, lipid metabolism and amino acid metabolism. Moreover, the expression levels of differentially-expressed miRNAs were negatively correlated with that of several lipid metabolism-related genes, such as,,,,and, indicating that these miRNAs might play an important regulatory role in virus-mediated lipid metabolism.

;virus; small RNA sequencing; microRNA; lipid metabolism

2021-04-30;

2021-07-14

国家自然科学基金面上项目(编号：42076086，31771972)和福建省自然科学基金(编号：2019J01696，2020J01676)资助[Supported by National Natural Science Foundation of China (Nos. 42076086, 31771972) and Fujian Province Natural Science Foundation of China (Nos. 2019J01696, 2020J01676)]

张恩权，在读硕士研究生，专业方向：海洋微生物组学。E-mail: chriszhangen@163.com

刘静雯，博士，教授，研究方向：海洋微型生物分子生态学。E-mail: ljwsbch@163.com

10.16288/j.yczz.21-164

2021/7/30 11:38:01

URI:https://kns.cnki.net/kcms/detail/11.1913.r.20210730.0837.001.html

(责任编委: 赵呈天)