表观遗传机器对真核生物的基因调控起到重要作用。但是，目前对各个物种中> 600个表观遗传学调控子（epigenetic regulators，ERs）包括负责组蛋白修饰和DNA修饰的阅读、改写和擦除的基因的系统认识尚不全面。本研究基于分子进化分析构建了一个包含25个物种的ERs的列表资源，这是当前最为全面的ER列表。进一步的人、小鼠以及其它羊膜类动物的比较转录组学分析表明多物种间一致的组织来源的器官具有相近的基因表达特异模式，表明了ER基因进化和表达特异性以及基因功能之间的关系。此外，物种进化过程中的ER基因复制事件可以产生针对相同表观遗传学底物的组织特异性ER基因，这有助于它们在组织发育中发挥调控特异性功能。该工作有助于从分子进化角度更好地理解ER基因的组织特异性表达模式以及组织特异性表观遗传学调控的功能意义。

表观遗传学调控子（epigenetic regulators，ERs）中的主要组分是染色质调控子（chromatin regulators，CRs）,包含可催化、调控、擦除DNA或组蛋白上的化学修饰的一类蛋白。例如，组蛋白乙酰转移酶和去乙酰转移酶可以从组蛋白N末端尾巴上增加或去除乙酰化基团。除了染色质调控子外，RNA修饰调控子也是最近研究的热点。具有表观遗传学修饰结构域(例如含溴结构域)的基因可以从功能上分为以下四个类别：reader（识别者）, writer（写入者）, eraser（擦除者）和remodeler（重塑者）。ER基因涉及下游基因表达的正确调控，对于细胞发育和组织功能至关重要。ER基因的扰动是引起癌症和一些其它疾病包括神经系统疾病的关键因素。之前有几个ER的数据库，包括EpiFactors，dbEM以及CR2Cancer，都是汇总人类癌症的ER基因及其功能基因组注释。此外，WERAM数据库包含了8个物种的组蛋白乙酰化和甲基化调控子的信息。然而，对ER基因的功能研究大部分集中于人类，在其它物种中并没有全面的基因或蛋白列表。因此，构建一个囊括真核生物的全面的ER列表的数据库资源十分必要。

       于是，国科大温州研究院癌症诊断分子靶标识别研究团队（网站：https://www.x-mol.com/groups/bioinfolv）PI吕杰研究员等人构建了一个关于ER的全面的蛋白质组资源。该数据资源包括23个物种以及2个细菌物种的ER的蛋白信息，这些物种基本包含了真核生物进化树中的代表性物种。进一步构建了系统发生树，并识别了进化中复制基因（duplicated gene）和复制事件（duplication events）。进一步的分析表明复制基因和祖先基因（被复制基因）具有显著相似的基因表达和功能。对人和小鼠的ER基因进行比较转录组学分析结果显示ER基因在两个物种中的表达模式相近；进一步分析羊膜类动物的ER基因表达也发现同样的模式。此外，研究还发现ER基因的复制事件可以产生具有相同表观遗传学底物的ER基因，但在不同组织中具有特异性表达模式，这可能是不同组织器官的特异性调控功能所必需的。

数据集：(1) GTEx计划的人类表达数据; (2) 小鼠ENCODE计划的表达数据; (3) 10个羊膜类动物的基因表达数据（GSE30352）。基因表达数据通过TPM (transcripts per kilobase million)进行量化。（4）人类ERs基因列表来自FACER论文，共690个人类ERs。

复制事件识别方法：选择23个具有代表性的真核生物以及2个细菌物种(bacteria和archaea)作为外类群进行系统发育分析。首先，从UniProt数据库下载25个物种的蛋白质FASTA文件。其次，从pFAM数据库中取得45个具有代表性的表观遗传学结构域的谱文件。再次，使用HMMER软件中的hmmsearch命令对从蛋白FASTA文件中搜索每个Pfam谱中可能的人类同源基因（默认阈值）。最后，汇总ER的同源蛋白表格。为了防止计算扫描过程遗漏同源蛋白，还从UniProt数据库中获取同源蛋白信息进行补充。使用OrthoFinder软件预测基因复制事件，进化树是通过FigTree软件进行可视化。

