高温胁迫下2个棉花品种转录组可变剪切差异分析
丁　宇，王马寅，唐敏强，等．高温胁迫下２个棉花品种转录组可变剪切差异分析［Ｊ］．江苏农业科学，２０２３，５１（５）：１－１１．
ｄｏｉ：１０．１５８８９／ｊ．ｉｓｓｎ．１００２－１３０２．２０２３．０５．００１
高温胁迫下２个棉花品种转录组可变剪切差异分析
丁　宇，王马寅，唐敏强，李子阳，谢尚潜
（海南大学，海南海口５７０２２８）
  摘要：棉花（ＧｏｓｓｙｐｉｕｍｈｉｒｓｕｔｕｍＬ．）是一种重要的经济作物，是世界上第二大的天然纺织纤维来源和重要的食用油来源。棉花对生物和非生物胁迫高度敏感，尤其是高温胁迫极易影响花粉的活力和花药的开裂。为了解析棉花在高温胁迫下基因转录水平的相应变化机制，开展了热敏和耐热２个棉花品种的全长转录组响应高温胁迫处理变化的研究。通过转录本的可变剪切分析发现，２个棉花品种在高温胁迫下可变剪切的总数显著增加。热敏品种Ｃｈｅ６１－７２中发现了２９００个差异表达基因，并且差异基因在加热前样本（Ｒ０）和加热１２ｈ后的样本（Ｒ１２）中分别识别到了１３２５１个和２５２９６个可变剪切事件，其中内含子保留事件增加得最多，有３８３７个。耐热品种新陆早３６号中发现了２４３７个差异表达基因，在加热前样本（Ｔ０）和加热１２ｈ后的样本（Ｔ１２）中鉴定到了１１
２４８个和１３７６９个可变剪切事件，外显子跳跃事件变化得最大，增加了４１４４个。富集分析发现，２个品种的差异基因都显著富集到了光系统Ⅰ的光捕获、叶绿体类囊体膜和光合作用－天线蛋白通路中，并筛选出５个关键基因（ＣＰＢ３、Ａ０Ａ１Ｕ８ＩＺＦ２、Ａ０Ａ１Ｕ８ＫＣＡ２、Ａ０Ａ１Ｕ８ＮＤＷ４和Ａ０Ａ１Ｕ８ＮＩ７０），均被注释为叶绿素ａ／ｂ结合蛋白，它们参与了调控棉花光合作用动态平衡。本研究为棉花在高温胁迫的调节机制的深入研究提供了理论依据，为后续耐高温的种质改良及新品种培育提供了数据支持。  关键词：棉花；高温胁迫；全长转录组；纳米孔测序；可变剪切
  中图分类号：Ｓ５６２．０１  文献标志码：Ａ  文章编号：１００２－１３０２（２０２３）０５－０００１－１０
