棉花HDAC 基因家族鉴定及其在黄萎病菌侵染下的表达分析

王艳情，郑杰，许艳超，蔡小彦，周忠丽，侯宇清，王坤波，王玉红，陈浩东，刘方*，李志坤

（1.河北农业大学/ 华北作物改良与调控国家重点实验室，河北 保定 071001；2.中国农业科学院棉花研究所/ 棉花生物学国家重点实验室，河南 安阳 455000）

棉花是重要的经济作物，其产量和品质经常受到黄萎病的影响[1]。 我国棉田中的黄萎病主要由大丽轮枝菌（Verticillium dahliae）引起[2]。 黄萎病也被称为棉花的“癌症”[3]，棉株感染大丽轮枝菌后，叶片严重失绿，形成黄色病斑，造成减产甚至绝收[4]。棉属（Gossypium）有4 个栽培棉种和49个野生棉种[5-6]。 栽培棉种对病虫害及非生物胁迫较为敏感，而野生棉种长期生长在恶劣的自然环境中， 具有丰富的遗传多样性与抗病虫害基因源[7]，是棉花遗传改良的优异资源。Dong 等[8]比较分析了瑟伯氏棉（Gossypium thurberi）、雷蒙德氏棉（G.raimondii）和三裂棉（G.trilobum）在大丽轮枝菌胁迫下的表型，发现瑟伯氏棉具有很强的黄萎病抗性。 因此，进一步挖掘野生棉，特别是瑟伯氏棉中的抗病基因对棉花抗黄萎病遗传改良具有重要意义。

组蛋白去乙酰化酶 （Histone deacetylase，HDAC）通过降低染色体的乙酰化水平从而抑制转录[9]。 植物中HDAC 家族分为3 个亚家族，分别是RPD3/HDA1（Reduced potassium dependency 3/Histone deacetylase 1）、SIR2 （Silent information regulator 2）和HD2（Histone deacetylase 2）。 RPD3/HDA1 作为该家族最大的亚家族，成员具有典型的HDAC 结构域[10]。 SIR2 又被称为Sirtuins，该亚族成员与其他HDAC 没有结构相似性，并且需要烟酰胺腺嘌呤二核苷酸（NAD）作为辅助因子[11]。 HD2 亚家族成员主要有3 个结构域：具有保守五肽序列的氮末端结构域、高电荷的中央酸性结构域和可变的碳末端结构域[12]；HD2 首先在玉米中被鉴定，仅存在于植物中[13]。

近年来，在拟南芥[14]、水稻[15]和大豆[16]基因组中分别鉴定出18、17 和28 个HDAC家族基因，并发现其在植株生长发育与响应逆境胁迫中具有重要功能。 拟南芥中AtHDA6和AtHDA19均表现为促进开花[17]，而AtHDA9抑制开花[18]。拟南芥AtHDA9过表达植株对盐胁迫更加敏感；相反，hda9突变体植株对盐胁迫的抗性增强[19]。 此外，HDAC 家族成员在病原菌防御中也起着关键作用，过表达AtHDAC19基因的拟南芥通过激活ERF1而增加病程相关蛋白（Pathogenesis-related proteins，PRs）的表达增强植株的抗性，沉默该基因则导致植株对病原菌更加敏感[20]。

棉花GhHDT4D 通过调节GhWRKY33启动子区H3K9 乙酰化水平， 从而抑制GhWRKY33的表达，在应对干旱胁迫响应中发挥作用[21]。 然而，HDAC 在棉花抗黄萎病方面的作用还未见报道。 本研究以棉花基因组数据为基础，鉴定并比较分析了陆地棉、亚洲棉和瑟伯氏棉中HDAC基因家族成员基因结构、保守结构域、进化关系及顺式作用元件的差异，并结合瑟伯氏棉转录组数据[8]和实时荧光定量聚合酶链反应（Quantitative real-time polymerase chain reaction，qRT-PCR）明确了在黄萎病菌侵染情况下GthHDAC家族基因的表达模式，研究结果为进一步挖掘野生棉中优异抗病基因提供了基础，并为棉花抗病育种提供了新的基因资源。

1 材料与方法

1.1 试验材料

所用的二倍体瑟伯氏棉（D1）由中国农业科学院棉花研究所野生棉课题组/ 国家种质三亚野生棉圃提供。

1.2 试验方法

