低温胁迫下海滨木槿的转录组分析

赵文静　郑旭　葛金涛　刘兴满　缪美华　孔祥伟

摘要：海滨木槿（Hibiscus hamabo Sieb.et Zucc.）是一种典型的耐盐植物，甚至能够在海水中生长，是沿海滩涂绿化的绝佳绿化树种，但在栽培试验中发现海滨木槿不耐寒，很难在北方地区越冬。为探究海滨木槿在低温胁迫下的分子应答机制，挖掘海滨木槿的抗寒基因，培育耐寒海滨木槿，通过转录组测序（RNA-seq）的试验方法分析海滨木槿在低温胁迫下的差异基因表达情况。试验结果表明，低温胁迫下的处理组与常温对照相比，共有23 689个差异表达基因，其中上调表达基因10 256个，下调表达基因13 433个。通过GO富集分析，发现差异基因大多富集在代谢过程、细胞过程、膜、催化活性等方面；通过 KEGG的富集分析，发现这些差异基因大多富集在代谢途径、次生代谢物的生物合成、光合作用、植物昼夜节律、植物激素信号转导等方面，本研究还筛选了一些潜在的耐寒差异基因，如FKF1、 GI、PYL、PP2C、CML、CPK、SNRK2、NCED和BIN2等。本试验结果可为研究海滨木槿的耐寒机制以及培育海滨木槿耐寒品种提供技术参考和基础。

关键词：海滨木槿；低温胁迫；转录组；差异基因；耐寒

中图分类号：S718文献标志码：A

文章编号：1002-1302（2023）05-0056-09

海滨木槿（Hibiscus hamabo Sieb.et Zucc.） 是锦葵科木槿属的小乔木或灌木，在我国主要分布于浙江舟山、宁波等沿海一带，在国外，日本和朝鲜也有发现报道［1-2］。海滨木槿极耐盐碱，能够在含盐量1.5%的泥滩生长，即使被海水淹浸，仍能正常生长和开花结实；同时对土壤的要求不高，在pH值为5.5～8.0之间，无论是泥滩、沙滩、丘陵地都能正常生长［3］。在以往的研究中发现，海滨木槿是一种典型的盐生植物，一种优良的园林绿化树种，能够应用到盐碱地改良和海岸滩涂绿化中，具有很高的生态价值和经济价值［4-8］。江苏地区滩涂面积较大，海滨木槿的应用前景广阔，但在引种栽培试验中，我们发现海滨木槿在江苏北部地区，如连云港地区，地上部分枝条基本无法安全越冬，只能在第2年从根部另发新枝。目前，国内外对海滨木槿的研究主要集中于驯化、形态及耐盐机制等方面［4-9］，对耐寒机制方面的研究较少。因此，基于转录组技术的海滨木槿耐寒机制的研究将为探索海滨木槿分子标记等方面提供新的思路，促进海滨木槿耐寒遗传育种的发展。

1 材料与方法

1.1 试验材料与试验地点

供试材料为1年生海滨木槿扦插盆栽苗，其母本为从浙江舟山引进，并在连云港经过5年低温驯化筛选出的耐寒品种。

试验地点为连云港市农业科学院内。

1.2 试验时间与方法

于2019年6月20日将长势一致的1年生海滨木槿扦插苗置于人工气候培养箱内，模拟低温胁迫。以20 ℃处理为对照（CK），胁迫处理（T）温度设置为4 ℃，低温处理48 h采集组织，每个处理3次重复，每组处理各取10张大小均匀的叶片，用锡箔纸包裹后迅速放入液氮中冷冻，然后放入-80 ℃超低温冰箱中保存，用于总RNA的提取。

1.3 总RNA提取

使用Omega Bio-Tek RNA提取试剂盒（货号R6827）提取RNA。利用琼脂糖凝胶电泳和NanoDrop，检验提取的RNA的完整性和纯度；最后用Agilent 2100生物分析仪确定RNA的完整性。

1.4 建库和转录组测序

建库、Unigene基本注释和序列分析方法参见赵春旭等的试验方法［10］。

1.5 差异表达基因的GO分析和KEGG分析

