一作深度解读-张一婧课题组合作完成面包小麦精细表观组图谱绘制及全基因组顺式作用元件鉴定
2019年7月15日，Genome Biology期刊在线发表中科院分子植物科学卓越创新中心/植物生理生态研究所张一婧研究组与南京农业大学张文利研究组合作完成的题为“The bread wheat epigenomic map reveals distinct chromatin architectural and evolutionary features of functional genetic elements”的研究论文。该工作生成并绘制面包小麦精细的表观组图谱，以此为基础针对性地开发整合计算流程，对全基因组顺式调控元件进行了系统的挖掘与鉴定,并初步探索了其作用机制,为小麦基因调控机制解析研究提供了重要的资源。
广泛种植的六倍体面包小麦（T. aestivum, 2n = 6x= 42, BBAADD）具有庞大而复杂的基因组，高质量的面包小麦全基因组序列于2019年初公布，大小约为16 Gb，是人类基因组的5倍。其中93%是非编码序列，蕴含着丰富的基因远端调控元件，在小麦全基因组水平准确鉴定顺式元件并解析其调控机制，是研究小麦多倍化及驯化过程中基因表达调控的关键步骤。由于表观修饰在基因调控过程中发挥了重要作用，有机整合表观组信息有助于在全基因组水平精准预测顺式调控区域。但是，与基因组序列相对简单的模式生物相比，庞大而复杂的小麦基因组对组学数据产生、数据分析及机制解析均带来巨大挑战。
1普通小麦基因组的染色质状态图谱
为了系统分析普通小麦中的表观基因组特征，合作团队生成并分析了以DNaseI过敏位点(DHS)为指征的染色质开放区域，基于亚硫酸盐测序的甲基化组数据，以及七种组蛋白修饰的ChIP-seq数据（图1a）。在所有的染色质标记中，H3K9me2和DNA甲基化主要分布在着丝粒近端区域，其他与激活功能相关的标记主要分布在染色体两端富含基因的区域（图1a）。我们定义了染色质标记的富集区(peak)和甲基化区域，发现超过一半的peak以及甲基化区域的都在基因间区 (图1b)，这一比例远远高于之前在拟南芥和水稻等具有较小基因组的模式植物的研究结果。同时这个结果也暗示着这些标记可能涉及基因活性的远程调控。
图1. 染色质特征图谱分析揭示了基因的表观遗传调控
a）Circos图展示了表观遗传标记在染色体上的分布，最外层代表每条染色体，第二层表示基因密度，红色和白色分表代表密度的最高值和最低值，中间9层分别表示七种组蛋白修饰的强度以及DHS和DNA甲基化水平，最内三层表示小麦中三种主要TE的密度分布。
b)表观遗传标记peak在基因组不同区域中的分布。TSS：转录起始位点；TES：转录终止位点
c)根据基因及启动子中组蛋白修饰强度分成的五组基因。采用各个修饰在基因中标准化的强度进行K均值聚类。
d)小提琴图展示每个分组中基因的表达强度。
e)箱线图展示每个分组中基因在各个组织中的表达特异性（CV）。
f)每个分组中基因的功能域的富集分析，根据富集的P-value进行排序。
g)保守基因（旧基因）和非保守基因（新基因）与五组基因的富集分析。
2 基因周围的染色质状态
我们首先分析基因和启动子区域的染色质特征并结合转录组数据分析染色质特征与基因转录活性的关系。根据组蛋白修饰强度将基因分成的五组（图1c）。发现高表达基因倾向于被H3K4me3、H3K9ac和H3K36me3标记（图1d），与已报道的动物和植物中的情况类似。第4组基因被H3K27me3标记，该组基因显示出较高的组织特异性表达（图1e），这与PcG家族蛋白负责催化H3K27me3来调控植物的发育的研究相符合。通过功能富集分析（图1f），我们发现在第四组基因中有30%的基因属于NB-LRR基因，暗示PcG可能具有调控小麦免疫能力的功能。
小麦多倍化过程造成了很多基因家族的扩张，同时也产生了大量的新基因。我们希望探究这些新产生的基因在表观层面如何被调控。将小麦的基因分成两个类别：将与二倍体和四倍体祖先种共线性区域高度相似的基因定义为“旧”基因即保守基因，将与祖先种相似性低或没有相似性区域的基因定义为“新”基因即非保守基因，观察到“新”基因与H3K27me3显著富集（图1g）。加上H3K27me3主要在染色体两端基因密度高的区域富集，推测H3K27me3很可能抑制特定基因组区域内的新基因。
3 染色质状态在亚基因组之间的动态变化
普通小麦是具有三个亚基因组的异源六倍体植物，亚基因组之间的多样性增加了小麦对环境的适应性。为了鉴定亚基因组中表观修饰的偏好性以及这种偏好性与基因表达调控的关系，我们采用三元图的方法对小麦中的三联组基因(亚基因组中基因个数比为1:1:1的同源基因组)进行分析，表观修饰的偏好性可以分为七类，包括一个在三个基因中平衡的类别和六个分别在某个基因中更高或更低修饰的类别（图2a）。在所有的组蛋白修饰中，H3K4me3，H3K9ac，H3K36me3的多样性更高，其他的修饰相对在亚基因组中相对平衡（图2b）。将亚基因组偏好的表观修饰与亚基因组的偏好表达做相关分析，我们发现H3K4me3，H3K9ac，H3K36me3的修饰偏好性与基因表达的偏好性高度相关（图2c）。
图2. 同源基因启动子中各种表观标记在亚基因组之间修饰偏好
a）三元图展示七种表观修饰标记在同源基因组中的相对强度。每个点代表一组三联同源基因，基于各个亚基因组的标准化的读数深度与所有亚基因组的读数深度的比值确定每个点的坐标。
b)同源基因组中具有不同修饰偏好性的比例分布。
c)表观修饰对亚基因组的结合的偏好性与亚基因组的偏好表达做富集分析，深蓝色代表显著富集。
4 根据染色质状态特征预测基因组中功能元件
普通小麦具有很大的基因间区，是人类基因组的5倍，这些间区序列中的很多序列被多种表观修饰所标记，很可能富集了基因远端调控元件。为了总结这些染色质标记在全基因组中组合模式，我们应用多变量隐马尔可夫模型（HMM）来进行染色质状态的分类，将所有的组蛋白修饰状态的组合分成了15种类别，即15种染色质状态（State 1 – State 15），每种状态有其对应的特点与功能。我们针对每种染色质状态进行了基因结构分布、染色质可及性、与祖先种相比的序列保守性、CpG 岛的分布、DNA甲基化程度的分析（图3）。在State 15中观察到所有修饰的超高覆盖，但其序列的保守性比较低，因此，判断该状态中的修饰信号更可能是一种背景噪音。
在State 1 - State4中的区域富集编码基因序列，占整个基因组的1.7%。这些States主要被H3K4me1和H3K36me3修饰，是两种典型的标记活跃转录基因的组蛋白修饰（图3a和3b）。在该区域中91%的序列中都检测到了转录活性，其中65%是已经注释的基因，24%是未注释序列但具有mRNA的转录，2%是未注释序列但有lnRNA的转录（图3c）。图3d中RNA-seq密度分布直观地展示了该区域序列的转录活性。综上所述，这些组蛋白修饰的组合能够预测转录活跃基因，较高的核酸酶敏感性、具有H3K9ac或H3K27ac修饰是基因组调节元件如启动子和增强子的典型特征。序列占全基因组1.5% 的 State 5 – State 7 中都富集了H3K9ac，部分States富集了H3K27ac（图3a），一些区域也具有高密度的H3K4me3和H3K36me3修饰。这些状态也有较高的染色质开放程度（图3e），表明其中富含能够被转录因子和其他调节蛋白结合的DNA区域。除了具有高覆盖的DHS和组蛋白乙酰化外，State 5 – State 7序列的保守程度相对较高，与基因区域相似（图3f），甲基化程度相对比较低（图3g）,这与染色质State 1 – State 4的高甲基化程度形成了鲜明对比。有趣的是，这些区域并具有较高密度的CpG 岛覆盖（图3h），动物中的研究表明CpG岛是调控元件比较富集的区域。以上分析结果都暗示State 5 – State 7包含了大量活跃的顺式元件。
图3. 功能元件染色质特征的分析与预测
a) 用多变量隐马尔可夫模型来确定染色质状态。热图中颜色表示不同染色质状态的发射参数，颜色越深，表示该区域某种标记的强度约强。每行代表一个状态，每列代表一个染色质标记，最后一列代表给定状态的在基因组中比例。重复数据表明结果有很高的一致性。
b) 条形图显示每个State中各种基因组结构的分布。
c）饼图显示State1-4中序列转录情况的分布。
d）基于State1中的染色质特征可以预测未注释的基因。
e-f）对于每种染色质状态，统计以DHS密度为特征的开放染色质区域的分布（e）和序列保守程度（f）。
g）每种染色质状态CG、CHG、CHH三种DNA甲基化比例的分布。
h) 每种染色质状态与CpG岛重叠区域的比例（g）
i）每种染色质状态与各种类型的TE的重叠区域的比例。
j）State12和State13中的TE在染色体远端（R1和R3）、着丝粒和近端（R2和C）区域的分布。
k）State12和State13中TE与最近基因的累积距离。
5 顺式元件周围表观修饰的亚基因组保守性
为了进一步评估基因组表观遗传水平上亚基因组之间的保守性和多样性，我们对共线性区域中亚基因组之间不同染色质状态的相似性进行了打分。除State14和State15外，染色质状态在三个亚基因组中都具有一定的相似性（图4），主要包含调节元件的State5和主要包含编码基因的State1、State2在三个亚基因组中高度相似，Jaccard相似性指数范围为0.64-0.67。这表明在这些区域最有可能在功能上保守。
图4 顺式元件周围表观修饰的亚基因组保守性
亚基因组之间染色质状态相似性（Jaccard指数）的一对一的比较
Genomebrowser track 展示同源顺式元件的表观修饰状态在亚基因组之间的保守性
6 启动子和增强子的染色质特征区分
State5中包括基因近端和远端调节元件。为了分析近端启动子和远端调节元件之间的差异，我们比较了这两种类型序列在State5中的序列特征和染色质特征。在人类中，增强子通常以高丰度的H3K4me1和H3K4me3修饰为特征。然而，小麦基因编码区富集H3K4me1修饰，在调节区域State5中更富集H3K4me3 (图5a-b)，在水稻中观察到类似的模式，可以推测H3K4修饰的调节作用似乎在植物和动物之间独立进化。近端和远端功能元件都具有相似水平的染色质开放性和H3K9ac修饰水平（图5c-d）。为了进一步区分启动子和类增强子序列，我们对顺式元件做了富集分析（图5e）。不同类别的转录因子结合基序在启动子和类增强子元件中富集，其中启动子富含GCC/GGC-rich 基序，而类增强子元件富含GA/TC-rich和AT/TA-rich 基序。这个结果与在人类的早期研究中类似，即GA二核苷酸重复DNA序列可以广泛地用于预测活性增强子。启动子和类增强子区域之间序列偏好的差异可以反映这两种类型的调节元件之间与转录因子结合的多样性。
图5 增强子预测及实验验证
a) H3K4me1（b），H3K4me3（c），DHS（d）和H3K9ac在染色质状态5（S5）的基因近端和基因远端、染色质状态区域1 （S1）和2（S2）周围的覆盖密度分布。
e) 顺式元件中启动子和类增强子序列的基序富集分析。
7 类增强子序列的功能鉴定
我们基于荧光素酶报告系统验证了类增强子元件的调控能力。选取了26个具有不同DHS密度且位于基因远端的类增强子序列（距离最近基因的距离> 20kb）进行实验验证（图6a）。将这些序列插入到35S启动子控制下的报告载体中，之后将重组质粒在本氏烟草中瞬时转染（图6b展示部分结果）。在检测的六个序列中，具有较高DHS信号的三个区域在报告实验中表现出稳定而较强的活性（图6c），其信号比单独携带35S启动子的对照组高2倍以上。通过定量比较表观信号的强度和实验获得的增强活性，我们发现DHS，H3K9ac与增强子活性相关性最强（图6d-e）。综上所述，我们的结果表明染色质特征状态可用于鉴定具有增强子功能的基因组区域。
该研究工作由中国科学院植物生理生态研究所,南京农业大学和中国科学院遗传与发育生物学研究所团队合作完成。中国科学院植物生理生态研究所的张一婧研究员和南京农业大学的张文利教授为论文的共同通讯作者；中国科学院遗传与发育生物学研究所的薛勇彪研究员与童依平研究员参与项目的设计与指导，博士生李子娟、王梅月、林堪德、谢忆琳为共同第一作者。该研究受到中科院战略科技先导专项、基金委和教育部项目的资助。
文章链接：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1746-8