2个四倍体栽培棉种冷胁迫响应基因的鉴定与比较分析

祁博文，郭梦兰，卢合均，梅欢，赵汀，方磊

（浙江大学农业与生物技术学院，杭州 310058）

棉花是天然可纺织纤维的主要来源，是一种重要的经济作物。陆地棉（Gossypium hirsutum）和海岛棉（G.barbadense）都是四倍体栽培棉种，由同一祖先驯化而来。陆地棉具有产量高、对逆境耐受能力好的特点，但是纤维品质一般；海岛棉的棉铃小，相对产量低，却能提供更细长强韧的优质纤维。陆地棉是最为广泛的栽培品种，占我国棉花种植面积的95%，而长绒的海岛棉只占1%[1]。棉花主要生长在热带和亚热带地区，是一种对低温敏感的作物[2]。随着生产需求的不断增加，棉花的种植区域渐渐转移到了气温较寒冷的高纬度地区[3]，更易受到冷害的影响。

冷胁迫对植物形态的破坏主要表现为植株明显矮化[4]，叶片失水萎蔫、表面出现黄斑或坏死，严重时可导致死亡[5]。植物遭受冷胁迫后，细胞膜的稳定结构被破坏，流动性受阻，渗透压发生改变，从而引起生长发育的失衡[6]。冷胁迫可以诱发植物细胞内产生大量活性氧（reactive oxygen species,ROS），使得植物内DNA、脂类等物质过氧化，破坏代谢平衡，因此，多数植物内部存在一套活性氧清除系统[7]。冷胁迫还会改变植物内源激素水平，从而引发一系列基因的表达形成反应机制，如脱落酸（abscisic acid,ABA）积累可以增加膜透性，诱导冷胁迫响应基因的表达[8]；生长素内源水平是植物耐冷性的另一考察指标，通常认为其内源含量与耐冷性呈负相关，这也许与生长素的极性运输有关[9]。

目前，已有许多针对棉花冷胁迫响应相关基因的报道。ZHENG 等的研究表明，冷胁迫可以使棉花纤维长度显著缩短，且可以诱导苹果酸代谢、可溶糖代谢、纤维素合成等途径基因的差异表达[10]。王俊娟通过对陆地棉抗冷品种的转录组测序，成功克隆到一个编码脱水素蛋白的基因GhDHNA1，该基因在抗冷棉花品种和冷敏感品种中有显著表达差异，且表达水平与耐冷性呈正相关[11]。LU等从棉花中鉴定到一个新的G蛋白受体基因，将其转入拟南芥后发现转基因株系在冷胁迫下表现出更高的萌发率[12]。HU 等对陆地棉TM−1 逆境处理材料进行转录组测序，共鉴定出8 725 个与耐冷性相关的基因，其中ABA信号转导途径相关基因，如ABA受体、SnRK蛋白和PP2C蛋白磷酸化激酶等在冷胁迫后显著上调表达[13]。目前的研究更多地是对单个棉花品种冷胁迫响应的基因进行鉴定，而对造成棉种间耐冷性差异的分子机制的研究较少。已有的研究结果显示，陆地棉和海岛棉具有冷害耐受性差异[13]，而其背后的分子机制依然未知。本研究通过对陆地棉TM−1 和海岛棉Hai7124 在冷胁迫处理后不同阶段的叶片进行转录组测序，旨在揭示棉花受冷胁迫诱导的基因，并比较二者冷胁迫响应基因的分子功能和代谢途径差异，以解析TM−1和Hai7124冷胁迫响应差异的分子基础，为棉花耐冷和优质新品种的培育提供基因资源。

1 材料与方法

1.1 棉花材料

2 种四倍体栽培棉花材料：陆地棉遗传标准系TM−1和海岛棉Hai7124。

1.2 棉花的冷胁迫处理与表型观察

将黏土和蛭石按照体积比2∶1 混合，分装于小号盆钵（直径10 cm）中。用水湿润后，选取长势均匀、饱满的棉花种子，在每个盆钵中播种3～4 粒。种子生长条件为温度22 ℃，14 h 光照/10 h 黑暗，相对湿度70%。待幼苗生长均匀并萌发至两叶一心时，分为2 组：一组为冷胁迫处理组，另一组为对照（CK）组。冷胁迫处理条件为生长于恒温培养箱中，温度4 ℃，14 h光照/10 h黑暗，相对湿度70%。对照组条件为生长于恒温培养箱中，温度28 ℃，其他条件与冷胁迫处理一致。