差异表达基因（differentially expressed genes，DEGs）分析采用DESeq2［11］进行分析，筛选阈值FDR<0.05 且|log2FC|>1。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在基因中显著性富集的Pathway。

2 结果与分析

2.1 測序数据及质量情况

通过对原始数据进行数据过滤，以减小无效数据所来带的分析干扰。对下机的 raw reads利用fastp［12］进行质控，去除不合格的reads得到clean reads。从表1可以看出，未过滤前的原始序列数据（raw datas）各处理均在3 900万以上，过滤后的数据（clean data）占比均超过99.8%，不合格的序列数据较少。从表2可以看出，Phred数值大于20的碱基占总体碱基的百分比（Q20）在97%以上，Phred数值大于30的碱基占总体碱基的百分比（Q30）在93%以上。碱基G和C的数量总和占总碱基数量的百分比（GC含量）在46%～47% 之间。说明测序的建库质量良好。

2.2 转录组的组装

采用组装序列的长度分布来衡量转录组组装质量。由表3可知，共组装了123 276个Unigene，平均长度为807 bp，N50长度为1 225 bp，最大长度为 17 146 bp，最小长度为201 bp。由图1可知，组装出来的 Unigene集中在0～1 000 bp之间，数据量及质量均能保证后续分析要求。

2.3 基因功能注释

本试验分别在NR、SwissProt、KEGG 和COG/KOG这4个数据库继续注释，结果如图2所示，4个数据库可注释的Unigene数量共76 175个，其中有34 723个Unigene能在所有数据库中得到注释。

通过NR数据库比对，结果（图3）表明，与海滨木槿同源基因相似度最高的是陆地棉（Gossypium hirsutum），相似率达到24.36%；其次是亚洲棉（G. arboreum），相似率达到21.94%；然后是雷蒙德氏棉（G. raimondil），相似率达到20.68%；与可可树（Theobroma cacao）的相似率达到18.60%。

2.4 差异基因表达分析

差异基因表达分析采用DESep2［11］，筛选条件为FDR<0.05 且 |log2FC|>1的基因为显著差异基因。如图4所示，低温胁迫下的处理组与常温对照相比，共有23 689个差异表达基因。其中上调表达基因有10 256個，下调表达基因有 13 433 个。

2.5 差异表达基因的GO分析

将差异表达基因进行GO富集功能分类，结果见图5，主要从生物过程（biological process，BP）、细胞组分（cellular component，CC）、分子功能（molecular function，MF）3个方面进行分析。

在生物学过程中，试验材料DEGs 富集最多的依次为代谢过程（metabolic process），有6 595 个Unigene得到注释，占比81.3%；细胞过程（cellular process），有6 019 个Unigene得到注释，占比 74.2%；单有机体过程（single-organism process），有4 653 个Unigene得到注释，占比 57.36%；生物调节（biological regulation），有2 493 个Unigene得到注释， 占比30.73%； 生物过程的调节（regulation of biological process），有2 284个Unigene得到注释，占比 28.16%；对刺激的反应（response to stimulus），有1 954 个Unigene得到注释，占比24.09%；本土化（localization），有1 556个Unigene得到注释，占比19.18%；发展过程（developmental process），有1 089个Unigene得到注释，占比13.42%；细胞成分组织（cellular component organization or biogenesis），有 1 197 个Unigene得到注释，占比14.76%；信号（signaling），有712 个Unigene得到注释，占比8.78%； 多细胞生物的过程（multicellular organismal process），有834 个Unigene得到注释，占比10.28%；多有机体过程（multi-organism process），有390个Unigene得到注释，占比4.81%；生殖（reproduction），有516个Unigene得到注释，占比6.36%；生殖过程（reproductive process），有511个Unigene得到注释，占比6.3%；免疫系统过程（immune system process），有175个Unigene得到注释，占比2.16%；生物过程的正向调控（positive regulation of biological process），有203个Unigene得到注释，占比2.5%；生长（growth），有159个Unigene得到注释，占比1.96%；生物过程的负调控（negative regulation of biological process），有189个Unigene得到注释，占比2.33%；节奏过程（rhythmic process），有62 个Unigene得到注释，占比0.76%；生物黏附（biological adhesion），有40个Unigene得到注释，占比0.49%；排毒（detoxification），有10 个Unigene得到注释，占比0.12%；运动（locomotion），有5个Unigene得到注释，占比0.06%；细胞死亡（cell killing），有1个Unigene得到注释，占比0.01%。

