張恩權，蔡偉聰，李桂玲，李健，劉靜雯

研究報告

赫氏顆石藻()響應病毒感染的microRNA轉錄組分析

張恩權，蔡偉聰，李桂玲，李健，劉靜雯

集美大學食品與生物工程學院，廈門 361000

海洋顆石藻病毒–宿主互作是影響海洋碳、硫生物地化循環(huán)及全球氣候變化的重要環(huán)節(jié)。作為大的雙鏈DNA病毒，顆石藻病毒進化出一種“病毒細胞代謝”模式，通過重編程宿主代謝途徑以滿足其代謝需求，但對這一代謝模式的調控機制尚缺乏足夠的認識。MicroRNA (miRNA)作為一種基因表達調控的重要因子，能夠通過調控代謝過程中的靶基因表達，從而調節(jié)相關代謝通路。本研究采用small RNA測序技術分析病毒感染顆石藻差異表達的miRNA及其靶基因功能，鑒定出26條成熟miRNA (包括2條病毒來源的miRNA)，均來自23條新的miRNA前體序列，其中5條miRNA顯著差異表達，包括4條上調，1條下調。實時熒光定量PCR (quantitative real-time PCR, qRT-PCR)驗證結果與miRNA-seq結果基本一致。功能富集分析顯示，5個差異表達的miRNA可能參與調節(jié)糖代謝、脂代謝和氨基酸等代謝。此外，差異表達miRNA的表達水平與脂質代謝相關靶基因如、、、、、等的表達水平呈負相關，說明這些miRNA可能在病毒感染過程中對宿主的脂質代謝發(fā)揮重要的調控作用。

赫氏顆石藻；顆石藻病毒；small RNA測序；microRNA；脂代謝

海洋顆石藻(Coccolithophores)是一種全球廣泛分布且具有重要生態(tài)功能的真核微型浮游植物[1～3]，其中赫氏顆石藻(, Eh)具有形成“球石?！焙透弋a(chǎn)二甲基硫化物(DMSP)的能力，且?guī)缀趺磕甓荚诖笱?尤其是高緯度海域)中形成大面積赤潮[4]，該藻赤潮的迅速大規(guī)模消亡被證實是特異性病毒(virus, EhV)感染和裂解所致[5,6]。因此，EhV-Eh的互作過程是影響海洋碳、硫生物地化循環(huán)及全球氣候變化的重要環(huán)節(jié)。顆石藻病毒與宿主間侵染和抵抗的博弈形成了一種新的、穩(wěn)定的協(xié)同進化代謝模式，是研究真核浮游植物宿主-病毒互作的理想模式系統(tǒng)[7]。

脂代謝是當今病毒-宿主互作研究的重要問題。作為大的雙鏈DNA病毒，EhVs進化出一種全新的“病毒細胞代謝(virocell metabolism)”模式，通過重編程宿主代謝途徑以滿足其更高的代謝需求[8,9]。Evans等[10]2009年首次報道在赤潮消亡過程中，病毒感染誘導宿主細胞脂肪酸的組成由多不飽和向單不飽和轉變，導致被病毒感染的宿主細胞減少多不飽和脂肪酸向食物鏈中更高級營養(yǎng)水平的傳遞，從而降低海洋生態(tài)系統(tǒng)的總生產(chǎn)力。本課題組前期研究發(fā)現(xiàn)，病毒感染早期宿主細胞高飽和三?；视?triacylglycerols, TAGs)合成積累、并形成脂滴聚集在病毒顆粒中促進病毒外殼疏水蛋白的存儲，滿足病毒復制和組裝的需求[11]。最新研究發(fā)現(xiàn)，病毒感染顆石藻能夠以外泌體形式富集并分泌TAGs，以此加速病毒感染過程[12]。以色列學者Carmit團隊2016年首次發(fā)現(xiàn)，EhV編碼的絲氨酸棕櫚酰轉移酶(viral serine palmitoyltransferase,：鞘脂從頭生物合成途徑中的第一個關鍵限速酶)在感染細胞中能夠催化合成病毒特有的新型鞘糖脂(viral-specific glycosphingolipids, vGSLs)，該物質被認為是病毒組裝的必要條件[9]，這表明病毒能夠利用編碼的輔助代謝基因(virus-encoded auxiliary metabolic genes, vAMGs)重塑宿主鞘脂代謝途徑。同時vGSLs還能作為重要的信號分子誘發(fā)ROS(H2O2)[13]和NO[14]的產(chǎn)生，進而啟動宿主細胞凋亡程序，最終裂解細胞釋放病毒粒子[15～17]。然而，目前對EhV感染重塑的脂代謝的miRNA調控機制尚缺乏足夠的認識。