人类和小鼠的组织特异性表达分析：组织特异性表达是基于ER基因的表达数据的TPM数据进行计算。为应对TPM数值为0的情况，对TPM加上伪计数并计算log表达值。组织特异性分数（Tissue specificity score，TS score）是通过R的rlm函数计算的残差进行计算。一个基因在一个组织中的TS score大于所有组织的2.5倍的残差的标准误即可被认为是组织特异的基因。该分数可以识别到被传统的2-倍差异法舍弃的基因，同时具有较低的假阳性率。

羊膜类动物的组织特异性表达分析：针对羊膜类动物进行ER基因的比较转录组学分析，包括10个物种：人类、黑猩猩、倭黑猩猩、大猩猩、婆罗洲猩猩、恒河猴、小鼠、短尾负鼠、鸭嘴兽和红原鸡。基于RPKM度量计算表达水平，并进行log转换（加上伪计数）和Z尺度变换，最终得到Z score。

基因集富集分析：使用ReactomePA R包进行ReactomePA通路富集分析（图5e,f）。Enrichr网站进行不同羊膜类动物内的组织差异基因的通路富集分析（图5d）。评估duplicated基因的功能富集性（图1f）的基因集和数据集包括（1）本文作者之前开发的预测癌症驱动基因的DORGE软件得到的1172个基因。（2）https://www.tau.ac.il/~elieis/HKG/下载的3804个看家基因Housekeeping gene (HKG)列表。（3）OGEE数据库下载的必需基因和duplicated genes列表。由于必需基因过多，为了使结果更可靠，只保留了2340个在原始列表中出现次数大于两次的基因。使用GOSemSim R包计算基因本体论相似性。单尾Fisher精确检验用来评估基因集在duplicated genes列表中的富集性（图1f）。

统计学分析：PCA分析基于R 4.2.0的prcomp函数进行分析，图片是由前两个主成分PC1和PC2绘制（图3c-g，图5a,b）。PCA图中的离群点是由最大的10%的PC1和PC2的马氏距离定义（R mt包）。皮尔森相关系统是由R的cor函数进行计算（图2b）。使用双尾Wilcoxon秩和检验评估箱线图的组间差异（图1c-e）。

1）基因复制事件的功能基因组学评估

本研究通过分子进化分子方法识别新的表观修饰酶蛋白结构域。首先，收集表观修饰相关的酶蛋白结构域和癌症组学数据。由于表观遗传学结构域数量较多，为了减少分析时间，版面简洁，只分析了45个之前研究表明比较重要的表观遗传学结构域。同源基因预测是通过基于hidden Markov model (HMM)模型的蛋白质序列搜索得到的。基因复制（gene duplication）是一个特定物种获得适应环境的能力所必需的一种机制。使用OrthoFinder软件识别了25个物种的45个结构域的2482个基因复制时间。汇总的基因复制数目显示：ER从脊椎动物开始集中出现（蓝色深浅和酶蛋白的数量正相关），暗示大量之前未报道的ERs可能是脊椎动物以及更高等的动物进化所必须的ERs（图1a）。BRD1（BRPF2）基因的进化树展示了该基因进化的2个复制点，以这两个复制点基因为起点（祖先），新生成的2个新的基因BRPF1和BRPF3（图1b）。之前的文献报道显示这3个基因具有较大差异的生化性质，即结合亲和性，表明基因复制事件产生的新基因和祖先基因的功能可能有所不同，甚至新基因之间也具有不同的功能。

进一步系统评估人类的duplicated genes和祖先基因之间的表达水平和功能关系。与预期一致的是，复制事件关联的基因之间确实存在统计学显著的表达相似性（图1c, p-value = 2.53e-131)。这里的对照数据是来自打乱的复制事件关联的关系。此外，同样发现基因功能相似性，包括基因本体论（Gene Ontology，GO）生物学过程（BP）和分子功能（MF）分支（图1d, p-value = 2.97e-124；图1e, p-value = 3.26e-255）。