在细胞组分中，试验材料DEGs 富集最多的依次为细胞（cell和cell part），有3 832 个Unigene得到注释，占比74.99%；细胞器（organelle），有3 173个Unigene得到注释，占比 62.09%；膜（membrane），有2 142个Unigene得到注释，占比 41.92%；膜部分（membrane part），有1 536个Unigene得到注释，占比 30.06%；细胞器部分（organelle part），有1 392 个Unigene得到注释，占比27.24%；大分子复合物（macromolecular complex），有742个Unigene得到注释，占比14.52%；细胞连接（cell junction），有228 个Unigene得到注释，占比4.46%；细胞外区（extracellular region），有134 个Unigene得到注释，占比2.62%；膜封闭腔（membrane-enclosed lumen），有74 个Unigene得到注释，占比1.45%；病毒粒子（virion），有20个Unigene得到注释，占比0.39%；类核素（nucleoid），有8 个Unigene得到注释，占比0.16%；细胞外区部分（extracellular region part），有4个Unigene得到注释，占比0.08%；细胞外基质（extracellular matrix），有4 个Unigene得到注释，占比0.08%；超分子纤维（supramolecular fiber），有3个Unigene得到注释，占比0.06%。

在分子功能中，试验材料DEGs 富集最多的依次为催化活性（catalytic activity），有5 259 个Unigene得到注释，占比68.86%；黏合物（binding），有4 632个Unigene得到注释，占比60.65%；转运蛋白活性（transporter activity），有454个Unigene得到注释，占比5.94%；核酸结合转录因子活性（nucleic acid binding transcription factor activity），有382个Unigene得到注释，占比5%；信号传感器活动（signal transducer activity），有119个Unigene得到注释，占比1.56%；结构分子活性（structural molecule activity），有85个Unigene得到注释，占比1.11%；分子换能器活性（molecular transducer activity），有81 个Unigene得到注释，占比1.06%；分子功能调节剂（molecular function regulator），有42 个Unigene得到注释，占比0.55%；抗氧化活性（antioxidant activity），有18个Unigene得到注释，占比0.24%；转录因子活性，蛋白质结合（transcription factor activity，protein binding），有14个Unigene得到注释，占比0.18%；电子载流子活性（electron carrier activity），有10个Unigene得到注释，占比0.13%；翻译调节活性（translation regulator activity），有1个Unigene得到注释，占比0.01%。

其中，按富集程度（qvalue＜0.05）由高到底排序依次为质体（plastid）、质体部分（plastid part）、类囊体（thylakoid）、质体类囊体（plastid thylakoid）、叶绿体（chloroplast）、葉绿体部分（chloroplast part）、质体基质（plastid stroma）、质体被膜（plastid envelope）等（图6）。

2.6 差异表达基因的KEGG富集分析

生物个体的功能表达通常需要大量基因的共同协作才能完成，因此通过对代谢通路的分析，能够更加清晰地表达 DEGs 的生物学功能。将差异表达基因进行KEGG分析，并在130个代谢通路中得到注释。富集程度（qvalue＜0.05）最高的20个通路分别是：代谢途径（metabolic pathways）、次生代谢物的生物合成（biosynthesis of secondary metabolites）、光合作用（photosynthesis）、植物的昼夜节律（circadian rhythm-plant）、卟啉和叶绿素代谢（porphyrin and chlorophyll metabolism）、光合作用-天线蛋白（photosynthesis-antenna proteins）、甘氨酸、丝氨酸和苏氨酸代谢（glycine，serine and threonine metabolism）、光合生物中的碳固定（carbon fixation in photosynthetic organisms）、二萜生物合成（diterpenoid biosynthesis）、类胡萝卜素生物合成（carotenoid biosynthesis）、碳代谢（carbon metabolism）、植物-病原体相互作用（plant-pathogen interaction）、半胱氨酸和蛋氨酸代谢（cysteine and methionine metabolism）、氮代谢（nitrogen metabolism）、植物激素信号转导（plant hormone signal transduction）等。表示富集程度的qvalue均＜0.001，其中得到注释的 DEGs数量最多的代谢通路是代谢途径，共有1 752个。其次，次生代谢物的生物合成有995个DEGs，植物的昼夜节律有121个DEGs，甘氨酸、丝氨酸和苏氨酸代谢有100个DEGs（图7）。