MicroRNA(miRNA)是一種長度分布在18～26 nt的內源性非編碼小RNA。動植物miRNA來自基因組位置略有不同，動物miRNA通常位于內含子區(qū)，而植物miRNA更多地來自基因組的基因間區(qū)[18]。miRNA通過結合在靶基因的不同位點抑制或激活靶基因的表達水平，作為豐度極高的基因調控因子參與調控多種生物學過程。最近在高等動植物研究中發(fā)現(xiàn)，miRNA能夠通過調節(jié)代謝過程中的靶基因，從而調控相關代謝通路，特別是在脂質合成、脂肪酸氧化及脂蛋白分泌等代謝網(wǎng)絡中起重要調節(jié)作用[19]。相對于高等動植物而言，由于浮游植物種內和種間的差異性特點，且目前只有少數(shù)種類基因組信息得到較完整的詮釋，因此對miRNA在海洋浮游植物生物及非生物脅迫中的重要作用了解較少。目前，主要針對氮、磷、硅及鐵等營養(yǎng)脅迫條件下miRNA的表達情況進行了初步研究[20～22]。有關浮游植物響應病原菌感染相關miRNA的研究國內外尚未報道。

本研究以Eh-EhV模式系統(tǒng)為研究對象，在病毒感染早期(6 h)和晚期(45 h)收集藻細胞樣品，采用高通量小RNA測序技術篩選病毒誘導的差異表達miRNA，分析miRNA在病毒重塑宿主脂代謝過程中可能的調控作用，為從表觀遺傳學角度深入理解EhV-Eh互作的分子機制提供新的理論認識。

1 材料與方法

1.1 赫氏顆石藻的培養(yǎng)及病毒感染

赫氏顆石藻BOF92及其特異性裂解病毒株系EhV99B1均由挪威卑爾根大學生物系微生物研究所Gunnar Bratbak教授饋贈并保存于本實驗室，藻株和病毒株均為純化株系。顆石藻培養(yǎng)采用70%海水配制的f/2-si加富培養(yǎng)基，培養(yǎng)條件為：溫度16℃±2℃，光照強度為60 μmol quanta m–2s–1，光照周期為14/10 (光/暗)。在2 L指數(shù)生長期的藻液中(～106cells/mL)，以1∶50 (EhV∶Eh)體積比加入濃縮病毒裂解液(病毒初始濃度約為107viruses/mL)，另外2 L添加等量高溫滅活的病毒作為對照組?；诒菊n題組前期mRNA轉錄組和脂質代謝組學分析結果，EhV感染早期(6 h)和中后期(45 h)，宿主轉錄組發(fā)生了顯著差異表達[23]，細胞脂質代謝產(chǎn)物的積累也發(fā)生了明顯變化[11]，因此本研究選擇病毒感染的6 h和45 h作為采樣時間點。分別于病毒感染6 h和45 h離心收集500 mL藻細胞樣品(4℃，7000 r/min，5min)，立即置于液氮中速凍，于?80℃保存?zhèn)溆?。每個樣本設置兩個生物學平行，共計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系統(tǒng)和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。后續(xù)所有分析基于clean reads。篩選過程基于Cutadapt 1.7.1[24]和Fastx toolkit 0.0.14[25]軟件。Clean reads經(jīng)過比對Rfam 11.0數(shù)據(jù)庫，過濾掉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數(shù)據(jù)庫收錄的38,589條miRNAs、其他微藻的2650條(數(shù)據(jù)庫尚未收錄的miRNAs)以及在CCMP1516中已鑒定出的18條miRNA[28]進行比對分析。鑒定保守miRNA的條件為：miRNA種子序列(第2到第8個堿基)不允許有錯配，且總錯配數(shù)≤2。由于顆石藻基因組具有明顯的種內變異，為了進一步了解BOF92的miRNAs及其前體的功能和進化情況，分析了miRNA及其前體序列與其他13株測序顆石藻株的保守性[29]。其中，92A、Eh2和Van556為深度測序的分離株，其他10株是低覆蓋度測序(覆蓋度可以達到91%至 95%)。對病毒來源的2條miRNA及其前體序列也與其他13株測序病毒株進行保守性分析。

