不同光照强度胁迫下的多型杜氏藻转录组分析

高 帆， 杨 鹏， 王 飞， 徐 璐， 汤丽群

(1. 山西大学 生命科学学院， 太原 030006； 2. 山西大学教务处， 太原 030006；3. 山西大学生物技术研究所， 太原 030006； 4. 山西大学 体育学院， 太原 030006)

光照胁迫是植物界一类重要的非生物胁迫，对植物的光合作用、生长发育、呼吸作用具有重要的调节作用[1]。目前，光照胁迫的相关研究主要集中在高等陆生植物，尤其是农作物，对低等藻类植物的相关研究较少。从藻类中挖掘与光照胁迫响应紧密相关的调控因子及关键基因，并对其区别于高等植物的独特调控方式进行分析，对藻类植物的适应性进化、生长发育及光合作用调控机制的研究很必要。

杜氏藻属绿藻门，为单细胞真核微藻，无细胞壁，具鞭毛可游动，胞内含较大的杯状叶绿体能进行光合作用；可在含有0.05～5 mol/L NaCl的盐水中生长，抗逆性极佳，是一种较为理想的抗逆境生物。杜氏藻中富含生物活性物质和藻类色素等高附加值产物，如多糖、甘油、多不饱和脂肪酸、类胡萝卜素等[3]，在医疗保健、食品加工、生物柴油等领域具有广阔的应用前景，是一类重要的特色经济微藻。多型杜氏藻(Dunaliellapolymorpha)是杜氏藻属的一种[2]，目前该藻株多在内陆盐湖、近海海域发现。

随着国内外高通量测序技术的快速发展，先后出现了以Ion Torrent、Roche 454、Thermo Fisher、SOLiD、HiSeq、BGISEQ为代表的高通量测序平台[4]。对杜氏藻株进行转录组测序，已成为获取其在不同胁迫处理条件下关键基因、调控因子的重要手段之一。近年来，研究人员分别在高盐、重金属等非生物胁迫下对不同杜氏藻株进行了转录组测序，获取抗逆响应相关基因[8-10]。然而，上述研究均基于杜氏藻属的其他种类，对多型杜氏藻，目前仅限于其形态学鉴定及系统分类方面[2]，对该藻株光照胁迫下关键应答基因的筛选、光合调控机理机制的相关研究仍未有报道。光照是影响杜氏藻室内培育、生物活性物质富集、代谢产物提取重要的影响因素[5-7]。研究在细胞、生化、转录组水平上，对不同光强胁迫下多型杜氏藻株的形态、光合效率、基因表达及代谢通路进行分析，填补了该藻株光合作用机理机制研究的空白，研究结果为多型杜氏藻不同光强胁迫下分子应答机制的深入研究，高品质优势藻株的遗传改良奠定了理论基础，同时也为其光胁迫响应紧密相关代谢产物的工业应用提供了借鉴。

1 材料与方法

1.1 材料搜集与培养

多型杜氏藻种(CCAP 19/14)购于英国CCAP藻种保藏中心，用优化的D.m培养基配方[11]进行光照静态培养，温度为25 ℃±5 ℃，光周期为16 h∶8 h，培养室无菌通风。起始光照强度设置为5 000 lx，光源采用白炽灯，每隔10 d采样进行显微镜检，确保无染菌。根据实验室前期对杜氏藻生长周期的测定，第60天进入成熟期，此时将培养藻株平均分为3份，一份进行暗室遮光处理(0 lx)，另两份分别进行中强度(10 000 lx)和高强度(20 000 lx)光照胁迫处理，处理时间均为30 min，遮光处理作为对照组(DP_CO)、中强度光照和高强度光照胁迫分别作为实验组DP_ML和DP_HL。

1.2 形态学观察与叶绿素荧光参数测定

利用D-SW50光学显微镜(光学仪器一厂，中国上海)观察处理后的多型杜氏藻单细胞显微形态。利用FMS-2叶绿素荧光测定仪(Hansatech公司，美国诺福克)对处理后的多型杜氏藻液的PSII原初光能转化效率(Fv/Fm)和潜在光合活性(Fv/Fo)进行测定，测定方法参照仪器使用说明及文献[12]的方法进行。

1.3 总RNA提取、文库构建及测序