2.7 抗寒差异表达基因

本试验在对差异表达基因功能分析的基础上，对试验材料特异表达的基因进行进一步的分析，挑选出一些与植物耐寒相关的基因，发现这些基因主要分布在代谢途径、次生代谢物的生物合成、卟啉和叶绿素代谢、氨基酸代谢、类胡萝卜素生物合成、光合作用与碳代谢、昼夜节律、植物激素信号转导等通路中。

如表4所示，与代谢途径相关的表达基因有E1.6.5.4、NCED、E2.4.1.13；与次生代谢物的生物合成相关的表达基因有katE，CAT，catB，srpA、G6PD，zwf；与卟啉和叶绿素代谢相关的表达基因有hemD，UROS、hemA；与氨基酸代谢相关的表达基因有AGXT2、serC，PSAT1；与类胡萝卜素生物合成相关的表达基因有CYP707A、AAO3、ZEP，ABA1；与光合用与碳代谢相关的表达基因有petH、psaG、LHCB6、MDH2、IDH3、PGD，gnd，gntZ；与植物昼夜节律相关的表达基因有GI、FKF1；与植物激素信号转导相关的表达基因有CML、CPK、SNRK2、BIN2、IAA、PP2C、PYL。

3 讨论与结论

温度是植物生长发育过程中最受环境制约的因素之一，限制了植物的生长发育和地理分布［13］，以至于有了“江南冬景似春华，江北冬岭只枯丫”的对比。植物在受到低温胁迫的时候，细胞膜系统首先遭到破坏，膜透性改变，继而扰乱细胞内正常代谢，如光合作用、呼吸作用等，严重时导致细胞凋亡［14-19］。为了研究植物体内应答低温胁迫的调控机制，随着近年来技术的发展，RNA-seq、Chip-seq、iTRAQ、TMT、LC-MS/MS 等新一代技术应运而生，分别从转录组、蛋白组、代谢组及转录翻译调控等不同方向对植物的抗寒机制进行解析［20］。其中，转录组测序（RNA-seq）是对某一物种的mRNA进行高通量测序的技术，反映在特定条件和特定时间下该物种的基因表达情况［21-22］，因此，通过转录组测序，有利于揭示作物在冷胁迫下的生理及分子机制。目前为止，转录组测序手段已经揭示了小麦、水稻、茶树、薰衣草、草地早熟禾、烟草、拟南芥、西葫芦等经济作物的耐寒分子机制，为培育其耐寒品种提供了信息技术参考［10，13，15，23-24］。

本研究利用转录组测序方法，研究了经过低温驯化的耐寒海滨木槿在低温胁迫处理下的基因表达情况。经过转录组测序、数据组装和注释后，通过NR数据库比对与海滨木槿同源基因相似度最高的是陆地棉（G. hirsutum），相似率达到24.36%；其次是亚洲棉（G. arboreum），相似率达到21.94%；然后是雷蒙德氏棉（G. raimondil）的相似率达到20.68%；可可树（Theobroma cacao），相似率达到18.60%。其中，陆地棉、亚洲棉、雷蒙德氏棉均为锦葵科木槿族棉属，与海滨木槿同为锦葵科植物，可可树也是锦葵目植物，说明测序准确性是比较高的。采用DESep2进行差异基因表达分析，筛选出FDR<0.05 且 |log2FC|>1 的基因为显著差异基因，共得到23 689个差异表达基因，其中上调表达基因有10 256个，下调表达基因有13 433个。通过GO富集分析，发现DEGs主要富集在代谢过程、细胞过程、膜、催化活性、黏合物等方面。这与相关报道较为一致，植物在遭遇低温胁迫后，首先细胞膜透性会改变，传递信号，调控细胞代谢，以应对低温带来的伤害［25-26］。通过 KEGG的富集分析，发现差异基因主要富集在代谢途径、次生代谢物的生物合成、光合作用、植物昼夜节律、卟啉和叶绿素代谢、甘氨酸、丝氨酸和苏氨酸代谢、二萜生物合成、类胡萝卜素生物合成、碳代谢、植物-病原体相互作用、半胱氨酸和蛋氨酸代谢、氮代谢、植物激素信号转导等，说明这些途径可能在海滨木槿低温胁迫响应中发挥作用。

