棉花GR基因家族的全基因组鉴定及分析

张传义，许艳超，蔡小彦，王星星，侯宇清，王玉红，王坤波，刘方，周忠丽

（棉花生物学国家重点实验室/ 中国农业科学院棉花研究所，河南安阳455000）

植物在生理活动中经过体内氧化还原反应会产生过氧化氢、超氧阴离子、羟自由基和单线态氧等活性氧物质[1-2]。正常条件下，植物体内产生的活性氧保持动态平衡状态[3]。然而在干旱、高温、盐碱等非生物胁迫下，活性氧物质的动态平衡被打破，植物会产生大量的活性氧物质，对植物造成氧化胁迫[4]，进而破坏细胞膜完整性、导致膜脂过氧化[5]，改变蛋白活性[6]，甚至导致细胞死亡，严重影响作物的生长发育[7-9]。

谷胱甘肽是含有巯基的抗氧化物质，具有储存和转运还原态硫、调节酶活性[10]、传导氧化还原信号的功能[11]。研究表明谷胱甘肽能够有效地防止植物中活性氧 （Reactive oxygen species，ROS）积累，植物对非生物胁迫忍耐性与体内谷胱甘肽的含量及其氧化还原状态紧密相关[12]。谷胱甘肽还原酶（Glutathione reductase，GR）是 1 种黄素蛋白氧化还原酶，主要借助烟酰胺腺嘌呤二核苷酸磷酸酯 （Nicotinamide adenine dinucleotide phosphate，NADPH）提供电子进而将氧化型谷胱甘肽（Glutathione oxidized，GSSG）还原为还原型谷胱甘肽（Glutathione，GSH），维持植物体 GSH/GSSG的比例，提高植物对非生物胁迫的耐受性[13-14]。GR 普遍存在于叶绿体中，少量存在于细胞质、线粒体、过氧化物酶体中[15-16]。对其结构研究发现GR 主要包括黄素腺嘌呤二核苷酸 （Flavin adenine dinucleotide，FAD） 结构域和 NADPH 结构域，其中以同型二聚体的形式存在为主，但在衣藻中却以单体存在[17]，在玉米中以异二聚体存在[18-19]。GR 具有多种由不同基因编码的同工酶，对提高生物胁迫和非生物胁迫抗性具有重要的意义[20-21]。在盐胁迫下，柑橘[22]、豌豆[23]、水稻[24]、大豆[25]、番茄[26]、小麦[27]、谷子[28]等物种的 GR 活性提高；同样在干旱胁迫下，玉米[29]、小麦[30]、水稻[31-32]等物种的GR 活性也会提高。在低温胁迫下，水稻幼苗胞质GR基因的转录水平提高[33]。在强光下，对抑制GR表达的转基因拟南芥研究发现，光系统Ⅱ活性被抑制，植物体内H2O2含量升高，植物生长发育受到影响[34]。在臭氧和百草枯引起的氧化胁迫下，过表达GR基因能够提高烟草光合作用[35-36]。在低温胁迫下，抑制番茄GR基因表达，植物体内H2O2含量升高，增强了对冷害胁迫的敏感性[37]。敲除水稻GR3 基因增加对盐胁迫的敏感性[38]。过表达GR基因的棉花对短期低温胁迫的耐受性提高[39]；在低温和强光照条件下，叶绿体中过表达GR基因的棉花与野生型比较，具有更高的光合速率[39-40]；在盐胁迫下，棉花的GR 活性提高，抗氧化性增强[41]。

目前雷蒙德氏棉 （Gossypium raimondii）、亚洲棉（G.arboreum）、陆地棉（G.hirsutum）、海岛棉（G.barbadense）基因组序列已先后被解析。本研究通过生物信息学方法，对GR基因家族成员进行鉴定，进一步分析其理化性质、序列特征、系统进化关系、染色体位置，并对陆地棉中GR基因的表达模式进行研究，进而为后续研究GR基因功能提供基础。

1 材料与方法

1.1 棉花GR基因家族的鉴定