在冷胁迫施加后的不同时间点（0、1、3、6、12、24 h），分别取冷胁迫和对照组的5 株植株，进行快速拍照和表型观察，每个时间点的植株在表型观察后及时丢弃，不放回培养箱。

1.3 棉花的取样与RNA 提取

在冷处理开始前，分别从处理组和对照组选取生长较为一致的健壮幼苗作为0 h的材料。冷处理开始后，在1、3、6、12 和24 h 时分别从处理组和对照组取3 株长势一致且苗期棉花的第2 片真叶（即为3个生物学重复），立即放置在－80 ℃液氮中，用于RNA 的提取，以进行转录组测序。测序采用Illumina HiSeq 4000测序平台（美国Illumina公司），生成150 bp 的双端测序读长。总RNA 的提取使用RK16−50T 型植物RNA 快速提取试剂盒（南京钟鼎生物技术有限公司），具体操作按使用说明书进行。

1.4 原始数据的处理和冷胁迫下受诱导的基因的鉴定

得到fastq格式的原始数据后，使用软件fastp[14]过滤以获得高质量读长。使用HISAT×2（v2.1）[15]将高质量读长比对到TM−1 参考基因组（v2.1）或Hai7124 参考基因组（v1.1）（http://cotton.zju.edu.cn/download）上，并基于基因转移格式（gene transfer format, GTF）的注释文件和DESeq2 软件包鉴定差异表达基因（differentially expressed genes, DEGs），接着使用featureCounts 2.0.0[16]对基因表达进行定量，计算每百万读长中来自某一基因每千碱基长度的片段数目（expected number of fragments per kilobase of transcript sequence per million base pairs sequenced,FPKM），以评估基因的表达水平。差异基因的筛选标准为取伪发现率（false discovery rate,FDR）矫正后的P＜0.05 且差异表达倍数≥2（上调）或≤0.05（下调）的基因，上调和下调的DEGs均认为是受诱导的响应基因，并以FPKM=1 进行过滤，即FPKM＞1 的基因视为表达，FPKM＜1 的基因视为不表达。基因本体（gene ontology,GO）富集分析和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）富 集 分 析[17]使 用TBtools 1.046软件[18]。转录组原始数据在本课题组之前的研究中已上传至NCBI 数据库（https://www.ncbi.nlm.nih.gov/）[13]，数据编号为PRJNA490626。

1.5 实时荧光定量反转录聚合酶链反应（real time quantitative reverse polymerase chain reaction,RT-qPCR）验证

以Histone 3作为内参基因，引物为Y8991（F，5´−CGGTGGTGTGAAGAAGCCTCAT−3´；R，5´−AA TTTCACGAACAAGCCTCTGGAA−3´）。每个样品设置3 个技术重复。使用Primer−Blast 工具（http://www.ncbi.nlm.nih.gov/tools/primer−blast）设计引物。引物序列如表1 所示。RT−qPCR 使用的是HiScript®ⅡQ RT SuperMix 试剂盒（南京诺唯赞生物科技股份有限公司），反应在不含RNase 的PCR反应管中进行：向反应管中加入RNA 模板（500 ng）、gDNA wiper 混合物（4 μL），并用无RNA 酶的双蒸馏水补至16 μL，于42 ℃条件下反应2 min 后，加入HiScript®ⅡqRT SuperMix Ⅱ（4 μL），然后经过25 ℃引物结合10 min，50 ℃延伸30 min，85 ℃酶失活5 min 的程序，以此来合成cDNA 的第一条链。使用引物Y347 对得到的cDNA 进行PCR 扩增和电泳检测，以验证cDNA样品的质量。

表1 用于定量检测、基因克隆和载体构建的引物序列Table 1 Primer sequences for quantitative detection,gene cloning and vector construction

实时荧光定量PCR 使用的荧光染料是SYBR Green Master Mix（南京诺唯赞生物科技股份有限公司），在StepOnePlus 型PCR 仪（美国Thermo FisherScientific 公司）上进行反应，方法为2 步法PCR 循环，基因的相对表达水平通过2－ΔΔCT方法计算得到。

2 结果与分析

2.1 TM-1和Hai7124在冷胁迫下的响应表型观察

