基于表观遗传学探索玉米双向启动子对转录的影响

朱倩洁，曹亚欣，马宇馨，张林旺，邢继红,2，张 康,2，董金皋,2

(1.河北农业大学 生命科学学院，河北 保定 071000；2.河北省植物生理与分子病理学重点实验室/河北农业大学真菌毒素与植物分子病理学实验室，河北 保定 071000；3.河北省荣军医院，河北 保定 071000)

双向启动子 (Bidirectional promoter，BDP) 是位于两个相邻且转录方向相反的基因之间，能够同时启动其两端结构基因转录的一类特殊启动子。具有双向启动子的基因普遍存在与酵母［1］、果蝇［2］、人类［3-4］、植物［5］等真核生物基因组中。研究发现，具有双向启动子的基因主要参与人的DNA修复、细胞周期、物质代谢和人类疾病等过程［6-10］。与单向启动子相比，双向启动子可以结合更丰富的RNA聚合酶II，具有更丰富的H3、H3K9和H3K27的乙酰化作用、H3K4me 2/3的甲基化作用［11-12］，而H4乙酰化作用较少，这表明双向启动子可能具有独特的染色质特征。随着全基因组测序和植物转录数据的发布，植物双向启动子已经受到了相当大的关注。到目前为止，已经开展了拟南芥［13］、水稻［14］、玉米［15］和杨树［5］中双向启动子的研究。有研究发现哺乳动物和植物基因组［14,16］中，双向启动子具有保守的序列特征。然而，目前还不清楚植物中基因对的双向转录和共表达的表观遗传调控机制。

目前的研究发现ATAC-seq (Assay for Transposase-Accessible Chromatin using Sequencing) 数据与基因表达水平呈正相关，其所代表的开放染色质为转录因子或增强子等调控因子的结合部位，主要富集在转录起始位点 (Transcription start site，TSS) 的上游。组蛋白修饰则主要富集与TSS下游或基因区域，与基因表达水平之间的关系较为复杂，例如H3K4me3、H3K9ac、H3K27ac与基因的表达水平呈正相关，主要富集于TSS下游；H3K36me3与基因表达水平呈正相关，主要富集于整个基因区域；H3K27me3与基因表达水平呈负相关，主要富集于整个基因区。通过开放染色质及不同组蛋白修饰的富集位置及其与基因表达水平的关系可以为双向启动子的研究提供很好的研究基础。

目前，玉米中关于表观遗传修饰对双向启动子及相关基因表达模式的调控作用尚不清楚。本研究拟利用公共数据平台中的数据，对易接近染色质区域(Accessible chromatin region，ACR)和组蛋白修饰对玉米双向启动子相关基因转录的影响进行分析，明确玉米中表观遗传调控与双向启动子及相关基因之间的关系，为进一步深入研究双向启动子对双向转录基因对的调控功能及其机制奠定基础。

1 材料与方法

1.1 玉米双向启动子的识别

玉米基因组相关信息从玉米基因组注释数据库(MaizeGDB) 中获取。确定具有双向启动子的基因对中的双向定向，每个基因对TSS之间的基因间区域长度为0～1 kb的双向启动子 (BDPs)。BDPs分为三类：0～250 bp (BDPs I)、250～500 bp (BDPs II) 和500～1000 bp (BDPs III)。所有被标注为蛋白质编码基因的基因对都用于下游分析。相比之下，单向启动子 (Unidirectional promoters, UDPs) 是从单向基因中挑选出来的，其表达水平类似于具有双向启动子的基因对，用于与BDPs进行并行分析。利用ATAC-seq数据，对具有BDPs的ACRs进行分析，获得ACRs的位置信息。根据ACRs在BDPs中的分布情况，将所有的玉米BDPs分为四个类别：one mid-ACR，只显示在BDPs中间的一个峰；bi-ACRs，表示位于BDPs中的两个单独或部分重叠的峰；one amesial ACR，只表示在BDPs中不对称的一个峰；和no ACR，表示在BDPs中没有发现峰。