通過进一步的发掘，发现耐寒相关的差异基因主要分布在代谢途径、次生代谢物的生物合成、卟啉和叶绿素代谢、氨基酸代谢、类胡萝卜素生物合成、光合作用与碳代谢、昼夜节律、植物激素信号转导等通路中。梁雨晨等认为，在薰衣草昼夜节律通路中，植物受光敏色素（PHY）调控，诱导产生相关抗性基因，通过蛋白网络互作分析和筛选，发现枢纽基因 FKF1和GI，在昼夜节律通路中上调表达，参与调控植物昼夜节律信号通路，表明生物钟调控在薰衣草应答寒冷胁迫中发挥了重要作用［24］。本研究中，筛选出昼夜节律中的差异基因GI和FKF1均为上调表达，与梁雨晨的研究结果［24］一致。Chen等通过对经过低温胁迫处理后耐寒基因型菠萝的转录组测序分析，DEGs 通过 KEGG 分析，发现与细胞壁特性、气孔关闭以及脱落酸（ABA）和活性氧（ROS）信号转导相关的基因在菠萝耐寒性中起重要作用，其中，所有与ABA信号转导相关的冷胁迫响应基因表达均上调，包括PYR/PYL同源基因和PP2C同源基因［27］。本研究与之相一致，在本研究中，筛选出的差异基因PYL和PP2C均上调表达。赵春旭等在探究青海野生草地早熟禾低温胁迫下的分子应答机制时，发现在低温胁迫下，差异基因主要富集在光合作用、氧化还原反应过程、碳水化合物代谢、细胞膜系统、转运蛋白以及次生物代谢中，另外通过耐寒和不耐寒2种材料转录组的对比，筛选出一批在耐寒材料中上调表达的潜在的抗寒基因，如CML、CPK、CALM、DHAR、GST、NCED、SNRK2、BSK、CKX、BIN2、ARF和PEK等［11］。本研究的结果与之一致，在本研究中，筛选出的CML、CPK、SNRK2、NCED、BIN2等与赵春旭研究中重合的基因，也均上调表达。

本研究利用高通量转录组测序技术建立了海滨木槿低温胁迫前后叶片的转录组数据库，筛选出差异表达基因，为研究海滨木槿的耐寒机制以及培育海滨木槿耐寒品种提供了技术参考和基础。

参考文献：

［1］张朝芳. 浙江植物志：第一卷 蕨类植物·裸子植物［M］. 杭州：浙江科学技术出版社，1993.

［2］张若蕙. 浙江珍稀濒危植物［M］. 杭州：浙江科学技术出版社，1994：330-339.

［3］唐丽红，马明睿，韩 华，等. 上海市景观水体水生植物现状及配置评价［J］. 生态学杂志，2013，32（3）：563-570.

［4］肖祖飞，赵 姣，张 杰，等. 海滨木槿繁殖技术［J］. 北方园艺，2017（12）：74-77.

［5］於朝广，殷云龙，芦治国，等. 海滨木槿研究进展及应用前景［J］. 现代农业科技，2016（19）：140-143，145.

［6］李翠云，姜彦成，乔桂荣，等. 海滨木槿叶片蛋白质双向电泳体系的建立［J］. 西南林业大学学报，2012，32（4）：30-35.

［7］王宗星. 海滨木槿和蜡杨梅对海水胁迫的光合生理响应［D］. 北京：中国林业科学研究院，2011.