在冷处理前，陆地棉TM−1叶片生长状态良好。冷胁迫处理1 h后，未产生明显的表型变化；3 h后，TM−1 的子叶表面较为平整，真叶叶尖出现轻微下垂；6 h 后叶片逐渐失去光泽，子叶在萎蔫的同时出现轻微卷缩，真叶下垂程度加剧，叶片出现向背面卷缩；12 h后整个植株萎蔫明显，整个子叶从叶柄处完全下垂，真叶叶面的卷缩程度进一步加剧，但真叶仍能勉强保持大部分叶面平整；24 h后，整个植株萎蔫，真叶叶面不再平整，皱缩严重，幼嫩主茎开始出现较为严重的倒伏现象（图1）。

海岛棉Hai7124在4 ℃冷胁迫下的表型变化趋势和TM−1 一致，但更快地出现了表型变化。在冷胁迫下处理3 h后，相比TM−1，Hai7124真叶萎蔫下垂的程度较为明显；6 h后二者的表型差异进一步加剧，且Hai7124的子叶卷缩更加明显，子叶下垂程度更加严重；12 h 后Hai7124 幼嫩的主茎出现轻微倒伏，真叶皱缩严重，但TM−1的幼嫩主茎仍能维持直立状态，子叶仍能保持较为平整的叶面；24 h后整个Hai7124 植株的子叶叶柄下垂程度、叶片皱缩程度均明显超过TM−1。通过对冷胁迫处理不同时间的表型比较发现，TM−1比Hai7124对冷胁迫的耐受力更强，冷胁迫表型更晚出现（图1）。

图1 TM-1和Hai7124棉花幼苗在冷胁迫下的表型比较Fig.1 Comparisons of phenotypes of TM-1 and Hai7124 cotton seedlings under 4 ℃treatment

采用Illumina HiSeq 高通量测序技术对冷胁迫处理和对照组的TM−1 和Hai7124 的幼苗进行转录组测序，结果显示：平均每个样本产生54.69百万个读长，转录组比对率在93.9%～97.8%之间。Q20值（碱基错误识别概率小于1%）均在98.5%以上，Q30值（碱基错误识别概率小于0.1%）均在96%以上。使用DESeq2软件包比对处理组和对照组的测序数据，得到TM−1和Hai7124在冷胁迫处理后各个时间点的DEGs（图2）。根据不同处理时间进行统计，得到TM−1 在冷胁迫处理后的1、3、6、12、24 h 分别有745、853、1 349、3 804、4 693 个DEGs，TM−1 响应冷胁迫的基因数随处理时间延长呈现上升趋势；Hai7124 在冷胁迫处理后的1、3、6、12、24 h 分别有113、3 936、3 085、5 674、4 301个DEGs，说明Hai7124在冷胁迫处理初始的反应比TM−1缓慢，而在3 h时有大量基因受诱导，说明3 h时Hai7124已经对冷胁迫做出强烈的反应，有较多基因差异表达。

图2 冷胁迫下各个时间段TM-1和Hai7124的差异表达基因数Fig.2 Numbers of differentially expressed genes (DEGs) of TM-1 and Hai7124 at different stages under cold stress

一个基因如果在1、3、6、12、24 h中的任意一个时间点受到诱导，该基因即被视为总体上受诱导的基因。以此统计了TM−1 和Hai7124 中受诱导的基因数目。从图3A可知：冷胁迫处理下TM−1中共有8 713 个基因受诱导，63 649 个基因不受诱导；Hai7124中共有12 461个基因受诱导，62 065个基因不受诱导。通过对TM−1 和Hai7124 之间的直系同源基因比较发现，两者之间有3 884 个共同受诱导的基因，TM−1 拥有1 518 个特异受诱导的基因，而Hai7124拥有3 577个特异受诱导的基因（一些基因并没有对应的直系同源基因）（图3B）。

图3 受诱导基因数Fig.3 Numbers of induced genes

对TM−1中1 518个特异受诱导的基因进行GO富集分析，结果显示这些基因多富集在信号受体活性、光合作用、蛋白质代谢过程、细胞蛋白修饰过程、对光刺激反应、前体代谢产物和能量产生、脂质代谢过程等功能中（图4A）；KEGG 富集分析显示，这些基因主要集中在光合作用−天线蛋白、昼夜节律——植物、光合作用、信号转导、植物激素信号转导、环境信息处理、半胱氨酸和蛋氨酸的代谢、磷酸肌醇代谢、三羧酸循环等代谢途径中（图4B）。