1.2.1棉花HDAC 基因家族成员的鉴定。二倍体亚洲棉（G.arboreum，A2）和瑟伯氏棉、四倍体陆地棉（G. hirsutum，AD1）基因组数据来自中国农业科学院棉花研究所（http://cgp.genomics.org.cn/page/species/index.jsp）。 按照李晓斐等[22]的方法从Pfam[23]（http://pfam.xfam.org/）数据库中下载HDAC结构域（PF00850）、SRT 结构域（PF02146），从TAIR（http://www.arabidopsis.org）网站下载拟南芥HD2 的氨基酸序列后，以PF00850、PF02146以及HD2 的序列为种子序列，利用本地Blast 搜索3个棉种的蛋白质数据库，得到候选HDAC 家族蛋白序列。 通过在线网站Pfam、NCBI CDD 数据库（http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.sg）网站验证候选HDAC 蛋白中的保守结构域。

1.2.2基因结构、保守基序和顺式作用元件分析。HDAC家族基因的基因结构信息，包括开放阅读框（Open reading frame，ORF）长度、蛋白质长度、外显子数均从棉花功能基因数据库（https://cottonfgd.org/）获得。 使用MEME（http://memesuite.org/）在线网站分析家族成员的保守基序，并使用TBtools[24]软件对其进行可视化处理。 使用ExPaSyServer（http://web.expasy.org/computer_pi） 预测蛋白质的相对分子质量和等电点。 使用TBtools 软件获得HDAC基因起始密码子上游的2 000 bp 序列， 然后将获得的序列提交数据库PlantCARE（http://www.dna.affrc.go.jp/PLACE/signalscan.html），预测其中的顺式作用元件。

1.2.3染色体分布和亚细胞定位预测。 根据HDAC家族成员的基因ID，使用TBtools 软件在3 个棉种的基因组注释文件中进行搜索， 以确定HDAC基因的染色体位置信息并对其进行可视化处理。 使用在线工具TargetP 1.1（http://www.cds.dtu.dk/services/TargetP）对HDAC基因编码的蛋白进行亚细胞定位预测，并通过WoLFPSORT（https://wolfpsort.hgc.jp/） 网站中植物模式进行验证。

1.2.4系统发育分析。为了更好地了解HDAC基因家族成员的进化关系， 从TAIR 网站下载拟南芥HDAC基因家族的蛋白质序列，与鉴定出的陆地棉、亚洲棉和瑟伯氏棉中HDAC基因编码的氨基酸序列进行多序列比对。 利用MEGA X[25]软件，采用最大似然法（Maximum likelihood method）构建系统发育树，Bootstrap 值设为1 000，其余参数均为默认值。

1.2.5胁迫处理。 瑟伯氏棉幼苗长到两叶一心时，采用伤根法[25]接种大丽轮枝菌，取接种后0 h、24 h、48 h 的根、茎、叶于液氮速冻，放入超低温冰箱（－80 ℃）保存，用于后续试验。 每个处理设3次生物学重复。

在棉花两叶一心期分别使用5 mg·L－1乙烯（Ethylene, ET）、10 μmol·L－1茉莉酸（Jasmonic acid,JA） 和10 μmol·L－1水杨酸 （Salicylic acid,SA）均匀地喷施在叶片上，在处理后0 h、6 h、12 h和48 h 收集根、茎、叶。所有样品在液氮速冻后放入超低温冰箱（－80 ℃）保存，用于后续基因表达分析。 每个处理设3 次生物学重复。

1.2.6瑟伯氏棉HDAC 基因家族抗病表达模式分析。 基于大丽轮枝菌胁迫下瑟伯氏棉根、茎、叶的转录组数据， 将HDAC家族基因的表达量用FPKM（Fragments per kilobase of exon per million fragments mapped）值表示，经log2（1＋FPKM）标准化处理， 采用TBtools 软件绘制热图并进行归一化处理，最终得到HDAC家族成员的热图。

1.2.7RNA 提取和qRT-PCR。 采用RNAprep Pure Plant Plus Kit 提取试剂盒提取样品RNA，通过琼脂糖凝胶电泳以及分光光度计检测样品质量。使用TranScript-All-in-One First-Strand cDNA Synthesis SuperMix 反转录试剂盒获得每个样品的cDNA。 使用Primer Premier 5 设计GthHDAC家族基因的特异性引物（表1），以棉花GthUBQ7（Genebank：No. AY189972） 作为内参基因进行qRT-PCR 分析，按照2－ΔΔCt法计算基因的相对表达量，具体方法见参考文献[26]，每组数据来自3个生物学重复以及3 次技术重复。