进一步的功能基因组学基因集的富集性表明duplicated genes富集了癌症相关基因（Fisher’s exact test, p-value = 1.86e-98）、必需基因（Fisher’s exact test, p-value = 1.74e-60）和看家基因（Fisher’s exact test, p-value = 3.19e-74）（图1f）。癌基因的富集性和之前报道的体细胞突变率可以由进化研究中经常使用的测度进行建模以及癌基因早在真核生物的早期进化过程中就已经出现（Lyu J, et al., 2020）的观点相吻合。此外，OGEE数据库中的duplicated genes同样富集这里识别的duplicated genes（Fisher’s exact test，p-value = 3.85e-12），但是本研究识别了更多的候选duplicated ER基因。

图1. ER基因的分子进化分析和功能基因组学表征。

2）人与小鼠的ER基因的比较转录组学分析

进一步，拟考察不同物种间被进化出的新的ER基因是否具有与祖先基因不同的表达模式。基于公共数据库中的基因表达数据集研究了人和小鼠的不同组织的ER基因的基因表达模式（图2A）。对于每个基因和组织，计算特异性分数TS score，通过在一个给定组织中的基因表达相对于所有组织的平均表达水平而计算出来。有趣的是，利用该分数，人类和小鼠的ER基因表现出明显一致的相同组织来源的表达模式(图2a,b)，即使是不同的类型的ER基因（图2a）。

图2. 人和小鼠的ER基因的比较转录组学分析。

结果表明睾丸组织是所有ER基因中的最与众不同的组织，在人和小鼠中都是如此（图3 a,b左）。相反，像心脏这样的器官，几乎没有表现出组织特异性表达模式（图3 a,b右）。这里，还使用主成分分析（PCA）方法分析人（图3c）和小鼠（图3e）组织的ER基因表达数据，结果发现了组蛋白写入者（writer）和重塑者（remodeler）在PCA图谱中具有更多的离群（outlier）基因（图3d）。PCA loading结果指出ER基因在生殖系统（睾丸、卵巢和胎盘）和神经系统（大脑、小脑和垂体）具有更大的基因表达差异，这和图3c和3e一致。

图3. 人和小鼠的ER基因的组织特异性表达分析。

3）具有相同组蛋白修饰底物的ERs具有不同的表达模式

基于之前研究定义的ER的组蛋白修饰底物的信息（包括H3K4/27/36甲基化修饰写入者/擦除者, H3K9 甲基化修饰写入者/擦除者, H3K27乙酰化写入者, 以及H3K9乙酰化擦除者），在人类组织表达数据中进行PCA分析。结果显示，H3K4甲基化写入者例如SMYD1, PRDM9以及SETD7在人类不同组织中具有明显不同的表达模式（图4a,b）。类似地，H3K27甲基化擦除者例如UTY和PHF8同样表现出差异较大的表达模式，并且这种表达差异不能由基因表达检测引起的噪声所解释（图4b）。总之，结果暗示了有相同组蛋白修饰底物的ERs在不同组织具有不同的表达模式，这种组织特异性表达模式和组织特异性功能可能是密切相关的。

图4. 具有相同组蛋白修饰底物的人类ER基因的表达分析。

4）羊膜类动物的ER基因的比较转录组学分析

接下来，试图将分析扩展到羊膜类动物中。这里分析的羊膜类动物包括10个物种。主成分分析表明少数ER基因在羊膜类动物的组织中展示出较高的基因表达特异性模式（图5a）。这里，使用主成分分析（PCA）方法分析ER基因表达数据，结果发现了DNA甲基化擦除者（erasers）在PCA图谱中具有更多的离群（outlier）基因（图5a）。由于这里使用的基因表达数据包含的基因没有人和小鼠的多，因此这里的统计学分析结果在很大程度上是被低估的。和图3类似，睾丸组织同样表现出较高的组织特异性（图5b）。比较有趣的是，其它组织在主成分分析图中也具有很高的区分度，表明组织间差异要远大于物种间的差异（图5b）。特别值得注意的是，大脑和小脑组织在PCA loading图中距离较近，这和这两个组织中大部分细胞类型是神经细胞的事实是一致的。