1.2 数据收集及整理

收集公共数据平台Gene Expression Omnibus(GEO) 和Sequence Read Archive (SRA)数据库中玉米相关的转录组及表观基因组数据，并进行整理和挑选。

1.3 数据分析

对于公共数据平台下载的数据，首先利用SRA ToolKit将原始数据从SRA格式转化为fastq格式，然后对数据进行质量检测和去接头处理获取高质量的序列。对于ChIP-seq（Chromatin immunoprecipitation sequencing） 和ATAC-seq数据，利用Bowtie2将其比对到基因组上。随后，利用MACS软件对ChIP-seq和ATAC-seq数据中组蛋白修饰或易接近染色质的富集位置进行预测。在得到组蛋白修饰或开放染色质的富集位置后，利用软件对这些富集的位置进行下一步分析。

利用TopHat软件将RNA-seq数据中的序列比对到基因组上，然后用Cufflinks软件计算基因表达FPKM (Fragments Per Kilobase of exon per Million fragments mapped reads)。

1.4 数据可视化

利用Integrative Genomics Viewer基因组浏览器将ATAC-seq、ChIP-seq及RNA-seq数据上传至IGV进行可视化，观察ATAC-seq和ChIP-seq数据富集位置与双向启动子的关系，进而分析基因表达水平与双向启动子的关系。

1.5 共表达分析

玉米相关基因之间的皮尔逊相关系数来自于玉米共表达网络数据库 (MCENet)［17］。皮尔逊相关系数计算两两基因间表达的相关性。使用Perl语言中的Statistics::Basic qw (:all) 模块的correlation函数，依据两两基因的FPKM值进行皮尔逊相关系数的计算。

1.6 显著性检验

