广西蔗糖产量约占全国食糖产量的60%，蔗糖产业是广西的传统优势支柱产业^{[1, 2]}。在经常出现冬季持续低温及春季“寒潮”带来的低温阴雨天气下，广西现有的甘蔗(Saccharum officinarum L.)主栽品种新台糖22号(ROC22)出现叶片萎蔫、生长点损伤、生长迟缓、糖分降低等情况^{[3, 4]}。因此，选育耐低温甘蔗品种是广西乃至我国目前甘蔗育种关注的重点之一。

甘蔗对低温的耐受性往往是多基因控制的复杂性状，响应低温胁迫的基因通常涉及低温信号感知、传导与调控、抗氧化酶系统等方面。在甘蔗响应低温胁迫的研究中，Menossi等^[5]将甘蔗幼苗在4 ℃胁迫处理后，分离鉴定得到25个表达受到抑制的基因和34个表达上调的基因，这些基因与干旱胁迫诱导产生的差异基因相同。陈香玲等^[6]利用起始密码子多态性(Start Codon Targeted Polymorphism, SCoT)技术分析了甘蔗在低温胁迫下稳定出现的可能与甘蔗的抗寒基因表达相关的5个cDNA差异片段。黄杏等^[7]应用同源克隆和逆转录PCR (Reverse Tranion PCR, RT-PCR)技术克隆了5个与甘蔗抗寒性相关的基因的cDNA，并发现这些基因被脱落酸(Abscisic Acid, ABA)和低温协同调控，ABA预处理不同程度地诱导其表达。Khan等^[8]通过抑制性消减杂交技术，从抗寒甘蔗品种中分离出465个低温响应基因，这些低温响应基因通过几个相关途径协同调节其交互作用，从而实现甘蔗的抗寒应答反应。成伟等^{[9, 10]}利用RT-PCR技术成功克隆甘蔗Sc CBF1基因，并发现该基因在低温胁迫下表达迅速上调，推测该基因在逆境信号传导及植物抗逆过程中扮演重要角色。甘蔗是高度异源多倍体，遗传背景复杂，其对低温的应答机制研究还相对滞后。转录组测序技术能够全面、快速地获得某一物种特定组织或器官在某一状态下几乎所有的转录本序列信息，以从整体水平研究其基因功能及结构，进而揭示其特定的生物学过程。该技术已广泛应用于植物候选基因发掘、功能鉴定及遗传改良等领域，为甘蔗响应低温胁迫的转录组研究提供了新手段。

本研究在前期田间农艺性状观测和生理生化研究的基础上^{[11, 12]}，筛选出3个抗寒能力差异比较明显的甘蔗品种作为研究材料，采用转录组测序技术研究不同抗寒能力的甘蔗品种在低温胁迫下的基因表达差异，挖掘与甘蔗抗寒性相关的基因，为进一步研究甘蔗抗寒机理和选育抗寒性强的优良甘蔗新品种提供理论依据和基因资源。

1 材料与方法 1.1 材料

抗寒性强甘蔗品种：广西壮族自治区农业科学院甘蔗研究所选育的桂糖28号(GT28)、广西壮族自治区亚热带作物研究所选育的桂热2号(GR2)。

低温敏感品种：新台糖22号(ROC22)。

1.2 方法 1.2.1 试验设计与管理

选择无病虫害、蔗茎大小均匀的甘蔗种茎切成单芽段，先用清水冲洗种茎并用干布擦洗干净，接着用50%多菌灵可湿性粉剂1 000倍液浸泡12 h消毒，再用蒸馏水浸洗种茎1 min后，用蒸馏水浸洗拧干的棉布将种茎包好并做好标记，然后用橡皮筋扎好放进温度为25 ℃的恒温箱催芽。当种茎萌芽并长出幼根时，将其移植到装有营养土的塑料盆并做好标记。每盆1段种茎，盆高17.5 cm，盆宽16.0 cm。育苗期间，每株施用完全营养液2次，每次10 mL。当幼苗两叶一心时，选择长势、大小均一的幼苗进行低温处理。低温胁迫处理组(T)温度为4 ℃，光照强度为5 000 lx，低温处理3 d；对照组(CK)温度为28 ℃，光照强度为5 000 lx。每个品种的处理组和对照组均设置3个重复。