分别取处理后的样本300 mL，每份样本再平均分为3份(进行3个生物学重复)，低温离心后用液氮速冻并快速研磨，按InvitrogenTRIzol试剂盒[赛默飞世尔科技(中国)有限公司，中国上海]说明分别提取藻株的总RNA，Agilent 2100 Bioanalyzer[安捷伦生物(杭州)有限公司，中国杭州]对其质量和浓度进行评估，质量检测合格(OD260/OD280比值在1.8～2.0，RIN值>7.0)的RNA样本[13]用Promega反转录试剂盒[普洛麦格(北京)生物技术有限公司，中国北京]进行反转录(操作步骤按照说明书进行)，按照SMART cDNA文库构建试剂盒[宝生物工程(大连)有限公司，中国大连]说明书分别对样本进行文库构建。构建好的文库在BGISEQ-500测序平台(华大基因，中国深圳)进行转录组测序。

1.4 测序质量分析及片段组装

利用SOAPnuke2.0软件[14]删除低质量reads，即测序片段5′和3′端含接头，测序片段上poly结构中N占比大于5%(N表示未确定碱基)，筛选获得的clean reads(Q20值>95%，Q30值>85%)，利用Trinity2.4.0软件对clean reads进行denovo从头组装。利用Bowtie2.0软件[15]预测clean reads中的Unigene，利用Tgicl 2.1软件[16]分析Unigene的数量和长度，利用BUSCO 4.0.2软件[17]评估组装质量。

1.5 数据注释

将Unigene与Gene ontology(GO)、RefSeq nonredundant protein(NR)、Nucleotide sequence(NT)、SwissPort、Pfam、Kyoto encyclopedia of genes and genomes(KEGG)及euKaryotic orthologous group(KOG)数据库分别进行数据比对，参数设置：E值≤1.00E-15，查询序列长度百分比≥90%。比对获得的数据在Linux系统下运行uniq命令，对其中的重复序列进行去冗余处理。去冗余后的数据在UNIX系统下执行HMMER程序，查询Pfam和COG数据库，预测Unigene功能及潜在的转录因子[18]。

1.6 基因表达分析

Fragments per kilobase of transcript per million fragments(FPKM)计算Unigene表达量[19]，RSEM软件筛选差异表达基因(Differentially expressed gene，DEG)，DEseq2软件构建DEGs分布火山图(Fold change≥|5.00|，Q值≤ 0.05)[20]，显著差异基因采用Fold change≥|25|，Q值≤0.01的标准进行筛选；相同标准下，用Limma/voom软件再次统计DEGs数量，比较分析两者共有DEGs筛选结果及其基因表达相关性；利用层次聚类法构建不同处理组间的DEGs聚类热图[21]。

1.7 GO和KEGG通路富集

运行Blast2GO对Unigene进行GO分析，并对GO条目进行富集(FDR≤0.01)。运行KAAS程序对Unigene进行KEGG通路分析，并对KEGG通路进行富集(FDR ≤0.01)[22]。

2 结果与分析

2.1 藻细胞形态与PSII光化学效率分析

由图1(a)可知，遮光处理后的成熟期藻细胞呈球状，浅绿色，有鞭毛，细胞直径为(16±2)μm；中强度光照处理后的藻细胞呈椭球状，深绿色，有鞭毛，细胞直径为(23±1.5)μm；高强度光照处理后的藻细胞呈椭球状，黄绿色，有鞭毛，细胞直径为(24±2)μm。表明随着光刺激的加强，藻细胞色素尤其是叶绿素的合成加强，细胞形态也发生了变化。由图1(b)可知，遮光处理后的多型杜氏藻光合效率较低，Fv/Fo和Fv/Fm的变化趋势均是先降低后逐渐升高，这表明遮光处理极大地影响藻株的光化学效率，但随着光照的恢复，光化学效率逐渐增强。中光强处理后的多型杜氏藻光合效率最高，Fv/Fo和Fv/Fm的变化趋势均是先降低后升高，后趋于平稳，表明中光强(10 000 lx)适合该类型藻株进行光合作用，当光照强度降低后，光化学效率也降低，但随着藻株的生理适应，效率又逐渐恢复。高光强处理后的多型杜氏藻光合效率低于中光强组，Fv/Fo的变化趋势是先降低后升高，再降低，而后趋于平稳；Fv/Fm的变化趋势是先降低后升高，这表明强光(20 000 lx)不利于该藻株进行光合作用，影响了其光化学效率的提高，但随着光强的降低及藻株的生理适应，其光化学效率虽得以恢复，但仍然无法达到原始水平。

2.2 转录组测序及组装分析

多型杜氏藻总RNA平均浓度为95 ng/μL，OD260/OD280在1.9～2.0(平均值为1.93)，RIN平均值为7.2，符合转录组测序标准。由表1和图 2(a)可知，3个实验处理组DP_CO(0 lx)、DP_ML(10 000 lx)、DP_HL(20 000 lx)共9个测序文库中共鉴定到52 853条Unigenes，共检测到282条完整的Contig片段(完整度占比约92.01%，contig总长约4.25×107bp，平均长度约804 bp)。N50长1 072 bp(N70长785 bp，N90长396 bp)，GC占比约52.25%。Unigene中共预测到29 025条CDS(总长为2.05×107bp)，14 220条SSR(总长为18 643 bp)。由图2(b)可知，9个测序样本的Q20和Q30百分比平均值分别为98.95%和89.72%。Denovo组装的BUSCO质量评估结果显示，组装完整的测序片段、单拷贝及双拷贝片段总占比大于90%，组装质量较高。