1.5 miRNA靶基因預測及功能富集分析

miRNA的靶基因通過miRanda(v3.3a)[30]軟件預測，參數(shù)設置為：-sc 140 -en -20 -scale 4 -strict。靶基因的GO (Gene Ontology)和KEGG (Kyoto Ency-clopedia of Genes and Genomes)富集通過DAVID[31]在線軟件分析，經(jīng)過超幾何分布檢驗后，<0.05的term為顯著富集。使用R語言程序繪制柱狀圖及氣泡圖。

1.6 miRNA-脂代謝靶基因互作網(wǎng)絡分析

選取差異表達miRNA靶基因中與脂代謝相關的基因進行miRNA-靶基因互作網(wǎng)絡分析，包括脂肪酸合成、脂肪酸降解、甘油脂代謝、甘油磷脂代謝以及鞘脂代謝等。使用Cytoscape(v3.3.0)軟件繪制互作網(wǎng)絡圖。

1.7 qRT-PCR驗證miRNA及其靶基因

對6條差異表達的miRNA進行莖環(huán)熒光定量qRT-PCR驗證。為了進一步分析病毒誘導的miRNA是否可能參與調節(jié)宿主的脂質代謝，從miRNA靶基因中選取了6個與脂代謝相關的關鍵酶基因，即乙酰輔酶A羧化酶-1 (acetyl-CoA carboxylase-1,)、絲氨酸棕櫚酰轉移酶(serine palmitoyltransferase,)、?；o酶A氧化酶(acyl-CoA oxidase,)、乙酰輔酶A?；D移酶(acetyl-CoA acyltransferase,)、鞘氨脂堿N-棕櫚轉移酶(sphingoid base N-stearoyltransferase,)以及丁酰輔酶A脫氫酶(butyryl-CoA dehydrogenase,)進行常規(guī)熒光定量qRT-PCR檢測，并分析它們與miRNA表達水平之間的關系。用Primer 5軟件設計引物(表1)，其中通用莖環(huán)引物序列為：5′-GTCGTATCCAGTGC-AGGGTCCGAGGTATTCGCACTGGATACG-3′，均由鉑尚生物技術(上海)有限公司合成。提取的RNA經(jīng)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 統(tǒng)計學分析

采用SPSS 17.0軟件進行統(tǒng)計分析。qRT-PCR的實驗數(shù)據(jù)以平均值±標準差表示，顯著性差異采用檢驗，<0.05表示差異顯著，<0.01表示差異極顯著。

2 結果與分析

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

經(jīng)Illumina HiSeqTM 2500平臺高通量測序后，每個small RNA文庫至少獲得了兩千萬個原始reads，其中高質量的reads占比>97%，經(jīng)過篩選后每個樣本中符合條件的clean tags均在70%左右，滿足后續(xù)分析的需要。原始數(shù)據(jù)已上傳至NCBI的SRA數(shù)據(jù)庫(SRP108676)。

通過對small RNA進行過濾、質控、分類以及與參考基因組比對后，共鑒定出26條成熟miRNA (其中包括兩條病毒來源的miRNA)，23個前體miRNA，且病毒的2條miRNA來自同一個前體miRNA。66.7%的miRNA來自基因間區(qū)，29.2%的miRNA來自外顯子區(qū)，只有1條miRNA來自內含子區(qū)。miRNA的長度在18～30 nt之間，峰值為21 nt (圖1)。這些miRNA的拷貝數(shù)從8到765,680不等，表明在病毒感染條件下miRNA表達水平差異較大。此外，赫氏顆石藻miRNA的首位堿基對C有著較強的偏好性(54.17%)。成熟miRNA的序列、測序豐度以及在基因組的定位見附表1，前體miRNA序列見附表2。

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

用序列比對工具將預測的26條miRNA序列與miRBase數(shù)據(jù)庫收錄的miRNA以及未被錄入miRBase數(shù)據(jù)庫的微藻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在本實驗條件下發(fā)生了差異表達(表2)。這5條miRNA莖環(huán)熒光定量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作為靶基因預測的數(shù)據(jù)庫(NCBI登錄號為SRP189555)[29]，26條miRNA共獲得了26,380個靶基因。5條差異表達的miRNA 靶向了11,408個靶基因，其中3986個靶基因發(fā)生了差異表達(在感染早期和晚期分別為1246和3043個，交集為303個基因)。差異表達的靶基因GO富集結果主要分為3個類別(圖3)：生物過程(biological process)、細胞組分(cellular component)和分子功能(molecular function)。生物過程中富集到基因數(shù)目最多的是代謝過程(metabolic process)，還包括刺激應答(response to stimulus)、信號通路(signaling)等；細胞組分中膜(membrane)、膜組分(membrane part)以及細胞(cell)等富集最顯著；分子功能中，更多的靶基因與催化活性(catalytic activity)、結合(binding)以及轉運活性(transporter activity)相關，此外還包括信號轉導活性(signal transducer activity)和抗氧化活性(antioxidant activity)。GO富集結果表明miRNA對顆石藻各種功能基因的調控范圍十分廣泛。