收稿日期：２０２２－０５－１７
基金项目：国家自然科学基金（编号：３２０６０１４９、３１７６０３１６）；海南省自然科学基金（编号：３２０ＲＣ５００、３２１ＲＣ４６９）。
作者简介：丁　宇（１９９６—），女，海南海口人，硕士研究生，主要从事转录组、基因组等生物信息学研究。Ｅ－ｍａｉｌ：ｗｏｓｈｉｄｉｎｇｙｕｙａ＠１６３．ｃｏｍ。
通信作者：谢尚潜，博士，教授，主要从事生物信息学和基于三代测序的基因组学和转录组学研究。Ｅ－
ｍａｉｌ：ｓｑｉａｎｘｉｅ＠ｆｏｘｍａｉｌ．ｃｏｍ。  棉花（ＧｏｓｓｙｐｉｕｍｈｉｒｓｕｔｕｍＬ．）是一种重要的经济作物，是世界上第二大的天然纺织纤维来源和重
要的食用油来源［１］。近年来，由于温室效应，地表温度升高对植物生长产生不同程度的影响［２］。棉
花在高温胁迫下花粉活力会受严重影响并导致花
药的开裂［３－４］
。因此，从不同水平了解棉花抗高温
胁迫分子机理，选育耐高温棉花材料是十分必要的。
随着测序技术的发展以及基因组研究的兴起，棉花抗高温胁迫已经开展了相关的研究。已有研究表明，棉花的酶活性、萌发、生长、根系发育、营养
发育和开花最适温度在２５～３１℃［５］
，高温环境通
过减弱棉花光合作用和新陈代谢，引起授粉能力减
弱［６］
，影响棉花花朵的定型、生产率和所产生的纤维的质量［７］。Ｌｉａｎｇ等通过纳米孔测序转录组分析
新疆棉花事件是怎么回事了棉花的抗热胁迫机制，揭示了富集在核糖体、激素信号转导以及内质网蛋白质加工途径的差异表达基因、叶绿体ａ／ｂ结合蛋白基因和ｌｎｃＲＮＡｓ对棉
花耐热机制的调节作用［８］
。Ｒｅｄｄｙ等研究发现，在
２０～３０℃的最适温度下，棉花生物产量增加，超过
最适温度则导致产量严重下降［９］
。高于最适温度会抑制棉花光合作用［７－１０］，同时伴随着碳水化合物
的进一步消耗和呼吸作用的增强［
１１］
。Ｂａｒｕａ等发现，热休克蛋白与棉花等作物的光合作用能力密切
相关［１２］
。Ｍｉｎ等发现，在高温敏感品系种，棉花花
药在高温诱导下ＧｈＣＫＩ高度表达，该基因可能在高
温胁迫中有重要作用［１３］。Ｗａｎｇ等研究发现，
ｍｉＲＮＡｓ通过调节与生长稳态、活性氧、叶绿体功能、植物－病原菌相互作用和植物激素信号转导途
径相关的基因来发挥作用［
１４］
。Ｍａ等揭示了ｌｎｃＲＮＡ对高温下棉花雄性生殖的作用，并发现了
ＧｈＨＲＫ１基因在高温胁迫反应中的负调控作用［１５］
。
这些研究结果为理解棉花抗高温胁迫的基因水平的调控机制提供了支持，但高温胁迫条件下关键调控基因转录本水平的调控解析仍有待进一步研究。
目前，三代测序技术已经被广泛运用于植物的全长转录本研究，主要包括鉴定复杂的可变剪切事
件、全长转录变异体、融合转录本和选择性多聚腺苷事件。可变剪切事件的发生是调控基因表达的一种重要调控机制，也是促进转录组和蛋白质组多样性的动力［１６］。三代测序相较于二代测序在序列长度、转录本识别方面的优势，更有助于识别复杂的剪接异构体，研究植物响应各种胁迫的基因调控机制［１７］。前人利用三代测序技术已经在高粱（Ｓｏｒｇｈｕｍｂｉｃｏｌｏｒ）［１８］、玉米（Ｚｅａｍａｙｓ）［１９］、异源多倍体棉花（ａｌｌｏｐｏｌｙｐｌｏｉｄｃｏｔｔｏｎ）［２０］和花竹（Ｐｈｙｌｌｏｓｔａｃｈｙｓｅｄｕｌｉｓ）［２１］中分别鉴定到可观的可变剪切事件。多倍体棉花中与耐热机制相关的可变剪切事件仍有待进一步研究。
本研究基于纳米孔测序的全长转录组数据，分析了热敏（Ｃｈｅ６１－７２）和耐热（新陆早３６号）２个品种的棉花在热刺激处理前后的可变剪切差异，发现高温胁迫增加了可变剪切事件，并进一步分析
了在高温胁迫中差异表达基因的功能，以及可变剪切事件的解析。本研究为棉花高温胁迫的调节机制的深入研究提供了理论依据，为后续耐高温的种质改良及新品种培育提供了数据支持。
１　材料与方法
１．１　数据收集与处理
收集来自ＮＣＢＩ的２个棉花品种（热敏品种Ｃｈｅ６１－７２和耐热品种新陆早３６号）的纳米孔测序数据，此数据由石河子大学在２０２１年３月４日上传发布，登录号为ＰＲＪＮＡ７０６６０３。测序样品使用光照培养箱（温度为４０℃）进行热应激处理，在处理０ｈ和１２ｈ后测序。使用ｐｉｐｅｌｉｎｅ－ｎａｎｏｐｏｒｅ－ｒｅｆ－ｉｓｏｆｏｒｍｓ程序（ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／ｎａｎｏｐｏｒｅｔｅｃｈ／ｐｉｐｅｌｉｎｅ－ｎａｎｏｐｏｒｅ－ｒｅｆ－ｉｓｏｆｏｒｍｓ）对纳米孔测序全长转录组数据进行处理，再用ｍｉｎｉｍａｐ２ｖ２．２４软件将测序数据比对到棉花的参考基因组（ｈｔｔｐｓ：／／ｗｗｗ．ｎｃｂｉ．ｎｌｍ．ｎｉｈ．ｇｏｖ／ｇｅｎｏｍｅ／？ｔｅｒｍ＝Ｇｏｓｓｙｐｉｕｍ＋ｈｉｒｓｕｔｕｍ）［２２］，参数设置为－ａｘｓｐｌｉｃｅ，－ｓｅｃｏｎｄａｒｙ＝ｎｏ－Ｃ５。
１．２　注释与可变事件的鉴定
将ｃＤＮＡ＿Ｃｕｐｃａｋｅ（ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／Ｍａｇｄｏｌｌ／ｃＤＮＡ＿Ｃｕｐｃａｋｅ）的ｃｏｌｌａｐｓｅ＿ｉｓｏｆｏｒｍｓ＿ｂｙ＿ｓａｍ．ｐｙ脚本用于比对数据去除冗余的转录本和合并同种转录本，使用ＳＱＡＮＴＩ３ｖ４．３软件［２３］对处理后的转录本数据结果与基因组注释进行比较分析。可变剪切事件通过ＳＵＰＰＡｖ２．３软件［２４］使用默认参数进行识别鉴定，ＳＵＰＰＡ可以从注释生成转录事件和局部选择性剪切事件，把可变剪切事件分为以下７种类型：外显子跳跃（ＳＥ），外显子互斥（ＭＸ），５′端选择性剪切（Ａ５），３′端选择性剪切（Ａ３），内含子保留（ＲＩ），替代第１个外显子（ＡＦ），替代最后１个外显子（ＡＬ）。
１．３　基因表达水平的量化和差异表达分析
通过Ｓａｌｍｏｎ定量软件［２５］将全长ｒｅａｄｓ比对到参考基因组上，并对匹配质量大于５的ｒｅａｄｓ进行量化。通过每百万计数（ＣＰＭ）估计基因表达水平（ＣＰＭ＝映射到转录本的ｒｅａｄｓ／样本中对齐的总ｒｅａｄｓ×１０６）。使用ＤＥＳｅｑ２Ｒ包［２６］对２个棉花品种不同时间点的样品进行差异表达分析。在该分析中，由ＤＥＳｅｑ２确定的Ｐ＜０．０１、｜ｆｏｌｄｃｈａｎｇｅ｜≥２被认为是差异表达基因（ＤＥＧ）。
１．４　差异表达基因的功能富集分析
使用ｅｇｇＮＯＧｖ２．０．０软件［２７］进行差异基因注释，从差异基因的ＧＯ数据库中获得注释
背景信息。利用Ｒ的ｃｌｕｓｔｅｒＰｒｏｆｉｌｅｒＲ包［２８］对差异基因进行ＧＯ富集分析，Ｐ值小于０．０５的ＧＯ条目被认为是差异表达基因显著富集。进一步使用ＫＯＢＡＳｖ３．０［２９］分析ＫＥＧＧ通路中ＤＥＧｓ的富集情况（α＝０．０５）。
２　结果与分析
２．１　数据统计分析
本研究从ＮＣＢＩ公共数据库中选取了棉花的耐热品种（新陆早３６号）和热敏品种（Ｃｈｅ６１－７２）的纳米孔测序的转录组样本进行分析，其中热敏品种Ｃｈｅ６１－７２对照组（加热前样本，Ｒ０）和试验组（加热１２ｈ后样本，Ｒ１２）总的ｒｅａｄｓ分别为１９３０１６７７和２１４５６６１７，总碱基数为１５．６Ｇ和１９．４Ｇ，ＧＣ含量为４３．２％和４２．５％。序列的平均长度为８４９ｂｐ和９５１ｂｐ，Ｎ５０长度为９１８ｂｐ和１０６４ｂｐ。在耐热品种新陆早３６号的对照组（加热前样本，Ｔ０）和试验组（加热１２ｈ后样本，Ｔ１２）中，分别获得了１９２１９３６０和１６５６４９４５条ｒｅａｄｓ，总的碱基数为１４．７Ｇ和１３ ６Ｇ，ＧＣ含量分别为４３．３％和４２ ６％，序列Ｎ５０长度为８６８ｂｐ和９３４ｂｐ，平均长度为８０１ｂｐ和９６４ｂｐ（表１）。
２．２　数据的处理与注释
比对２个不同品种棉花获得转录本，使用ｃＤＮＡ＿Ｃｕｐｃａｋｅ对全长转录组序列进一步去除冗余转录本。４个样品的转录组序列比对到基因组上的
表１　纳米孔测序数据统计概要
样本Ｒｅａｄｓ数碱基数
（Ｇ）ＧＣ含量
（％）平均长度
（ｂｐ）Ｎ５０长度
（ｂｐ）Ｒ０１９３０１６７７１５．６４３．２８４９９１８Ｒ１２２１４５６６１７１９．４４２．５９５１１０６４Ｔ０１９２１９３６０１４．７４３．３８０１８６８Ｔ１２
１６５６４９４５
１３．６
４２．６
８６４
９３４
比对率分别为９９．９３％、９９．９３％、９９．９１％和９９ ９３％（表２）。
  使用ＳＱＡＮＴＩ３ｖ．４．３对转录组进行分析和注