在 Cottongen 数据库（https://www.cottongen.org/） 下载陆地棉 （AD1，HAU）、 海岛棉（AD2，HAU）、雷蒙德氏棉 （D5，JGI_v2.1）和 亚洲棉（A2，CR1）基因组数据。通过蛋白质家族比对和隐马尔可夫模型 （Hidden Markov Model，HMM）数据库 （Protein families database of alignments and HMM，Pfam 数据库）网站（http://pfam.xfam.org/）下载GR基因的种子文件PF07992 和PF02852。使用HMMER 3.0 软件和本地BLASTP 程序比对搜索含有GR 蛋白结构域的序列。将筛选出的蛋白序列使用在线软件Pfam、简单模块构架搜索工具 （Simple modular architecture research tool，SMART，http://smart.embl-heidelberg.de） 和保守结构域数据库（Conserved domain database，CDD，http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi）进行比对验证。

1.2 序列分析与系统进化树构建

利用ClustalW 软件对鉴定得到的陆地棉、海岛棉、雷蒙德氏棉、亚洲棉GR 蛋白序列进行比对，使用Mapman 绘制棉花多序列比对结果。利用 MEGA 6.0 软件对拟南芥（Arabidopsis）、可可（Theobroma cacao） 及3 个棉花基因组中提取的GR 家族成员进行多序列比对，使用邻接法（Neighbor-joining）构建系统发育树，Bootstrap 设置为1 000。利用在线工具ProtParam（http://web.expasy.org/protparam/）对GR 蛋白质的理化性质进行分析。使用在线工具ProtComp9.0（http://www.softberry.com/berry.phtml?topic=protcomppl&group=programs&subgroup=proloc）预测 GR 家族成员的亚细胞定位。使用InParanoid 软件预测陆地棉、海岛棉、雷蒙德氏棉、亚洲棉GR基因家族成员同源基因对。使用KaKsCalculator 计算同义突变率（Ks）和非同义突变率（Ka）。

1.3 染色体定位与基因结构分析

使用在线软件 MEME（http://meme.sdsc.edu/meme/）对保守基序进行分析，主要参数设置：number of unique motifs: 10；maximum and minimum search widths：50。在 4 个棉种（G.hirsutum、G.barbadense、G.raimondii、G.arboreum）的基因组 GFF（General feature format）注释文件中提取GR基因家族成员位置及结构信息，使用MapChart2.2 软件进行分析绘制染色体定位，利用基因结构显示系统 （Gene Structure Display Server，GSDS，http://gsds.cbi.pku.edu.cn/）绘制GR基因家族基因结构。

1.4 GR家族基因的表达模式分析

在 NCBI 的 SRA （Sequence read archive）数据库（http://www.ncbi.nlm.nih.gov/sra）下载陆地棉的转录组数据 （PRJNA248163）。使用 Tophat和Cufflink 软件包处理数据并计算表达量，使用R 语言热图绘制包（Pheatmap）绘制陆地棉GR基因表达谱，表达量为FPKM（Expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced，每百万片段中来自某基因每千碱基长度的片段数）。

2 结果与分析

2.1 棉花GR基因家族的鉴定

使用本地 BLASTP 程序和 HMMER 3.0 软件对种子文件进行比对搜索后得到候选序列，并且利用Pfam、SMART、CDD 软件验证保守结构域后，将不含完整GR 结构域的序列舍弃，最终鉴定得到18 个GR基因，陆地棉、海岛棉、亚洲棉、雷蒙德氏棉中分别有 6 个、6 个、3 个、3 个，分别 命 名 为GhGR1 ～GhGR6、GbGR1 ～GbGR6、GaGR1～GaGR3 和GrGR1～GrGR3（表 1）。GR蛋白质长度为 421～577 个氨基酸残基 （aa），相对分子质量在 45.807～63.072 kDa，等电点为5.9～8.459。亚细胞定位预测结果表明，GR 在4个棉种的叶绿体和细胞质均有分布，且主要定位于叶绿体（表1）。

2.2 序列分析与进化分析

GR基因家族系统发育分析发现GR基因主要分为2 个亚组，即Ⅰ组和Ⅱ组，分别包括12 和6 个成员。对GR基因的结构分析，发现各组内成员间具有相似的基因结构并且编码相似的结构域。对GR基因的外显子（Exon）数目统计发现，2个亚组的外显子数目具有明显差异：Ⅰ组基因均具有 16 ～18 个外显子；Ⅱ组除基因 Gbar_D11G013900 和 Gorai.007G142200 具 有 9 个 外显子外，其余基因均具有11 个外显子。GR基因序列比较发现，Ⅰ组基因结构序列较Ⅱ组明显变长（图1）。同一组内具有相似的外显子数目且具有相似的基因序列结构，表明系统进化分类的可靠性。