表2 差異表達miRNA的測序差異倍數(shù)對數(shù)值及qRT-PCR的–ΔΔCt結果

*表示<0.05；**表示<0.01；a表示Con_45 h樣本測序結果中未檢測到豐度值，TPM值以0.01代替；b表示Con_45 h樣本qRT-PCR結果中值達到最大閾值50。

圖2 5條差異差異表達miRNA的前體二級結構

在每個發(fā)夾結構中，紅色表示位于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可能在宿主抗病毒防御過程中具有一定的作用。衣殼蛋白的作用是包裹病毒的遺傳物質，能夠協(xié)助病毒感染，且具有免疫原性。而宿主編碼的miRNA可能通過靶向病毒衣殼蛋白基因以抑制其翻譯，進而影響病毒的裝配過程。表明ehx-miR2-5p和ehx-miR3-3p的上調表達可能是宿主的抗病毒策略。由于轉錄組測序數(shù)據(jù)中沒有檢測到病毒的轉錄本，因此有關宿主miRNA如何影響病毒的生理過程本文將不做進一步分析。

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

橫坐標表示不同類別的GO term，縱坐標表示每個GO term中基因的數(shù)量，每個柱子上面的數(shù)字代表miRNA靶基因占該通路中總基因數(shù)的百分比(%)；不同的顏色對應不同的顯著性值。

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)也富集顯著，絲氨酸可作為神經(jīng)酰胺合成的原料參與鞘脂代謝。在病毒感染晚期，富集最顯著的代謝通路包括ABC轉運蛋白(ABC transporters)和氨酰tRNA合成(aminoacyl-tRNA biosynthesis)，鞘脂代謝(sphingo-lipid metabolism)、脂肪酸降解(fatty degradation)以及糖酵解/糖異生(glycolysis/gluconeogenesis)等。特別是感染早期和晚期，差異表達miRNA的靶基因均能夠顯著富集到宿主脂代謝過程，暗示病毒感染重構的脂代謝過程可能存在轉錄后水平調控。

2.6 miRNA與脂代謝靶基因的互作網(wǎng)絡

miRNA-脂代謝靶基因的互作網(wǎng)絡圖能夠幫助人們更好地了解miRNA與脂代謝相關基因之間的調控關系。從圖5中可以看出，ehx-miR3-3p和ehv- miR1-5p靶基因數(shù)量最多，并且這兩個miRNA的靶基因在5個脂代謝相關通路中都有分布，表明這兩個miRNA可能在脂代謝調控中起到核心作用。同時，部分基因也能被多個miRNA靶向，比如(長鏈酰基輔酶A合成酶4)能夠同時被ehv-miR2-3p和ehx-miR3-3p靶向。此外，病毒miRNA和宿主miRNA也能同時靶向某一代謝途徑，表明病毒可能會利用宿主miRNA來共同調節(jié)脂代謝過程。

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

富集倍數(shù)為：(某通路中差異基因數(shù)目/有注釋的差異基因總數(shù))/(某通路中的有注釋的基因總數(shù)/背景基因中有注釋的基因總數(shù))；氣泡的大小對應該通路中基因數(shù)量的多少，顏色表示富集的顯著性值。

圖5 miRNA-脂代謝靶基因互作網(wǎng)絡圖

圖中紅色圓形表示miRNA，三角形表示脂代謝相關靶基因，不同顏色對應不同脂代謝過程。

2.7 miRNA靶基因qRT-PCR驗證