2 结果与分析

2.1 棉花HDAC 家族基因的鉴定

在陆地棉、亚洲棉以及瑟伯氏棉中分别鉴定出30、15 和13 个HDAC基因， 按照其在染色体上的位置分别命名为GhHDAC1～GhHDAC30、GaHDAC1～GaHDAC15、GthHDAC1～GthHDAC13（附表1～3）。

理化性质分析显示，3 个棉种HDAC家族成员的ORF 序列差异较大，长度从480 到2 100 bp，大多数为1 000～2 000 bp。陆地棉HDAC基因的平均外显子长度 （212.62 bp） 明显大于亚洲棉（145.88 bp）和瑟伯氏棉（138.21 bp）HDAC基因（图1A）； 亚洲棉HDAC成员的平均内含子长度大于陆地棉和瑟伯氏棉HDAC家族基因 （图1B）；瑟伯氏棉HDAC家族成员的外显子数量与亚洲棉和陆地棉HDAC 家族成员的外显子数量差异较大，而亚洲棉和陆地棉间无明显差异（图1C）；3 个棉种的HDAC基因ORF 的GC 含量差异较小，在40%～48%（图1D）。 亚细胞定位预测结果显示，3 个棉种大多数HDAC 蛋白定位在细胞核中，部分定位在细胞质，个别定位于线粒体和叶绿体（附表1～3）。

图1 3 个棉种中HDAC 家族基因结构分析Fig. 1 Structural analysis of HDAC genes in three Gossypium species

表1 qRT-PCR 所用引物序列Table 1 Primer sequences used in qRT-PCR

2.2 棉花HDAC 家族基因的染色体分布

根据亚洲棉、瑟伯氏棉以及陆地棉的基因组注释信息， 对58 个HDAC基因进行了染色体分布分析。由结果（图2）发现，亚洲棉HDAC家族的15 个成员分布在9 条染色体上，其中在Chr13 染色体上有较多的分布，分别是GaHDAC13、GaHDAC14和GaHDAC15。 陆地棉30 个GhHDAC基因分布在17 条染色体上，每条染色体上有1～3 个GhHDAC基因；瑟伯氏棉的13 个GthHDAC基因分布在8 条染色体上，每条染色体上有1～3 个基因， 其中13 号染色体上最多， 分别是GthHDAC11、GthHDAC12和GthHDAC13。

图2 3 个棉种HDAC 基因的染色体分布Fig. 2 Chromosome distribution of HDAC genes in three Gosspium species

2.3 HDAC 基因的系统发育分析

为了探究HDAC基因的进化规律，将3 个棉种中的HDAC 蛋白与拟南芥HDAC 蛋白共同构建系统发育树（图3），结果显示，3 个棉种和拟南芥中的HDAC 蛋白在RPD3/HDA1、HD2 和SIR2 亚族均有分布， 表明HDAC 蛋白在棉花和拟南芥不同物种间具有保守性。

图3 陆地棉、亚洲棉、瑟伯氏棉和拟南芥中HDAC 蛋白的系统发育分析Fig. 3 Phylogenetic analysis of HDAC proteins in three Gosspium species and Arabidopsis thaliana

2.4 棉花HDAC 蛋白质结构域及其保守基序分析

保守基序及结构域分析（图4）显示，SIR2 中的保守基序为motif3、motif9；RPD3/HDA1 的保守基序最多， 分别是motif1、motif2、motif4、motif5、motif6、motif7、motif8 和motif10；HD2 中的成员仅含有motif10。 由此可见，同一亚族成员含有相同的保守基序， 推测它们可能具有相似功能，不同亚组间的保守基序差异较大，推测在不同亚族的HDAC 蛋白可能存在功能差异。

图4 3 个棉种HDAC 的结构分析Fig. 4 Structural analysis of HDAC proteins in three Gossipium species

2.5 顺式作用元件分析