图4 TM-1中特异受诱导基因的GO和KEGG富集分析Fig.4 GO and KEGG enrichment analysis of specifically induced genes in TM-1

对Hai7124 中3 577 个特异受诱导的基因进行GO富集分析，结果显示，这些基因多富集在水解酶活性、催化活性、叶绿体、质体、类囊体、脂质代谢过程等功能中（图5A）；KEGG 富集分析显示，这些基因主要集中在各种代谢过程中，如脂质代谢，甘氨酸、丝氨酸和苏氨酸代谢，鞘脂代谢，β−丙氨酸代谢，组氨酸代谢，精氨酸和脯氨酸代谢，赖氨酸降解等（图5B）。

图5 Hai7124中特异受诱导基因的GO和KEGG富集分析Fig.5 GO and KEGG enrichment analysis of specifically induced genes in Hai7124

2.2 冷胁迫下脂肪酸代谢途径种间差异比较分析

脂肪酸是植物细胞的必需组分之一，是植物激素等许多重要生物分子的中间合成体，也是细胞能量储存的基础。脂肪酸的含量和组成能够影响细胞膜的流动性，从而影响植物对低温环境的耐受能力。根据现有研究报道，我们从TM−1 和Hai7124里鉴定出了脂肪酸合成途径关键家族中的乙酰辅酶A羧化酶（acetyl−CoA carboxylase,ACC）、脂肪酸合成酶（fatty acid synthase,FAS）、十八烷−ACP去饱和酶（stearoyl−ACP desaturase, SAD）、脂肪酸去饱和酶（fatty acid desaturase, FAD）、长链酰基辅酶A合成酶（long chain acyl−CoA synthetase, LACS）、酰基−ACP硫酯酶（acyl−ACP thioesterase,FAT）的基因并绘制了棉花脂肪酸代谢途径的网络图（图6～7）。在脂肪酸合成途径关键基因家族中，在TM−1 中共鉴定到91 个基因，其中有17 个受冷胁迫诱导；Hai7124中共鉴定到84个基因，其中有24个受到冷胁迫诱导。在TM−1中，我们发现16个ACC家族基因中有2个基因受冷胁迫诱导；26个FAD家族基因中有10 个受冷胁迫诱导；6 个FAT 家族基因中有2个受冷胁迫诱导；24 个LACS 家族基因中有3 个受冷胁迫诱导。共有11 个基因在整个冷胁迫过程中呈现上调表达趋势，也有一些基因展现出波动变化的表达模式（图6）。

图6 TM-1在冷胁迫下脂肪酸代谢通路图和基因表达热图Fig.6 Diagram of fatty acid metabolism pathway and heat maps of expressed genes in TM-1 under cold stress

相比之下，在Hai7124中，我们发现了4个ACC家族基因、11个FAD家族基因、2个FAT家族基因、4个LACS家族基因以及3个SAD家族基因受冷胁迫诱导。在冷胁迫下，只有一小部分基因在某些时间段上调表达（如FAD家族基因），而有10个基因在整个冷处理过程中呈下调表达趋势，也有一些基因呈波动表达模式（图7）。可以推测，在冷胁迫下，TM−1中有更多的脂肪酸代谢途径关键基因上调表达；而Hai7124 中则相反，包含更多下调表达的基因。这些受诱导的脂肪酸代谢途径基因的差异表达模式可能造成了TM−1和Hai7124对冷胁迫的耐受差异。

图7 Hai7124在冷胁迫下脂肪酸代谢通路图和基因表达热图Fig.7 Diagram of fatty acid metabolism pathway and heat maps of expressed genes in Hai7124 under cold stress

从TM−1 和Hai7124 脂肪酸代谢通路中挑选了5 个冷胁迫响应基因进行RT−qPCR 验证，结果（图8）表明，GH_A13G2126（GhFAD6−1）、GH_D06G0458（GhLACS9）、GB_D13G2209（GbFAD6）、GB_D03G0951（GbLACS4）表 达 整 体 上 调，GH_D13G2107（GhSFD6−2）表达先上调后下调。

图8 5个脂肪酸代谢途径基因差异表达的RT-qPCR验证Fig.8 RT-qPCR validation of five differentially expressed genes involved in fatty acid metabolism pathway