為了進一步了解miRNA與脂代謝相關靶基因的靶向關系，采用qRT-PCR方法驗證了miRNA靶基因的表達水平，并分析miRNA與其靶基因表達之間的相關性。選取6個脂代謝關鍵靶基因(附表8)，包括參與脂肪酸合成的，參與鞘脂代謝過程的和，以及參與脂肪酸降解過程的,和酶基因。熒光定量結果顯示，除了基因在病毒感染早期表達水平顯著升高外，其他幾個脂代謝相關酶基因的表達水平在病毒感染早期(6 h)和晚期(45 h)均不同程度下調(附圖1)；而且這些酶基因的表達水平與miRNA的qRT-PCR結果均為負相關關系，符合一般的miRNA對靶基因的調控規(guī)律，表明這些靶基因的表達水平很可能受到miRNA的調控。

3 討論

miRNA是近年來RNA生物學領域中的重大發(fā)現(xiàn)，其通過結合在靶基因的不同位點抑制或激活靶基因的表達水平，作為豐度極高的基因調控因子參與調控多種生物學過程。特別是新近研究發(fā)現(xiàn)miRNA與它作用的靶基因是響應生物及非生物環(huán)境脅迫的主要調控因子[34]。例如，miRNA作為關鍵因子在復雜的病毒–宿主互作網(wǎng)絡中調控病毒的復制、免疫逃逸及宿主抗病毒等過程，并可能成為病毒感染診斷及治療的新手段[35]。miRNA的另一個新功能是通過控制代謝過程中的靶基因，從而調控相關代謝通路，特別是在脂質代謝網(wǎng)絡中起重要調節(jié)作用[19,36]。越來越多的研究表明，病毒感染赫氏顆石藻能夠顯著重塑宿主脂肪酸代謝、鞘脂代謝及甘油脂等脂代謝過程，且這些相關脂代謝與病毒的侵染、復制、組裝及病毒顆粒的釋放等過程密切相關[9～11]。miRNA是否參與顆石藻病毒感染重塑的脂代謝過程尚不清楚。本研究通過對EhV99B1感染的赫氏顆石藻BOF92細胞樣品進行小RNA測序，鑒定出26條新的miRNA，分析其基本特征和保守性，并結合miRNA差異表達譜特征、靶基因功能富集以及qRT-PCR實驗等，初步分析了病毒感染誘導的差異表達miRNA在宿主脂代謝過程中可能發(fā)揮的作用。

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

本研究經(jīng)測序分析共鑒定出了24條顆石藻miRNA，長度主要集中于18～21 nt，峰值為21 nt，符合一般植物miRNA的特征。BOF92的miRNA首位堿基對C有著較強的偏好，這一點與營養(yǎng)脅迫條件下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系統(tǒng)發(fā)育樹的分析結果顯示[28]，包括顆石藻在內的各種微藻的AGO蛋白都各自聚集成一類，表明這些微藻的AGO蛋白在進化上同源性較低。因此顆石藻可能也存在特殊的RNA沉默機制，使得miRNA的首位堿基不偏好于U。

在病毒誘導的BOF92 miRNA中，66.7%來自基因間區(qū)，29.2%來自外顯子區(qū)，只有1條miRNA來自內含子區(qū)(附表1)。這與近緣株系CCMP1516的miRNA來源分布較為相似，CCMP1516的miRNA約77%來自基因間區(qū)，23%來自內含子區(qū)[28]。三角褐指藻()和團藻()的miRNA也都主要來源于基因間區(qū)，只有少部分miRNA比對到了外顯子區(qū)[38,44]，這些微藻的miRNA的來源更類似于陸生高等植物[18]。萊茵衣藻()和長囊水云()的miRNA來源類似，即大部分位于內含子區(qū)，少數(shù)位于基因間區(qū)[40,39]。綜上所述，不同浮游植物miRNA的來源具有多樣性，表明miRNA在浮游植物的轉錄后水平可能具有廣泛的調控功能。

3.2 miRNA保守性特點