为了确定BDP和UDP之间的基因表达，组蛋白标记和核小体占有率是否显著差异，进行了双样本student t检验(Student's t-test)。双样本t检验是检验2个样本平均数与其各自所代表的总体的差异是否显著，一般以P＜0.05为有统计学差异，P＜0.01为有显著统计学差异，P＜0.001为有极其显著的统计学差异。

2 结果与分析

2.1 易接近染色质区域在玉米BDPs中的分布

本研究对玉米基因组中有注释的基因进行分析，筛选到BDPs I、BDPs II和BDPs III类型的双向启动子分别有175、182和309个基因对。利用公共数据平台中玉米的ATAC-seq数据，对玉米双向启动子中易接近染色质区域 (Accessible Chromatin Region，ACR)的分布情况进行分析。结果发现，ACR在BDPs I、BDPs II和BDPs III类型的双向启动子中的分布存在明显的差别(表1)。其中，BDPs I类型的双向启动子中，ACR主要位于双向启动子的中间（one mid-ACR），有86个基因对，达到49.14%；其次是no ACR类型，有34个基因对，达到19.43%；bi-ACRs类型仅有13个基因对，为7.43%。BDPs II类型的双向启动子中，no ACR类型为主要形式，有54个基因对，为29.67%；其次是one mid-ACR类型，有43个基因对，为23.63%；bi-ACRs类型所占比例较少，有25个基因对，为13.73%。BDPs III类型的双向启动子中，no ACR类型也为主要形式，有100个基因对，达到32.36%；其次是bi-ACRs类型，有45个基因对，为14.56%；one mid-ACR类型最少，有32个基因，为10.36%。

表1 BDPs中ACRs的分布Table 1 Distribution of ACRs within BDPs

2.2 ACRs对具有双向启动子基因表达的影响

对BDPs与相应基因对的表达水平之间的关系进行分析，发现具有BDPs I型启动子的基因对的表达水平明显高于其它两种启动子类型的基因和随机选择的单向基因的表达水平，BDPs II中基因对的表达水平也明显高于BDPs III基因对的表达水平，BDPs III和UDPs之间没有显著的区别(图1a)。对ACR在双向启动子中的分布与基因对的表达水平之间的关系进行分析，发现one mid-ACR类型双向启动子的基因的表达水平最高，其次是bi-ACRs和one amesial ACR类型双向启动子的基因，而no ACR类型双向启动子的基因的表达水平最低(图1b)。

图1 不同类型双向启动子相关基因对的表达水平Fig.1 Expression levels of gene pairs associated with bidirectional promoters of different types

此外，在含有 one amesial ACR的BDPs中，靠近ACR的基因的表达水平明显高于ACR的远端(p值=4.4e-32，t检验)；在包含one mid-ACR基因对(p值=1.9e-9，t检验)和包含bi-ACRs基因对(p值=3.7e-5，t检验)中靠近ACR的基因的表达水平明显高于ACR的远端，虽然3种类型的BDPs都具有显著差异，但是含有 one amesial ACR的BDPs的差异仍明显高于其他两种(图2)。这样的结果说明BDPs的不同长度影响相关基因对的表达；相对于双向基因的TSS，ACRs的物理位置可能会影响具有双向启动子的基因对的表达。

图2 不同类型双向启动子相关基因对中差异表达水平Fig.2 Differential expression levels of different types of bidirectional promoter related genes

2.3 ACRs对具有双向启动子的基因对的共同表达的影响

ACRs的物理位置可能会影响BDPs的共表达方式，为了验证这种现象，本研究对ACRs物理分布和BDPs共同表达系数之间进行了相关性分析。结果表明，具有One mid-ACR的基因对共表达程度最高 (图3)，其中one mid-ACR基因对和one amesial ACR基因对共表达程度具有显著差异 (p值=0.01，t检验)。One mid-ACR 基因对和bi-ACRs 基因对 (p值=0.18，t检验)，bi-ACRs基因对和one amesial ACR 基因对 (p值=0.49，t检验) 不具有显著差异。这些结果表明，ACRs的物理位置可能会影响具有双向启动子的基因对的共同表达。

图3 ACRs对BDPs共同表达的影响.Fig.3 Effect of the ACRs profile on the co-expression of bidirectional promoter gene pairs

2.4 组蛋白修饰在双向启动子附近的分布影响基因的表达

为探索ACRs与基因表达之间的关系, 本研究利用公共数据平台中的RNA-seq 数据计算了玉米所有基因的表达值 (FPKM)，并按照表达水平将所有基因分成了六组, 其中 FPKM 大于零的基因被平均分成了五组 (表达水平最高的定义为第一组, 反之表达最低的定义为第五组， FPKM 等于零的基因统一定义为第六组)。然后，计算这六组基因周围ACRs的富集情况(图4)，发现ACRs在TSS上游的+1核小体位置的分布与基因表达水平呈正相关，并且ACRs富集水平最高点距离TSS越近，基因的表达水平越高(图4黑点表示)。

图4 ACRs与基因表达的关系Fig.4 Relationship between ACRs and gene expression

组蛋白修饰在真核生物基因表达过程中发挥非常重要的作用。为了分析双向启动子附近的组蛋白修饰富集情况，对相同组织的组蛋白修饰数据进行了分析，其中包括3个与转录激活相关的组蛋白修饰(H3K4me3、H3K27ac和H3K9ac) 和一个与转录抑制相关的组蛋白修饰 (H3K27me3)。在分析结果中可以发现，ACRs在双向启动子中具有明显的富集，同时，与转录激活相关的组蛋白修饰H3K4me3，H3K27ac和H3K9ac在双向启动子两侧具有明显的富集，而转录抑制相关组蛋白修饰H3K27me3在双向启动子及其附近区域的分布没有明显的富集现象 (图5)。

图5 BDPs附近ACRs和不同类型组蛋白修饰的富集水平Fig.5 Enrichment of ACRs and different types of histone modifications near BDPs

为探索双向启动子及其附近区域ACRs所处的物理位置与组蛋白修饰富集水平之间的关系，对不同类型的双向启动子附近的组蛋白修饰情况进行了分析，发现在one mid-ACR和bi-ACRs附近H3K4me3，H3K27ac和H3K9ac显著的富集在双向启动子的两端的TSS下游，但是one amesial-ACR (3′端TSS) 双向启动子附近，H3K4me3、H3K27ac和H3K9ac显著的与ACR共同富集，而没有no-ACR一端 (5′端TSS) 的 TSS下游H3K4me3、H3K27ac和 H3K9ac的富集水平显著降低 (图6)。H3K27me3则呈现出与H3K4me3、H3K27ac和H3K9ac截然相反的趋势。这些结果说明双向启动子对下游基因的调控作用需要开放染色质与多种组蛋白修饰构成的染色质状态共同发挥作用。

图6 不同类型BDPs附近组蛋白修饰的富集水平Fig.6 Enrichment of histone modifications near different types of BDPs

3 讨论与结论

易接近染色质区域的敏感性与真核生物基因组中单向基因的表达水平直接相关［18,19］。然而，ACRs与具有双向启动子的基因对表达的关系仍不清楚。在本研究中发现玉米基因组中含有one mid-ACR或bi-ACRs双向启动子的基因对具有明显的共同表达现象，这表明在玉米双向启动子中ACRs会影响具有双向启动子的基因对的表达模式。在启动子区域内的ACRs的对称位置可能是具有双向启动子的基因对共调节的关键参与者。在启动子区域中，ACRs通常为RNA聚合酶II和其他转录因子的结合区域，从而参与到基因转录的调控过程中［20］。在玉米BDPs中，ACRs (one mid-ACR或bi-ACRs)的对称分布在基因对的形成中扮演了2个可能的角色。其中一个作用是，ACRs的存在代表了开放染色质区域，这可能同时促进基因对的表达。另一个作用是由于共享相同的调控元件，基因对的表达可以由相同的转录机制控制，具有双向相等的效率。

在对玉米双向启动子进行鉴定时，通过与水稻双向启动子相关文献［14］进行比较后发现，由于玉米基因组要明显大于水稻基因组，并且玉米基因组中相邻基因之间的距离显著大于水稻相邻基因之间的距离。因此，在玉米基因组中鉴定到小于1000 bp的双向启动子数量也明显少于水稻。另外，还发现玉米中mid-ACR或bi-ACRs相关的基因对表达水平存在明显差异，而之前报道的水稻mid-ACR或bi-ACRs相关基因对之间的表达水平则没有明显的差异。猜测可能是基因组之间的差异性导致这种结果的出现，当然不排除可能是分析软件的差异对这种结果的出现造成了一定的影响［21］。

在真核生物中已经深入研究了核小体定位参与基因表达或基因调控的进化［22-24］。然而，关于染色质组织对植物中共表达基因对的调节的影响知之甚少。组蛋白修饰影响转录因子在ACR中的结合，染色质结构在调节哺乳动物中聚集基因的表达中起关键作用，有人提出基因转录可以在起始阶段或在延长过程中进行调控。两个步骤都可以受到在启动子和基因体区域中的组蛋白修饰的影响。因此，在转录起始步骤 (H3K4me3和H3K9/K27ac) 的活性或者在延伸步骤 (H3K27me3) 可以协调停滞或延长的RNA聚合酶II。BDPs中含有ACR，并且两侧核小体中包含激活相关的组蛋白修饰 (H3K4me3、H3K27me3、H3K9/27ac)，才能促进了相应的基因对的表达，表明ACRs的物理位置和周围的组蛋白修饰在玉米BDPs发挥功能的过程中起着重要作用。

结合与one mid-ACR和bi-ACRs相关基因的共表达现象，得出结论，转录激活相关的组蛋白修饰可能产生有利于基因对共表达的染色质结构。另一方面，发现ACR与TSS越接近，基因的表达水平越高。综上所述，本研究结果表明，ACRs的物理位置和组蛋白修饰的富集水平在基因对的共表达过程中起着重要的作用，为进一步探索双向启动子在玉米分子育种及抗病过程的功能研究奠定理论基础。