基于这10个羊膜类动物的ER基因的表达数据，TS score识别出72个组织特异性表达ER基因（图5c）。图5c同样说明了这些基因具有很明显的组织间差异性，特别是脑组织和睾丸组织和其它组织的表达差异很大。有趣的是，许多组织特异的ER基因实际上和癌症有关，例如TOP2A在众多癌症中频繁失调。一个可能的解释是TOP2A可以切断双链DNA并诱导有丝分裂过程中的基因转录。对这些差异表达基因进行功能富集分析发现它们在染色质组装、染色质重构和细胞周期调控等方面的GO条目富集，这和它们是表观遗传学调控基因的本质是一致的（图5d）。由于ER基因在神经系统的潜在重要作用，进一步分析脑和睾丸特异ER基因在Reactome通路中的富集性，结果显示脑特异ER基因富集于神经元相关的通路中（图5e），包括trafficking of AMPA receptors和RHO-GTPase effectors等；而睾丸特异ER基因富集于神经元相关的通路中（图5f），包括表观遗传调控、DNA修复、P53活性调控及染色质组装等通路，这可能归因于精子生成过程中所需要的大量的组织特异基因表达、DNA富集及修复（图5f）。总之，结果表明ER基因表达在不同物种的相同类型组织的中是相近的（图5b,c）。此外，组织特异基因倾向于在多物种中具有保守的功能，特别是和细胞周期调控和染色质组装相关的ER基因（图5d）。

图5. 羊膜类动物的ER基因的比较转录组学分析。

  过去30年间大量ER基因被广泛研究，但是尚缺乏一个全面汇总多物种ER基因的数据资源，这对于后续其功能基因组学和比较基因组学研究至关重要。在本文中，该研究在25个物种汇总了45个表观遗传学结构域相关的ER蛋白组信息，这和已有的ER基因资源是高度互补的。此外，该研究还从人、小鼠到10个羊膜类动物等维度分析了基因复制事件中ER基因的表达模式。

  总体上，ER基因在羊膜类动物中特别是人和小鼠之间，具有在相同组织中相近的表达模式；然而在同一物种中的不同组织中，一些特定的ER基因表现出差异较大的表达模式，这反映出ER基因的组织特异的调控模式，以及基因复制事件产生了一定程度的功能冗余的ER基因。从进化角度解释，真核生物为了适应复杂的生存环境，需要一定程度的基因冗余，以缓解负向选择的压力，而删除这种基因通常不会引起个体的致死。

  该研究还指出ER基因复制可以产生具有相同表观遗传学修饰底物的组织特异基因（图4）。例如，SMYD1, SETD7, PRDM9和KMT2C是四个H3K4甲基化写入者，它们在人类组织中表现出较大差异的表达模式（图4a,b）。其中SMYD1在心脏发育过程中具有表达特异性；而PRDM9在睾丸中特异表达，并涉及有丝分裂的调控；SETD7是一个涉及DNA损伤相应和染色质调控的基因，它可以将甲基化基团运送到底物的赖氨酸残基。这些例子表明具有相同表观遗传学底物的ER基因可以在不同组织中发挥特定的调控作用。

  总之，该研究基于蛋白质组序列搜索，汇编了多物种的ER蛋白列表。分析结果显示复制的ER基因与被复制的ER基因有相似的表达和功能。此外，比较转录组学分析还暗示了ER基因在转录组水平上的保守性以及组织间的差异性。进一步地，该研究还分析了具有相同表观遗传学修饰底物的ER基因的表达水平，暗示了基因进化和功能上的潜在联系，这为进一步理解表观遗传学修饰调控的细胞类型特异的表达及后续分子进化分析和功能研究奠定基础。本研究已在线发表在《Epigenetics》期刊（中科院升级版分区，生物学2区），该研究受到国家自然科学基金和国科大温州研究院启动经费的支持。