顺式作用元件预测发现，HDAC基因的启动子区域有大量与植物生长发育相关的顺式作用元件 （图5）， 如生长素反应元件（AuxRR-core，TGA-element）、分生组织表达相关的元件、胚乳表达相关的元件、MYB 反应元件。 其中，MYB 反应元件含量最多，GaHDAC基因中有120 个，GhHDAC基因中有217 个，GthHDAC基因中有125 个。 此外，在棉花HDAC基因的启动子中还发现了一些与激素反应相关的元件，例如脱落酸反应元件（ABRE）、茉莉酸甲酯反应元件（CGTCA-motif 和TGACG-motif）、SA 反应元件（TCA-element）。还发现了一些与胁迫相关的顺式作用元件，分别是干旱响应元件（MBS）、低温响应元件（LTR）、应急应答元件（TC-rich repeats）以及一些响应病原体的顺式作用元件（W-box，G-box，as-1 element 和P-box）等。

图5 棉花HDAC 基因的顺式作用元件分析Fig. 5 Cis-acting elements analysis of HDAC genes in cotton

2.6 瑟伯氏棉HDAC 基因家族在大丽轮枝菌胁迫下的表达模式分析

转录组数据分析表明， 瑟伯氏棉中13 个GthHDAC基因在黄萎病菌胁迫下具有不同的时空表达模式（图6）。黄萎病菌胁迫后12 h，大部分GthHDAC基因在叶片中差异表达，其中GthHDAC1、GthHDAC2、GthHDAC6、GthHDAC8、GthHDAC9、GthHDAC10、GthHDAC12和GthHDAC13上调表达，GthHDAC3、GthHDAC4、GthHDAC5和GthHDAC7下调表达， 说明GthHDAC基因可能在棉花抵御黄萎病菌中发挥作用。 有趣的是， 在黄萎病胁迫下，GthHDAC11在根和茎中高水平表达；而在黄萎病菌胁迫后48 h，在叶片中出现一定程度的上调表达，主要在茎和叶中表达。 在黄萎病菌胁迫后12 h，GthHDAC7在茎和叶中下调表达；48 h 时其在茎中表达量提升，但在叶片中依然下调表达。

图6 黄萎病菌胁迫下GthHDAC 基因转录组分析Fig. 6 Transcriptome analysis of GthHDAC genes under V. dahliae stress

利用qRT-PCR 技术进一步分析GthHDAC基因在黄萎病菌侵染后0 h、24 h 和48 h 的表达模式（图7）。 黄萎病菌侵染后，几乎所有GthHDAC基因在叶片中上调表达，GthHDAC4、GthHDAC7、GthHDAC11和GthHDAC12在侵染后24 h 和48 h 上调表达倍数均较高。 黄萎病菌侵染后24 h，GthHDAC7和GthHDAC13在根中显著下调表达；GthHDAC4、GthHDAC7和GthHDAC10在茎中上调表达，GthHDAC1、GthHDAC2、Gth-HDAC3、GthHDAC5、GthHDAC8、GthHDAC11、GthHDAC12和GthHDAC13在茎中下调表达。其中：GthHDAC4在黄萎病菌侵染后24 h 在根中下调表达，48 h 时上调表达； 黄萎病菌侵染后GthHDAC7在根中下调表达，在叶片上调表达。

图7 黄萎病菌胁迫下GthHDAC 基因在棉花不同组织的表达分析Fig. 7 Expression analysis of GthHDAC genes in different cotton tissues under V. dahliae stress

从响应大丽轮枝菌侵染的瑟伯氏棉HDAC基因家族成员中随机挑选GthHDAC4、GthHDAC5、GthHDAC11和GthHDAC13四个基因，分析其在不同激素处理下的表达模式 （图8）。GthHDAC4基因在ET 处理下表达水平逐渐上升， 在处理后48 h 达到峰值；SA 处理后48 h 表达量最高；JA 处理后12 h 达到峰值， 超过对照20 多倍。GthHDAC5基因的表达量在ET 处理后6 h 达到峰值；SA 处理后48 h，GthHDAC5表达量逐渐达到最高；JA 处理后12 h，GthHDAC5的表达量达到峰值。GthHDAC11和GthHDAC13基因的表达量在ET、SA 和JA 处理后6 h 达到峰值， 随后降低。 上述结果表明， 瑟伯氏棉中GthHDAC基因对黄萎病菌侵染产生应答， 并且其行使功能可能与ET、SA 和JA 等信号通路有关。

图8 不同激素处理下GthHDAC 基因的表达Fig. 8 Expression of GthHDAC genes under different hormone treatments

3 讨论