［8］王连吉. 新优树种海滨木槿在风景园林中的应用［J］. 中国园林，2010，26（4）：49-50.

［9］王秀丽，张 荻，刘红梅，等. 海滨木槿耐盐性的初步研究［J］. 上海交通大学学报（农业科学版），2010，28（3）：248-254.

［10］赵春旭，马 祥，董文科，等. 低温胁迫下不同青海野生草地早熟禾的转录组比较分析［J］. 草地学报，2020，28（2）：305-318.

［11］Love M I，Huber W，Anders S.Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2［J］. Genome Biology，2014，15（12）：550.

［12］Chen S F，Zhou Y Q，Chen Y R，et al. Fastp：an ultra-fast all-in-one FASTQ preprocessor［J］. Bioinformatics，2018，34（17）：1884-1890.

［13］张 健，唐 露，冉启凡，等. 植物响应低温胁迫转录组测序研究进展［J］. 分子植物育种，2020，18（6）：1849-1866.

［14］Wang H S，Yu C，Zhu Z J，et al. Overexpression in tobacco of a tomato GMPase gene improves tolerance to both low and high temperature stress by enhancing antioxidation capacity［J］. Plant Cell Reports，2011，30（6）：1029-1040.

［15］李 楠，郑勇奇，丁红梅，等. 低温胁迫下短枝木麻黄耐寒相关基因的差异表达分析［J］. 林业科学，2017，53（7）：62-71.

［16］刘文盈. 低温胁迫对黑果枸杞愈伤组织的影响［J］. 基因组学与应用生物学，2018，37（1）：408-412.

［17］Colinet H，Renault D，Roussel D. Cold acclimation allows Drosophila flies to maintain mitochondrial functioning under cold stress［J］. Insect Biochemistry and Molecular Biology，2017，80：52-60.

［18］Walker D J，Romero P，Correal E. Cold tolerance，water relations and accumulation of osmolytes in Bituminaria bituminosa［J］. Biologia Plantarum，2010，54（2）：293-298.

［19］宋丽莉. 紫花苜蓿抗寒转录组测序分析及MsNCX和MsCML基因的功能研究［D］. 哈尔滨：哈尔滨师范大学，2017.

［20］张 纯，唐承晨，王吉永，等. 转录组学在植物应答逆境胁迫中的研究进展［J］. 生物学杂志，2017，34（2）：86-90.

［21］慧 芳，刘秀岩，李宗谕，等. 转录组测序技术在药用植物研究中的应用［J］. 中草药，2019，50（24）：6149-6155.

［22］王丽鸳. 基于EST数据库和转录组测序的茶树DNA分子标记开发与应用研究［D］. 北京：中国农业科学院，2011.

［23］張宏亮. 温胁迫下西葫芦转录组分析与SSR分子标记开发［D］. 太谷：山西农业大学，2015.

［24］梁雨晨，蔺海娇，侯一航，等. 狭叶薰衣草叶片转录组耐寒相关差异表达基因分析［J/OL］. 分子植物育种：1-14［2022-11-10］. http：//kns.cnki.net/kcms/detail/46.1068.S.20210604.1045.010.html.

［25］Orvar B L，Sangwan V，Omann F，et al. Early steps in cold sensing by plant cells：the role of actin cytoskeleton and membrane fluidity［J］. The Plant Journal，2000，23（6）：785-794.

［26］范吉标. 狗牙根耐寒生理及分子机制解析［D］. 武汉：中国科学院研究生院（武汉植物园），2016.

［27］Chen C J，Zhang Y F，Xu Z Q，et al. Transcriptome profiling of the pineapple under low temperature to facilitate its breeding for cold tolerance［J］. PLoS One，2016，11（9）：e0163315.

收稿日期：2022-05-16

基金项目：江苏省林业科技创新与推广项目（编号：LYKJ-连云港［2021］01）；连云港市财政专项（编号：QNJJ1921）。

作者简介：赵文静（1991—），女，江苏连云港人，硕士，助理研究员，主要从事植物病理学和植物育种研究。E-mail：ZWJ6691@163.com。

通信作者：刘兴满，硕士，副研究员，主要从事园艺和植物育种。E-mail：13611551930@126.com。