表２　全长转录本数据
样本比对率
（％）转录本数
（个）基因总数
（个）已知基因数
（个）新基因数
（个）Ｒ０９９．９１２０６１６０５１３４４４５４２８５９１６Ｒ１２９９．９３３７７５１８６５３４７５３０５６１２２９１Ｔ０９９．９３２１６８４４５３４０２４６５１１６８９１Ｔ１２
９９．９３
１９５７９７
５０４３７
４４７１１
５７２６
释，转录本总数最多的是热敏品种Ｃｈｅ６１－７２的试验组（Ｒ１２），得到了３７７５１８个转录本，基因总数也
最多，其中已知基因数为６５３４７个，新基因数量为１２２９１个。而转录本最少的是耐热品种新陆早３６号的试验组（Ｔ１２），有１９５７９７个。基因总数为５０４３７个，其中已知基因数为４４７１１，新基因数有５７２６个。对新的基因进一步进行注释，并用于下游分析。
在热敏品种Ｃｈｅ６１－７２的Ｒ０中每个基因只有一个转录本的现象占的比重最大，有３２．２５％，其次是具有２～３个转录本，占了２８．４９％。而在Ｒ１２中具有６个以上转录本的基因的比重最高，有３６ ３％。在耐热品种新陆早３６号的Ｔ０和Ｔ１２中，
每个基因具有１个转录本的现象比重最大（图１）。经过热应激处理后，２个品种中６个以上转录本的数量增加，可能是生物体内基因需要产生更多的可变剪切体来应对热应激反应。这些结果表明，使用纳米孔测序平台进行大规模转录本测序是重建转录本的有效方法，
有利于复杂转录事件的研究。
２．３　热应激处理后可变剪切事件的广泛动态变化
在高等动植物基因组中，大部分基因在转录过程中存在可变剪切，从而引起同个基因的多个转录本。植物在胁迫条件下会产生更多的剪切变异来
微调其基因的表达模式［３０］
。本研究使用了ＳＵＰＰＡ
软件去识别可变剪切事件，可变剪切事件分为７种类型：外显子跳跃（ＳＥ），外显子互斥（ＭＸ），５′端选择性剪切（Ａ５），３′端选择性剪切（Ａ３），内含子保留（ＲＩ），替代第１个外显子（ＡＦ），替代最后１个外显子（ＡＬ）（图２－Ａ）。在热敏品种Ｃｈｅ６１－７２的Ｒ０样品中，从３５１６２个基因中检测到了６６６３８个可变剪切事件，在Ｒ１２样品中，从６１０９７个基因中检测
到了１３７２４８可变剪切事件（图２－Ｂ）。研究结果表明，经过４０℃热应激的外部条件刺激１２ｈ后，Ｒ１２与Ｒ０相比表达的基因变多，可变剪切事件增加，特别是内含子保留（ＲＩ）增加了２７５９１，占了总的可变剪切数量的３０．０６％，其次是３′端选择性剪切（Ａ３），占了２８．４６％（图２－Ｃ）。在耐热品种新陆早３６号的Ｔ０样品中，从３１６８７个基因识别到了５９６３９个可变剪切事件，其中内含子保留（ＲＩ）事件最多。在Ｔ
１２样品中，从３４９７７个基因中识别到了６９４８３个可变剪切事件，而３′端选择性剪切事件最多，占了总事件的３０．７３％（图２－Ｂ和图２－Ｃ）。在４０℃热应激刺激后，Ｔ１２与Ｔ０相比表达的基因
和可变剪切事件都有所增加，但耐热品种对热应激的反应没有热敏品种这般强烈，最多的３′端选择性剪切事件增加了４０２３个。３′端选择性剪切（Ａ３）是棉花中最普遍的形式［３１］。已有研究表明，内含子保留（ＲＩ）可能是植物应对胁迫的重要分子机制［３２］。棉花在热应激处理后，内含子保留（ＲＩ）事件发生了显著变化，这表明内含子保留（ＲＩ）可能在响应高温胁迫的转录后转录组调整中发挥作用。
２．４　响应热应激的差异表达基因与可变剪切的变化
为了探究热应激条件下基因表达的差异，分析了热敏品种Ｃｈｅ６１－７２和耐热品种新陆早３６号在温度为４０℃条件刺激１２ｈ后的基因表达情况。在Ｃｈｅ６１－７２中发现了２９００个差异基因（ＤＥＧｓ），其中１５９２个基因表达上调，１３０８个基因表达下调。差异基因在Ｒ０样本中共识别到了１３２５１个可变剪切事件，在Ｒ１２样本差异基因中，可变剪切事件大量增加，一共识别到了２５２９６个可变剪切事件，其中ＲＩ事件增加得最多，有３８３７个（图３－Ａ）。在新陆早３６号中发现了２４３７个差异表达基因，其中１３６９个上调基因，１０６８个下调基因。在Ｔ０样本中鉴定到了１１２４８个可变剪切事件，在Ｔ１２样本中鉴定到了１３７６９个可变剪切事件，在热应激后，可变剪切事件增加了２５２１个，其中外显子跳跃事件变化得最大，增加了４１４４个。
从２个不同品种的棉花样本的可变剪切变化中研究发现，热敏品种Ｃｈｅ６１－７２对热应激反应非常
强烈，不同类别的可变剪切事件都呈上升趋势，而耐热品种新陆早３６号的不同类别的可变剪切事件除了外显子跳跃显著增加和外显子互斥少量增加之外，其他５个类型都有减少的现象，这可能是与耐热品种中耐热机制对热应激的自身调控有关。
  为了防止因为材料本身的差异导致的结果误差，本研究继续分析了热敏品种（Ｒ０）和耐热品种（Ｔ０）之间的基因表达的差异情况。一共发现了３１８个差异基因，其中２０６个基因表达上调，１１２个基因表达下调。其与热敏品种的热应激的表达差异基因中有３７个共有基因（图３－Ｃ），与耐热品种的热应激的表达差异基因中有３８个共有基因（图３－Ｄ）。在后续的功能注释分析时，
剔除共有差异基因。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。