針對本研究中鑒定的24條宿主miRNA同源性較低的問題，推測可能有以下幾個原因。首先，不同微藻之間的miRNA可能是獨立進化的。Zhao等[40]首先對萊茵衣藻的miRNA進行了表征，作者將鑒定出來的萊茵衣藻miRNA與miRBase數(shù)據(jù)庫比對后未發(fā)現(xiàn)同源miRNA。隨后多位學者對不同種的微藻，如假微型海鏈藻()[21]，三角褐指藻()[44]及赫氏顆石藻CCMP1516[28]等在不同生長條件或環(huán)境脅迫條件下的miRNA進行了研究，也均未能在miRBase數(shù)據(jù)庫找到同源的序列。這暗示miRNA可能在浮游植物、動植物中分別有著獨立的進化機制。其次，病毒感染誘導產(chǎn)生特異性miRNA的表達。miRNA的表達是動態(tài)的，具有時空特異性，不同脅迫條件如營養(yǎng)脅迫、氧化應激、病毒感染等都會影響miRNA的表達，某些特異的miRNA可能表達水平很低，僅在特殊的生理生態(tài)條件下表達量才增加[28]。因此，病毒感染可能使赫氏顆石藻BOF92產(chǎn)生了特異性miRNA。另外，赫氏顆石藻miRNA株系間的變異較為明顯。根據(jù)本文24條宿主miRNA與其他13株赫氏顆石藻基因組的比對結果看，能夠全部比對上其他赫氏顆石藻株系的miRNA僅有10條，占miRNA總數(shù)的42%左右。赫氏顆石藻作為一種單細胞浮游植物，基因組的種內變異更有利于其適應動態(tài)變化的海洋環(huán)境以及在與病毒的“arms race”中獲得優(yōu)勢。本研究中鑒定的24條宿主miRNA保守性較低可能是這些因素共同作用的結果。此外，本研究也沒有找到EhV99B1 miRNA的同源序列。一方面，已鑒定的病毒miRNA非常有限；另一方面，與真核生物相比，病毒基因組有著更高的變異率，這意味著病毒擁有能快速適應宿主和環(huán)境條件的進化優(yōu)勢。

3.3 miRNA靶向宿主脂代謝途徑

miRNA靶基因功能富集結果顯示，病毒感染早期(6 h)和晚期(45 h)，差異表達的miRNA均能靶向宿主脂肪酸和鞘脂代謝通路(圖6)。就脂肪酸代謝而言，感染早期，宿主編碼的ehx-miR1-5p下調，其對脂肪酸合成過程的限速酶靶基因的抑制作用減弱，與此同時宿主編碼的ehx-miR4-3p上調，其對脂肪酸降解過程中的關鍵酶和靶基因的抑制作用增強，從而導致脂肪酸含量的積累(圖6A)。感染晚期，則是病毒編碼的ehv-miR2-5p 上調，并靶向抑制脂肪酸降解過程的相關酶基因(圖6B)。因此，推測感染過程中，宿主和病毒編碼的miRNA協(xié)同作用導致宿主細胞中脂肪酸含量的不斷積累。脂肪酸在病毒感染過程有著十分重要的作用，一方面大型雙鏈DNA病毒的復制和裝配有更高的代謝需求，另一方面脂肪酸可作為甘油脂和鞘脂的合成原料[7]。另外，病毒感染過程中，宿主鞘脂從頭合成過程中的兩個關鍵酶基因和的表達分別受到早期宿主上調的ehx-miR3-5p (圖6A)和晚期病毒上調的ehv-miR2-5p (圖6B)的靶向抑制，因而宿主鞘脂從頭合成途徑在整個病毒感染過程中被顯著抑制，并導致鞘脂類物質含量減少。上述結果可能在一定程度上補充解釋了本課題組之前的脂質組學數(shù)據(jù)，即EhV99B1感染早期(6 h)和晚期(45 h)，宿主赫氏顆石藻BOF92細胞中大部分脂肪酸代謝物顯著積累，而幾乎所有檢測到的16種鞘脂類物質的含量均顯著降低[11]。眾所周知，鞘脂是細胞膜上“脂筏”結構域的關鍵功能成分，也是病毒包膜結構的主要成分之一，對病毒的侵染、裝配及釋放過程十分重要，且作為信號分子，在感染晚期能夠誘導宿主程序性細胞死亡。對赫氏顆石藻CCMP1516及其特異性裂解病毒EhV86的全基因測序注釋結果顯示，通過基因水平轉移，病毒基因組“截獲”了一套鞘脂從頭合成途徑中的關鍵酶基因，并通過調節(jié)相關基因的表達在一定程度上掌控了宿主的鞘脂代謝過程[45]。在EhV201感染赫氏顆石藻CCMP374過程中發(fā)現(xiàn)，隨著感染時間的延長(感染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.

張健, 李佳芮, 楊璐, 岳娜娜. 球石藻及其生態(tài)功能. 海洋科學, 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)]

張恩權，在讀碩士研究生，專業(yè)方向：海洋微生物組學。E-mail: chriszhangen@163.com

劉靜雯，博士，教授，研究方向：海洋微型生物分子生態(tài)學。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

(責任編委: 趙呈天)