1.2.2 叶片总RNA提取、文库构建及转录组测序

用植物总RNA提取试剂盒[天根生化科技(北京)有限公司]分别提取18个样本的总RNA，并用NanoDrop-2000分光光度计[赛默飞世尔科技(中国)有限公司]测定其A₂₆₀、A₂₈₀值, 然后经15%的变性聚丙烯酰胺凝胶电泳检测质量合格后送深圳华大基因科技服务有限公司进行cDNA文库的构建。以野生种割手密(Saccharum spontaneum L.)的参考基因组(http://www.life.illinois.edu, Sspon.v20190103)进行比对, 最多允许5个错配。统计reads在参考序列上的分布情况及覆盖度。经检验合格后，运用Mortazavi等^[13]提出的FPKM(Fragments Per Kilobase per Million)法计算基因表达量，用于比较不同样品间的基因差异表达。参照Audic等^[14]描述的方法, 对差异基因进行显著性分析。筛选的阈值为FDR (False Discovery Rate) ≤0.001，且倍数差异(Ratio)在2倍以上的基因被定义为显著差异基因。

1.2.3 基因共表达网络分析

使用权重基因共表达网络分析(Weighted Gene Co-expression Network Analysis，WGCNA)软件包进行基因共表达网络分析。用R语言DESeq2和edgeR包构建差异表达基因矩阵，其中在至少两个样本中CPM (Counts Per Million)值小于0的基因将被过滤掉。通过计算WGCNA包中相关函数的权重值(Pick Soft Threshold)，以确定最佳软阈值参数，使网络无限趋向于无尺度网络的分布。根据基因间表达量的相关性构建聚类树，采用动态切割法将产生的聚类树切割，把表达模式相似的基因合并在同一分支上，每个分支代表 1个共表达模块，将模块与性状进行关联分析^[15-17]。

1.2.4 生物信息学分析及抗寒基因筛选

根据WGCNA软件包分析与抗寒相关的模块，对模块中在3个甘蔗品种中均存在差异表达的已知基因进行GO富集分析，同时计算特异性模块的特征基因与差异表达基因之间的相关性(Module Membership, MM)、基因与表型性状的显著性(Gene Significance, GS)，筛选出与模块、性状高度相关(|GS|>0.8和|MM| > 0.8)的基因用于GO富集分析，筛选关键抗寒基因。

抗寒基因的筛选主要针对特异表达基因，即在抗寒品种GT28、GR2中差异表达，但在不抗寒品种ROC22中无差异表达的基因。参考彭振等^[18]的研究方法，根据低温胁迫下各品种基因表达量差异，设置较高的筛选标准(|log2 Ratio|≥4)，且是在参考基因组中获得注释的已知基因。

2 结果与分析 2.1 测序数据及质量情况

本试验共有18个样品，共获得785.99 M raw reads，过滤后共获得765.88 M clean reads，平均每个样本获得42.55 M clean reads。Q20平均值≥97%，Q30平均值≥93%，其中样品比对基因组的平均比对率为74.93%，比对基因的平均比对率为71.06%(表 1)；共检测到表达的基因为115 537个，其中已知的基因为73 844个，预测的新基因为41 693个；共检测出163 137个新转录本，其中76 263个属于已知蛋白编码基因的新的可变剪接亚型，42 863个属于新的蛋白编码基因的转录本，剩下的44 011个属于长链非编码RNA。

2.2 差异基因表达分析

低温胁迫下，GR2有20 870个上调表达基因，18 875个下调表达基因；GT28有18 467个上调表达基因，27 796个下调表达基因；ROC22有20 670个上调表达基因，20 699个下调表达基因。进一步对3个品种的特有差异表达基因及共有差异表达基因进行分析，结果表明，有16 844个差异表达基因为3个品种共有；GT28、GR2和ROC22中特有的差异表达基因分别为17 207个、6 263个和7 811个(图 1)。

2.3 WGCNA分析结果 2.3.1 软阈值的选择

由图 2可知，当软阈值为18时，R²值大于0.8。故选18作为WGCNA后续分析的软阈值。

2.3.2 模块鉴定

本试验共鉴定出13个模块，由不同的颜色代表。其中black模块包含588个基因，blue模块15 818个，brown模块12 642个，green模块3 038个，green-yellow模块72个，grey模块4 437个，magenta模块429个，pink模块577个，purple模块405个，red模块2 836个，salmon模块41个，tan模块66个，turquoise模块27 578个，yellow模块7 636个(图 3)。

2.3.3 模块与性状的关联分析

将获得的模块与性状进行关联分析，结果如图 4所示。本试验中所选的差异表达基因与品种(breed：3个品种)的相关性不强，而与处理方式(condition)相关。其中位于brown、purple和blue模块中的差异基因与处理方式呈正相关；而位于magenta和yellow模块中的差异表达基因则与处理方式呈负相关。

2.3.4 显著性模块中差异表达基因选择

由图 4可知，blue和yellow模块中的差异表达基因与低温处理存在中等强度的相关性，因此将这两个模块作为抗寒相关特异性模块重点关注。本试验在blue和yellow模块中鉴定出250个和466个差异表达基因。其中关键基因的选择标准：GS值和MM值均大于0.8，且P值小于0.05 (图 5)。

2.4 生物信息学分析及抗寒基因筛选 2.4.1 特异性模块GO富集分析及关键基因筛选

在blue和yellow模块中，3个甘蔗品种中均差异表达的已知基因分别有118个和220个。对差异表达的已知基因进行GO富集分析，发现blue模块全部差异基因在GO数据库中均得到注释，并显著富集到转座子核蛋白(Retrotransposon nucleocapsid)、细胞核(Nucleus)、DNA整合(DNA integration)等条目(图 6)。yellow模块的差异基因主要富集到叶绿体(Chloroplast)、质体(Plastid)、类囊体(Thylakoid)、光系统Ⅱ稳定(Photosystem Ⅱ stabilization)等条目(图 7)，这些都与光合系统有关。进一步根据MM值和GS值筛选相关性较强的关键抗寒基因，结果如表 2所示，共筛选到6个基因。

2.4.2 特异表达基因GO富集分析及关键基因筛选

在两个耐寒品种中，共筛选到204个特异表达的已知基因，其中有174个基因上调和下调表达的趋势一致，一致率达到85.3%。为了探索特异表达基因主要参与的生物途径和分子功能，对上述174个基因进行GO富集分析，发现其显著富集到乙醛酸循环(Glyoxylate cycle)、乙醛酸代谢过程(Glyoxylate metabolic process)、苹果酸合成酶活性(Malate synthase activity)、蛋白磷酸化(Protein phosphorylation)等生物过程(图 8)。筛选出表达差异倍数较大、可能与抗寒功能有关的基因如表 3所示。

3 讨论

低温是影响植物生长、生存及分布的重要环境因子，是造成作物减产的主要原因^[19]，特别是热带和亚热带的许多重要作物和水果，如甘蔗^[20]、香蕉(Musa nana Lour.)^[21]、菠萝蜜(Artocarpus heterophyllus Lam.)等^[22]，低温(0-12 ℃)条件会使植株的生长受到抑制、损害，甚至会导致其死亡。在作物响应低温胁迫过程中，其体内会发生快速而广泛的分子水平调节，主要涉及低温信号转导、转录、转录后修饰和翻译等多层次调节^[23]。本研究在转录水平分析抗寒性不同的甘蔗品种在低温胁迫下差异表达的基因，并通过WGCNA鉴定出13个模块，结合相关性分析筛选到blue和yellow两个模块作为甘蔗抗寒机理研究的目标模块。分析在两个抗寒品种中特异表达的基因，发现有174个基因上调和下调表达的趋势一致，一致率较高，达到85.3%，这些基因可能与甘蔗抗寒性相关。

对各模块内的基因进行GO富集分析，发现这些基因富集到多个不同的GO条目，表明植物对胁迫的应答及防御机制十分复杂，往往涉及多种协同信号转导和代谢途径，如细胞壁(Cell wall)、类囊体、光系统Ⅱ稳定、乙醛酸代谢过程等，这与前人的研究结果^{[24, 25]}一致，后续研究应重点关注。

Blue模块上调表达的基因Sspon.07G0019810-4D为转录辅阻遏物SEUSS (Transcriptional corepressor SEUSS)，SEUSS的功能主要是参与调控花的形成及发育^[26-28]，有关其在低温胁迫中的研究较少，如SEUSS参与生长素信号途径^[29]。低温胁迫下植物体内的激素水平通常会发生改变，进而调节生理代谢活动，以增强植物抗寒性^{[30, 31]}。刘世红等^[32]研究油棕(Elaeis gunieensis Jacq.)抗寒性与内源激素的关系，发现抗寒性高的品种其叶片生长素含量显著高于抗寒性低的品种，且随着温度的降低，生长素含量快速升高，表明SEUSS与植物抗寒性有一定的相关性。Sspon.05G0016040-4D为乙醇脱氢酶(Alcohol dehydrogenase)，当温度降到0-12 ℃时，乙醇脱氢酶会使细胞膜脂过氧化，甚至致使某些对低温敏感的植物细胞死亡，即发生低温胁迫，而乙醇可抑制植物细胞膜脂的分解，增加膜脂流动性，提高低温耐受能力^[33]。乙醇脱氢酶是乙醇生成途径中的关键酶，关于乙醇脱氢酶与低温胁迫的关系前人也做了许多研究^[34-36]。yellow模块下调表达的基因Sspon.03G0036590-2C为类光传感因子蛋白(Phosducin-like protein 3)，类光传感因子蛋白在植物方面的研究较少。Castellano等^[37]研究发现类光传感因子蛋白参与拟南芥(Arabidopsis thaliana)中微管的形成。类光传感因子蛋白的下调表达可能是因为低温使甘蔗细胞受到一定损伤，且在不抗寒品种中下调更为明显，说明这一性状受低温胁迫影响较大。基因的表达调控具有明显的时空特异性，上调表达或者下调表达均反映出基因参与调控过程。

特异表达基因中Sspon.05G0013050-2C、Sspon.08G0030000-1D、Sspon.01G0014100-3C和Sspon.07G0021990-3D在两个抗寒品种中差异表达倍数较大。Sspon.05G0013050-2C为硬脂酰-ACP脱氢酶(Stearoyl-acyl carrier protein Δ⁹ desaturase，SAD)，是控制植物细胞饱和脂肪酸与不饱和脂肪酸比例的关键酶，该酶使细胞膜的流动性增强，进而提高植物的抗逆能力^[38]。Sspon.08G0030000-1D为ATP依赖性RNA解旋酶(Probable pre-mRNA-splicing factor ATP-dependent RNA helicase DEAH4 isoform X2)，其作为剪接体的组成部分参与前mRNA剪接。RNA解旋酶在植物的生长发育过程和逆境响应途径中起着重要的调控作用^[39]。Sspon.01G0014100-3C为锌指蛋白(Zinc finger protein)，是植物中一类重要的转录因子，在基因的表达调控、细胞分化等生命过程中发挥重要作用^[40]。Sspon.03G0021430-3C和Sspon.04G0006870-3D的生物信息学分析显示其具有寒冷应答(GO: 0009409)和温度刺激应答(GO: 0009266)功能，可能参与甘蔗的抗寒调控途径。在植物中研究较为透彻的低温响应信号通路为CBF/DREB1信号途径^{[41, 42]}。本研究筛选出的特异表达基因中，Sspon.02G0014860-2B参与到CBF代谢通路(K22484)中，可重点关注研究。

植物对胁迫的应答过程涉及多个代谢途径和调控通路，是一个复杂的过程。本研究通过转录组测序技术分析不同抗寒性甘蔗在低温胁迫下差异表达的基因，为后续选育抗寒性强的优良甘蔗新品种提供理论依据和基因资源。