利用MEME 软件对GR基因家族编码产物的保守基序进行分析，发现除Ⅱ组基因Ghir_D11G013530、Ghir_A11G013490 和 Gbar_D11G013900 外，其他GR基因的编码产物均含有8 个保守的模体（motif），表明 GR 蛋白结构域高度保守（图1 和图2）。亚细胞定位分析发现在Ⅰ组中，6 个 GR 定位于叶绿体，6 个 GR 定位于细胞质；在Ⅱ组中，6 个 GR 均定位于叶绿体。GR 高度保守的结构域表明，不同亚组间GR基因具有相似的功能，但其编码产物作用的细胞位置不同。

对陆地棉、海岛棉、雷蒙德氏棉、亚洲棉、拟南芥（包含2 个GR基因家族成员）和可可（包含2 个GR基因家族成员）GR基因家族编码的蛋白序列构建系统进化树，结果发现，可将GR基因家族分为2 类（图3），与棉花GR基因系统发育分析结果相同，表明在拟南芥、可可、棉花物种分化前，GR基因家族包含2 个亚组。此外发现可可基因Thecc1EG020353 和拟南芥基因AT3G24170 被归类于亚组Ⅰ；Thecc1EG046980和AT3G54660 归类于亚组Ⅱ，并且棉花GR基因与近缘种可可的亲缘关系明显高于与拟南芥。表明GR基因在近缘种中具有更高的相似度，基因的功能可能更相似。

表1 棉花GR基因家族及其编码产物信息Table 1 The information of GR gene family and their products in cotton

图3 棉花、可可和拟南芥GR 系统进化分析（A）与同源关系分析（B）Fig.3 Phylogenetic analysis of GR in cotton, cocoa and Arabidopsis (A) and analysis of homology of GR gene family in cotton (B)

2.3 基因定位与同源性分析

陆地棉的6 个GR基因定位在6 条染色上体，海岛棉的 6 个GR基因定位在 6 条染色体上，雷蒙德氏棉3 个GR基因定位在3 条染色体上； 亚洲棉3 个GR基因定位在 3 条染色体上（图4）。

通过InParanoid 软件分析GR基因同源基因对，共发现15 对同源基因对（图3）。同源基因的Ka/Ks值，主要分布于 0.001 0～0.784 9 （表2）。Ka/Ks值均小于1，表明陆地棉、海岛棉、亚洲棉和雷蒙德氏棉中的GR基因家族成员在进化过程中经历了较强的纯化选择过程。

2.4 陆地棉GR基因的表达模式分析

对在低温胁迫、高温胁迫、PEG6000 模拟干旱胁迫、盐胁迫处理条件下的陆地棉GR基因进行表达谱分析，结果发现在干旱胁迫和盐胁迫处理后，GR基因的表达模式相似 （图 5）。其中Ghir_D02G012650、Ghir_A11G013490、Ghir_D11 G013530 在干旱胁迫和盐胁迫处理后先上调表达，在 3 h 后出现明显下调表达；而 Ghir_A03 G010510、Ghir_A10G012520、Ghir_D10G015270在干旱胁迫和盐胁迫处理后上调表达。在低温胁迫、 高温胁迫处理后，GR基因表达模式较为一致：Ghir_A03G010510、Ghir_A10G012520、Ghir_D10G015270 先下调表达，在处理12 h 后又上调表达；Ghir_D02G012650、Ghir_A11G013490、Ghir_D11G013530 在低温胁迫处理1 h 时上调表达，而后下调表达； 在高温处理后，Ghir_D02G012650、Ghir_A11G013490、Ghir_D11G013530 随 着 处 理时间的增加均出现下调表达。Ghir_A11G013490和Ghir_D11G013530 作为同源基因，一个来自于A 基因组，另一个来自于D 基因组，在不同的胁迫处理后均表现出最为相似的表达模式，表明该同源对在后续的棉花进化中功能未发生变化。

3 讨论

图4 GR 基因在染色体上的分布Fig.4 Chromosomal distribution of GR genes

表2 GR 基因家族同源基因对Ka/Ks 值Table 2 The Ka/Ks ratio of homologous gene pairs of GR genes