3 讨论与结论

冷胁迫对棉花生产的影响巨大，我国有5%～12%种植面积的棉花是因为受冷害而补种或重播种的[19]。因此，解析陆地棉和海岛棉冷胁迫响应差异形成的分子机制，对培育高产、高品质且抗逆性好的棉花新品种极为重要。已有的研究证明，陆地棉对冷害的耐受性更好。在本研究中，经过冷胁迫处理后3 h 时，陆地棉TM−1 和海岛棉Hai7124 的表型变化开始出现差异，此时Hai7124 的真叶萎蔫更明显，到6 和12 h 时表型差异更大，Hai7124 的叶片萎蔫和主茎倒伏更加剧烈，说明TM−1 比Hai7124对冷胁迫的耐受能力更强。

目前，通过同源克隆等方法，已经在棉花中鉴定出了许多冷胁迫相关响应基因，如GhERF1[20]、GhAGP[21]、GhFAD[22−23]等。然而棉花的耐冷、耐旱等抗逆性状大多由多基因控制，无法通过单基因的变化来提升抗逆性，因此越来越多的研究人员致力于从全基因组或全转录组水平上对棉花的冷胁迫响应机制进行研究。ZHANG等从全基因组水平对陆地棉TM−1 进行正向选择基因研究发现，驯化后的TM−1中，A亚基因组在纤维发育方面受到正向选择，而D 亚基因组则在逆境耐受方面受到正向选择[24]。LU等通过对亚洲棉、雷蒙德氏棉和陆地棉的全基因组MATE 基因家族进行分析，分别鉴定出了70、72和128 个MATE 家族基因，将基因Gh_D06G0281在拟南芥中过表达后，显著提高了拟南芥植株对冷胁迫的抗性[25]。这些结果表明陆地棉经过正向选择和人为驯化而拥有了更好的环境耐受力。本研究中，我们从转录组水平鉴定出了冷胁迫处理后TM−1和Hai7124 中总体受诱导的基因，并通过直系同源基因的比较鉴定出了两者之中特异受诱导的基因。TM−1 和Hai7124 在冷胁迫处理后受诱导的基因随处理时间增加而增加，Hai7124 在3、12 h 时间点的受诱导基因数量增加明显。冷胁迫下，TM−1 中有1 518个特异受诱导的基因参与激素信号转导、昼夜节律等调节过程，而Hai7124中有3 577个特异受诱导的基因参与赖氨酸、精氨酸、丝氨酸、苏氨酸等氨基酸代谢，这反映了TM−1 和Hai7124 在冷胁迫下受诱导基因的分子功能和代谢途径的差异。

除了冷胁迫响应的分子机制以外，棉花的耐冷研究还包括耐冷性的测定以及冷胁迫下的各类生理指标等[26]。细胞膜的渗透性和流动性在一定程度上反映了植物在冷胁迫下的受损害程度，也可以从侧面反映植物对冷胁迫的耐受能力。而细胞膜的渗透性、流动性与膜脂中的不饱和脂肪酸比例密切相关。提高细胞膜内不饱和脂肪酸含量的一个重要方法就是饱和脂肪酸的去饱和反应，主要通过去饱和酶催化膜内脂肪酸来实现。在本研究中，筛选出TM−1 和Hai7124 里的冷胁迫相关基因后，我们着重关注了脂肪酸相关通路中的基因受冷胁迫诱导的表达模式差异：TM−1 脂肪酸途径相关受诱导基因中共有11 个基因在整个冷处理过程中呈上调表达趋势；而Hai7124 中共有10 个基因呈下调表达趋势。TM−1 和Hai7124 这些脂肪酸途径相关基因在冷胁迫下的差异表达，可能是造成其耐冷性差异的生物学基础。

陆地棉和海岛棉经过了同一起源、独立驯化的过程，形成了两者在产量、品质和适应性方面的差异。本研究从陆地棉TM−1 和海岛棉Hai7124 这2个栽培品种中鉴定出冷胁迫响应基因，为未来棉花耐冷机制的研究提供了基因资源；同时解析了TM−1和Hai7124冷胁迫响应基因在分子功能、代谢途径以及表达模式上的差异，这可能与两者适应性差异的形成密切相关，然而还需要在未来进行进一步的功能验证。本研究为创制转基因新材料并培育耐冷、高产、优质的棉花新品种提供了理论基础。