(a)藻细胞光胁迫处理30 min后的形态观察(放大倍数20×)；(b)藻液光胁迫处理30 min后的Fv/Fo及Fv/Fm比值。图1 多型杜氏藻形态学观察及PSII潜在光合效率及原初光能转化率Figure 1 Morphology of D. polymorpha and potential PSII acitivityand its primary light energy conversion efficiency

表1 多型杜氏藻转录组测序de novo组装基本信息Table 1 D. polymorpha transcriptome assembly information

(a)Unigene长度分布；(b)基于BUSCOs的de novo组装质量评估。图2 转录组测序Unigene长度及组装质量分析Figure 2 Analysis of transcriptomic unigenes length and de novo assembly quality

2.3 基因注释与功能预测

52 853条Unigenes在7大生物信息数据库中的交集注释数为21 980，占比41.59%(并集注释数为30 873，占比68.42%)，见表2，平均49.71%的基因在上述数据库中有注释信息。基于NR数据库的多型杜氏藻Unigne注释物种的分类结果显示[图3(a)]，约35.33%的Unigene可在藻类(Chlamydomonaseustigma，Goniumpectorale和Volvoxcarteri)和陆生植物(Dorcocerashygrometricum和Ricinuscommunis)中被注释到，剩余64.67%的数据未在已知植物物种基因库中检测到。由图 3(b)可知，COG分类中最多的10 377条Unigene功能与催化活性有关(占比20.22%)，其次为生物学调控功能(9 339条Unigene，占比18.20%)。由图3(c)可知，基因所属转录因子家族分类结果中预测最多的转录因子家族是MYB(31条Unigene，占比10.58%)。

(a)多型杜氏藻Unigene物种分类注释；(b)多型杜氏藻Unigene的COG分类；(c)多型杜氏藻Unigene的转录因子家族分类。图3 多型杜氏藻Unigene注释Figure 3 Unigene annotation of D. polymorpha

表2 多型杜氏藻转录组Unigene注释信息

2.4 基因表达分析

基于多型杜氏藻3个处理组Unigene表达的FPKM值，构建基因表达聚类图谱[图4(a)]。多型杜氏藻中光强处理组(DP_ML)和高光强处理组(DP_HL)的基因集可聚为一类，表明两者具有相似的基因表达趋势，遮光处理组(DP_CO)的基因集单独聚为一类，所有Unigene可划分为两大簇。

以DP_CO为对照，共检测到18 053条DEGs，其中显著差异表达基因1 005条。如图4(b)和(c)所示，DP_ML vs DP_CO共检测到上调基因5 825条，下调基因6 250条，非差异表达基因212条；DP_HL vs DP_CO共检测到上调基因2 221条，下调基因3 135条，非差异表达基因410条。如图4(d)所示，两对比较组间的共有DEGs 3 123条。如图 4(e)和(f)所示，DEseq2分析获得的差异基因结果稳定(Limma/voom与DEseq2分析获得的共有DEGs达17 849条，占比分别为98.86%和99.01%；两款软件获得的DEGs表达水平相关性较高，R2=0.987 6)。

(a)基因表达聚类热图；(b)DP_ML vs DP_CO DEGs火山图；(c)DP_HL vs DP_CO DEGs火山图；(d)两对比较组间共有DEGs维恩图；(e)Limma/voom与DEseq2分析所得共有DEGs维恩图；(f)Limma/voom与DEseq2分析所得共有DEGs表达相关性。图4 多型杜氏藻不同比较组间的DEGs分析Figure 4 Analysis of differentially expressed genes among different D. polymorpha comparison groups

2.5 GO与KEGG通路分析

多型杜氏藻组间共有DEGsGO分类结果，见图 5(a)，在生物学过程类别中，GO条目最多的是细胞过程(235条)；在细胞组分类别中，GO条目最多的是细胞(208条)；在分子功能类别中，GO条目最多的均是催化活性(320条)。图 5(b)所示组间共有DEGsGO富集top20结果，富集程度最高的GO条目是光合作用-天线蛋白(基因数量：119条，Q值：1E-120)。