组蛋白乙酰化修饰通过调节相关基因的转录在植物发育和抵御病原菌侵染中发挥重要作用，而HDAC 介导的去乙酰化是组蛋白乙酰化修饰中重要的一环[27]。 四倍体陆地棉（AADD）是由A 基因组供体棉种与D 基因组供体棉种通过远缘杂交形成的，亚洲棉属于A 基因组棉种，具有较强的黄萎病抗性[28-32]。 瑟伯氏棉是D 基因组棉种中的抗黄萎病的棉种[26]。 因此，研究亚洲棉、瑟伯氏棉和陆地棉HDAC基因家族的系统发育关系，有利于加深对HDAC基因家族的进化及其在黄萎病抗性中的作用的认识。

本研究在陆地棉、亚洲棉以及瑟伯氏棉中分别鉴定出30、15 和13 个HDAC基因。GthHDAC基因成员数量与GaHDAC基因的成员数量之和近似于GhHDAC家族基因成员数量， 这一结果进一步印证了Beasley[33]关于四倍体种结合了A染色体组和D 染色体组的遗传物质这一观点。系统发育分析显示，3 个棉种和拟南芥中的HDAC基因亲缘关系较近，表明HDAC基因家族在进化上高度保守。

HDAC 通过影响植物激素合成、信号转导等过程中相关基因的组蛋白乙酰化水平，调控相关基因的表达量，从而调节植物的生长发育、逆境响应等过程[22]。本研究发现在HDAC基因的启动子区域富含多种激素反应元件， 如SA 胁迫响应元件、JA 胁迫响应元件、 脱落酸胁迫响应元件等。 拟南芥AtHDA6 通过与F-box 元件相互作用进而参与到JA 信号通路， 在植株抵抗病原菌过程中发挥作用[17]；AtHDA19基因可以被ET 和JA 诱导， 过表达该基因可以诱导ERF1以及JA调控的病程相关蛋白基因的表达，进而增强拟南芥对芸薹生链格孢菌（Alternaria brassicicola）的抗性[20]。 本研究中，在ET、SA 和JA 的处理下，GthHADC4、GthHADC5、GthHADC11和GthHADC13这4 个基因均不同程度上调表达，这些结果表明GthHDAC基因可能参与了ET/JA和SA 信号通路，进而参与到棉花黄萎病抗性反应。

Brosch 等[34]发现玉米与圆斑病菌之间的相互作用可以抑制寄主HDAC 的活性，即HC 毒素通过抑制HDAC的表达充当诱导子，激活防御反应。 本研究利用转录组数据分析和qRT-PCR 分析验证了瑟伯氏棉HDAC基因家族成员在黄萎病菌胁迫下的表达谱。 但在2 种分析方式下HDAC基因的表达量存在一定差异，可能是因为两者表达量的计算方式不同。 其中GthHDA1、GthHDAC3、GthHDAC11和GthHDAC12基因显著上调表达，这些基因大多属于RPD3/HDA1 亚家族。 在拟南芥中，RPD3/HDA1 型HDAC基因（HDA19和HDA6）可能在ET 和JA 信号转导和病原体反应中发挥重要作用[17]。此外，邓琴霖等[35]发现RPD3/HDA1 家族成员可调控植物逆境应答， 推测GthHDAC家族基因中RPD3/HDA1 亚组成员可能参与植物对黄萎病的防御。

4 结论

通过全基因组分析，在陆地棉、亚洲棉和瑟伯氏棉中分别鉴定出30、15 和13 个HDAC家族基因。 通过进化关系分析，将其分为3 个亚家族。亚细胞定位预测发现，大部分HDAC 定位于细胞核。 表达模式分析发现，在黄萎病菌胁迫下，大多数瑟伯氏棉GthHDAC基因在叶片中上调表达，并且受ET、SA 和JA 的诱导。 本研究通过解析陆地棉、亚洲棉和瑟伯氏棉中HDAC基因的进化规律，明确了不同激素处理下HDAC基因在瑟伯氏棉中的表达模式，为研究HDAC 在抵御棉花黄萎病中的作用机理奠定了基础。

附表：

详见本刊网站（http://journal.cricaas.com.cn）本文网页版。

附表1 亚洲棉HDAC家族的理化性质

Table S1 Physicochemical properties ofHDACgenes inG.arboreum

附表2 陆地棉HDAC家族的理化性质

Table S2 Physicochemical properties ofHDACgenes inG.hirsutum

附表3 瑟伯氏棉HDAC的理化性质

Table S3 Physicochemical properties ofHDACgenes inG.thurberi