图5 陆地棉 GR 基因在低温胁迫（A）、高温胁迫（B）、PEG-6000 模拟干旱胁迫（C）、盐胁迫（D）处理条件下的表达模式Fig.5 The expression patterns of upland cotton GR genes under cold stress (A), high temperature stress (B),PEG-6000 simulated drought stress (C) and salt stress (D)

GR 作为植物体重要的抗氧化酶[42]，对提高生物胁迫和非生物胁迫抗性具有重要的作用[20-21]。大量的研究表明，在盐碱、干旱、低温等非生物胁迫下GR 活性增加，可明显提高植物的抗性[22-32]。GR基因的过表达提高了转基因植株对生物胁迫和非生物胁迫的抗性[33-37]。通过比对鉴定，本研究在陆地棉、海岛棉、亚洲棉、雷蒙德氏棉中分别鉴定出 6 个、6 个、3 个、3 个GR基因。陆地棉和海岛棉GR基因数目一致，亚洲棉和雷蒙德氏棉GR基因数目一致，且具有相似的基因结构，表明GR基因在陆地棉和海岛棉、 雷蒙德氏棉和亚洲棉之间高度保守并且具有相似的功能，同时表明亚洲棉和雷蒙德氏棉二倍体棉种的共同祖先可能存在相似的GR基因成员分布。对棉花系统发育分析发现，GR基因家族成员分为2 个亚组，各亚组的基因结构、外显子数目相似。异源四倍体陆地棉包含A、D 亚基因组，且棉花A 基因组的亚洲棉为A 亚组供体种，D 基因组的雷蒙德氏棉为 D 组的供体种[43]。本研究中陆地棉 6 个GR基因分别与亚洲棉和雷蒙德氏棉的3 个基因归为6对直系同源基因，表明在进化过程中异源四倍体陆地棉的GR基因没有经过单独的复制或明显的基因丢失事件，具有高度保守性。

系统进化分析表明，GR基因家族在可可、拟南芥和棉属分化的长期过程中保留了2 个亚组，也证明了这2 个亚组在植物生长发育与环境适应过程中的重要性。其结果还表明GR基因家族亚组Ⅰ和亚组Ⅱ的分化可能早于棉花、 可可、拟南芥的物种分化。Ka/Ks值可以用来判断是否有选择压力作用于基因。如果Ka/Ks＞1，则认为有正选择效应；如果Ka/Ks＝1，则认为存在中性选择；如果Ka/Ks＜1，则认为有纯化选择作用[44]。本研究同源基因对的Ka/Ks均小于1，表明GR基因在进化过程中高度保守，具有相似的功能。

在棉花中过表达GR基因可提高对短期低温胁迫的耐受性[39]，在盐胁迫处理下，谷胱甘肽还原酶活性提高，增强棉花的抗氧化性[41]。对陆地棉叶片中GR基因的表达模式进行分析，结果发现所有的GR基因在非生物胁迫下均积极响应非生物胁迫，表明在非生物胁迫中GR 起到关键的作用，这与上述前人的研究结果一致。但对于不同的非生物胁迫，基因的表达模式有明显差别。如：Ghir_A03G010510、Ghir_A10G012520、Ghir_D10G015270 在干旱胁迫和盐胁迫处理后，上调表达；而它们在低温胁迫、高温胁迫处理后先下调表达，在处理后期又上调表达。这说明不同的胁迫处理下，基因的表达具有显著差异性。对于GR基因在不同处理条件下的调控方式还需要进一步研究。

4 结论

在全基因组水平鉴定和分析了亚洲棉、雷蒙德氏棉、海岛棉、陆地棉中GR基因家族成员特征及潜在功能。发现GR基因家族在棉属分化前已存在2 种类型。在异源四倍体棉花中，随着亚洲棉和雷蒙德氏棉的偶然杂交和随后的全基因组加倍事件而发生家族扩张。GR基因家族成员在进化过程中自身和编码产物的结构保守，且受到纯化选择的作用。表达分析显示GR基因家族成员均参与棉花非生物胁迫应答，但其在不同类型的胁迫条件下表达模式略有差异。棉花GR基因家族的鉴定及分析可为深入解析GR基因功能提供参考，为棉花抗逆育种提供基因资源。