多型杜氏藻组间共有DEGsKEGG通路分类结果，见图 5(c)，所有的通路可划分为细胞过程、环境信息处理、遗传信息处理、新陈代谢和有机体系统等5类。其中，通路条目最多的是新陈代谢(148条)。组间共有DEGsKEGG通路富集top20结果，见图 5(d)，富集程度最高的KEGG通路是光胁迫应答(基因数量：38条，Q值：1E-40)。

(a)组间共有DEGs GO分类；(b)组间共有DEGs GO富集；(c)组间共有DEGs KEGG通路分类；(d)组间共有DEGs KEGG通路富集。图5 多型杜氏藻组间的差异表达基因GO及KEGG分类与富集Figure 5 GO and KEGG classification and enrichment of the common DEGs among different D. polymorpha comparison groups

3 讨论

利用BGISEQ-500测序平台[21]对不同光强胁迫下的多型杜氏藻进行转录组测序，测序结果可靠。单藻株clean reads的GC百分比(51.93%)高于Panahi等[22]利用Illumina GAIIx测序平台获得的盐生杜氏藻(D.salina)数据(clean reads GC%：50%)。多型杜氏藻单藻株测序文库质量指标Q20和Q30平均值高于Polle等[23]报道的盐生杜氏藻CCAP 19/18转录组测序质量。多型杜氏藻中预测的CDS数(29 025条)低于He等[10]在Illumina HiSeq 2000测序平台获得的盐生杜氏藻CCAP 19/30CDS数(49 751条)。多型杜氏藻转录片段从头组装后的N50长度(1 072 bp)短于Hong等[13]从盐生杜氏藻435所得N50长度(1 727 bp)。

多型杜氏藻Unigene中仅有35.33%可在已报道植物基因库中注释到，相较于其他藻类植物如莱茵衣藻，基因的物种注释信息量较低(一般高于45%)[24]。推测有两方面原因：(1)绿藻的系统进化时间早于高等陆生植物，晚于蓝细菌，甚至与某些原生动物系统进化关系较近[25]；未将该藻基因的物种注释范围扩大到细菌及原生动物，导致可比对数据范围缩小。(2)研究测序获得的多型杜氏藻Unigene数据未经开放型阅读框整合处理，造成非功能测序片段对物种注释结果的干扰。

多型杜氏藻PSII光合效率指标显示，尽管中高强度光照胁迫后的藻株Fv/Fo及Fv/Fm均有波动，但中强光胁迫组的光合效率在短期内可实现平稳增加，转录组测序结果显示DP_ML vs DP_CO的DEGs数量较高(12 075条)，这暗示随着光强的增加，藻株光合作用也逐渐增强，相关功能基因表达出现波动。然而，多型杜氏藻高强度光照胁迫组的PSII光合效率在短期内无法达到原初水平，转录组测序结果显示，DP_HL vs DP_CO的DEGs数量较低(5 356条)，这暗示当光强达到一定阈值(～20 000 lx)，藻株可能出现了代谢废物与有毒物质积累，甚至细胞凋亡，为调控外界刺激的影响，DEGs数明显下降，部分基因沉默，这可能是藻株为应答不利环境而做出的响应。此外，相对遮光处理组，多型杜氏藻中高光强胁迫组下调基因数均高于上调基因数，表明对该藻株光合作用、生长发育及新陈代谢等过程起关键作用的光胁迫响应基因，其总体调控方式是负向的。该结果与杜氏藻在其他非生物胁迫下(如高盐胁迫、高温胁迫、强酸强碱等)的响应基因调控方式不同(一般为正向调控)[10,26]。

多型杜氏藻组间共有DEGsGO分析表明，与光照胁迫响应相关基因功能可能主要与藻细胞内与光合作用相关功能组分的信号识别、催化过程有关(GO富集分类结果最高的分别是细胞过程、细胞和催化活性，GO富集程度最高的是光合作用-天线蛋白)，光胁迫下藻细胞的形态变化(球形延展为椭球形，叶绿素含量增加)也暗示了上述预测结果[27]。多型杜氏藻组间共有DEGsKEGG通路分析表明，与光胁迫响应相关基因的功能可能主要参与了与光合调控相关的代谢通路(KEGG分类结果最高的是新陈代谢，通路富集程度最高的是光胁迫应答)。

4 结论

多型杜氏藻在受到外界环境变化刺激后，会启动自身内在调控机制来适应新环境。研究基于多型杜氏藻转录组测序结果，分析其在不同光照强度胁迫下的差异表达基因及其代谢通路，结合其形态与光合效率指标，进一步分析该藻株潜在的光胁迫应答过程。研究结果为该藻株光胁迫响应过程中关键调控基因的挖掘与功能验证及光合调控分子机制的研究奠定基础，同时也为特色藻株的分子改良及其工业应用提供借鉴。