辐射诱导的染色质拓扑关联结构域层级变化及其在细胞辐射响应中的作用

陶 欢， 伯晓晨， 陈河兵*， 郑晓飞*

(1)军事科学院军事医学研究院辐射医学研究所，放射生物学北京市重点实验室， 北京 100850；2)军事科学院军事医学研究院卫生勤务与血液研究所， 北京 100850)

染色体构象捕获技术的发展与应用加速了基因组三维结构的探索。细胞核内每条染色体都占据着一个独特且不重复的区域，称为染色体疆域。染色体疆域中包括与常染色质和异染色质分别相关的A、B染色质区室，染色质区室由拓扑关联结构域(topologically associated domain，TAD)组成，并可进一步细分为染色质环，例如调控基因表达的增强子-启动子环等[1, 2]。环挤压模型认为TAD是由染色质环挤压形成的，环状挤压因子——黏连蛋白沿着染色质移动，当遇到另一个挤压因子或两个相反方向的边界元素——CTCF(CCCTC binding factor，CTCF)时，将其挤压成环，继而形成TAD[3, 4]。

TAD不仅是DNA复制、基因转录的基本功能单位，同样是DNA损伤修复的功能单元[5]。已有研究表明，DNA双链断裂(DNA double strand break，DSB)发生时，DSB位点的黏连蛋白不断积累，诱导断裂位点的单侧环挤出和扩大，并最终被阻滞在TAD边界，当含H2AX的核小体通过DSB位点锚定的黏连蛋白时，会迅速被DNA损伤修复关键激酶毛细血管扩张性共济失调突变(ataxia telangiectasia-mutated，ATM)蛋白磷酸化，使TAD内部以正确的方式建立γH2AX-53BP1(P53-binding protein 1)染色质域，该区域可作为DNA损伤修复位点的中心，使γH2AX迅速扩散至TAD边界，在TAD内部形成“易于DNA损伤修复”的环境，维持基因组三维结构的稳定性[5-7]。此外，招募到DSB位点的黏连蛋白将未受损基因隔绝在γH2AX募集区域外，保护基因免受染色质结构变化的影响，保证基因的正常转录，这种黏连蛋白依赖的DNA损伤响应机制从酵母到人类均保守，且与DSB修复途径无关，而募集到TAD边界的CTCF可以促进断裂的DNA双链进行同源重组修复[5, 6, 8, 9]。

TAD在不同的细胞系和物种中具有较高的保守性，其边界有大量的CTCF和黏连蛋白富集，能够形成相对独立的局部环境[6]。近期研究表明，TAD并非是完全独立的染色质结构，其内部常呈现多层级结构，与染色质的表观修饰及基因表达密切相关[10, 11]。根据TAD是否具有多层级结构，可将其分为两类，一类是有多层级的TAD结构，一类是无多层级的TAD结构，两类不同的TAD结构在生物学功能上具有差异[12]。比较而言，多层级的TAD具有更活跃的表观遗传状态和较高的基因表达水平[12, 13]。此外，多层级的TAD间常共享同一TAD边界，这些边界常表现出非对称性，较非共享边界具有更活跃的表观遗传修饰和基因表达，且富集大量的启动子，这些TAD边界在基因组的结构维持及基因表达调控中具有重要作用[12, 13]。

已有文献报道，TAD在DNA损伤修复中具有重要调控作用，TAD内染色质环挤出能够确保DNA损伤修复结构的建立，并避免基因表达受到染色质结构变化的影响[5, 6]。那么辐射诱导DNA损伤时，TAD多层级结构的拓扑变化是否在细胞辐射响应过程中发挥作用？本文基于前期已发表文献中5Gy X射线照射后的淋巴细胞、成纤维细胞以及ATM缺陷的成纤维细胞的Hi-C测序数据，探索了不同细胞系，辐射处理后不同时间，TAD多层级结构的变化规律，提出了辐射诱导TAD多层级结构变化调控基因表达参与辐射响应的模型，为理解染色质高级结构如何参与细胞辐射响应提供了新视角。

1 材料与方法

1.1 Hi-C数据集

本研究使用的26个Hi-C数据来自于GEO数据库(GSE136899)，包括淋巴细胞(GM12878)、成纤维细胞(BJ和MRC5，其中BJ细胞包含两个克隆：BJ1和BJ5ta)和ATM缺陷的成纤维细胞(ATM基因突变的BJ细胞，ATM抑制剂KU-55933处理的BJ细胞和ATM缺失病人的原代成纤维细胞GM02052)，5Gy X射线处理后30 min，24 h和5 d的Hi-C数据[5]。每个样本的Hi-C数据来源于5×106个细胞的测序结果，每个细胞包含约100～200个DSB位点[5]。本研究中使用的Hi-C数据为Sanders等在GEO数据库(GSE136899)中提供的经过质控和归一化处理后的Hi-C相互作用矩阵[5]。

1.2 拓扑关联结构域多层级结构的识别

TAD多层级结构通过OnTAD算法识别，该算法是一种优化的巢式TAD结构识别算法，可以从Hi-C数据识别多层级TAD结构[12]。OnTAD通过两步实现对TAD多层级结构的识别：首先通过自适应局部最小化搜索算法，用WxW的菱形滑动窗口扫描并识别候选TAD边界；随后被识别的候补TAD边界以基于评分的递归动态规划算法组装成多层级的TAD结构[12]。OnTAD方法将TAD分为5个层级，其中最外侧的TAD为层级1，层级2为层级1直接嵌套的子TAD，以此类推。以40 kb分辨率的Hi-C矩阵作为OnTAD算法的输入，将最大的TAD尺度定为50个bin(识别的最大的TAD为2 Mb)，将最小的TAD尺度定为1个bin(识别的最小的TAD为40 kb)，将“penalty”参数设置为 0.1。为保证样本间的可比性，对所有样本的计算使用同样的参数。同时，对TAD边界进行定义：根据某一边界两侧TAD的数量，以TAD数量多的一侧作为该边界的层级，TAD边界同样分为5级。

1.3 拓扑关联结构域边界层级的定义

通过OnTAD识别各样本不同层级TAD，得到TAD两端的位置信息，将TAD起点和终点的上下游60 kb定义为TAD边界。对26个样本进行TAD边界识别，将重复样本的TAD边界层级进行整合，取两重复样本中TAD边界层级数高的作为该样本的TAD边界层级。

1.4 拓扑关联结构域边界的聚类分析及可视化

在获得TAD边界信息后，将26个样本的TAD边界合并分析，采用ward层次聚类方法对26个Hi-C数据中保守的TAD边界进行聚类，并使用R包corrplot对不同样本TAD边界间的相关性分析，并进行可视化。

1.5 拓扑关联结构域边界层级变化的显著性检验

通过独立样本的两比例Z检验判断淋巴细胞、成纤维细胞以及ATM缺失的淋巴细胞辐射后30 min、24 h的不同层级TAD边界变化比例是否存在显著差异。N.S.:P>0.05;*P<=0.05;**P<=0.01;***P<=0.001;****P<=0.0001。

1.6 基因的拓扑关联结构域嵌套层级分数计算及统计学分析

TAD嵌套层级分数(TAD hierarchical score，TH score)的计算由bedtools软件的coverage功能实现，采用hg19基因注释文件进行注释，具有重复样本的数据，取基因在2个重复样本中TH分数(score)的平均值，并对ΔTH score(实验组-对照组)>1.5的基因进行筛选[14]。采用R包ggpubr中的非参数Wilcoxon符号秩和检验，并计算配对样本间的基因TAD层级分数差异。

1.7 基因功能富集分析

通过R包enrichR对GM12878细胞辐射后30 min TH score变化较大的基因进行功能富集分析[15]。

1.8 基因组位置注释

基因组位置信息通过UCSC基因组浏览器(http://genome.ucsc.edu)注释[16]。

2 结果

2.1 Hi-C数据的收集与分析

为了探究辐射后细胞TAD层级结构变化规律及其在细胞辐射响应中的作用，本研究选用了前期已发表的淋巴细胞，成纤维细胞和ATM缺陷的成纤维细胞，经5Gy X射线处理后30 min，24 h和5 d 的Hi-C数据(Fig.1A)。获得26个样本的Hi-C数据，对各样本中的TAD层级、TAD边界层级进行识别并计算了基因的TAD层级嵌套分数(Fig.1B)，从细胞系、ATM及辐射后的时间三个方面，对TAD边界层级和基因的TAD层级嵌套分数变化进行分析，比较不同辐射敏感性的细胞，在辐射后不同时间TAD层级结构的变化规律，探索ATM缺失对辐射后细胞TAD层级结构变化的影响及TAD层级结构变化在细胞辐射响应中的作用(Fig.1C)。

Fig.1 Scheme of analysis (A) 26 Hi-C data obtained from GEO database including irradiated lymphoblasts, fibroblasts and fibroblasts deficient in ATM with 5 Gy X-ray. (B) Identification of hierarchical TAD and TAD boundary and calculation of TH score. (C) Data comparison from three aspects: ① cell line，② ATM and ③ time after radiation

2.2 多层级拓扑关联结构域及边界的识别

使用OnTAD算法对26个样本中的多层级TAD进行识别，并定义TAD边界层级。结果显示，TAD多层级结构普遍存在于淋巴细胞和成纤维细胞中，以GM12878细胞系1号染色体的18.09 ～ 18.28 Mb区域为例，该区域为多层级的TAD嵌套结构，外部为大的1级TAD，内部嵌套着多层级的TAD结构，并包含有不同层级的TAD边界(Fig.2A、2B)。对各样本不同层级的TAD和TAD边界数量进行统计，平均每个样本识别出了 4 054(±206)个 TAD 结构和3 805(±149)个 TAD 边界，3层级及以上的TAD和TAD边界数量较少(Fig.2B)。对26个样本中的保守的TAD边界进行相关性分析并聚类，可以区分淋巴细胞，成纤维细胞和ATM缺陷的成纤维细胞(Fig.2C)，说明TAD边界层级具有一定细胞特异性。

Fig.2 Identification of the hierarchical TAD and TAD boudary (A) Multilevel TAD and TAD boundary of the subregion on chromosome 1 of GM12878 cell visualized by 40 kb-resolution Hi-C data. (B) Numbers of TAD and TAD boundary in 26 samples. (C) Hierarchical clustering of samples using the similarity of TAD boundary positions. Color bar shows the correlation coefficient between different samples

2.3 辐射诱导的拓扑关联结构域边界层级变化

为了确定辐射对TAD多层级结构的影响，对5Gy X射线照射后30 min和24 h的淋巴细胞(GM12878细胞)和成纤维细胞(BJ细胞)不同层级的TAD边界与相应对照组进行比较。与未照射细胞相比，辐射后30 min和24 h，GM12878细胞和BJ细胞的各层级TAD边界均发生了较大变化，其中发生变化的1层级和2层级TAD边界有30% ～ 40%，发生变化的3层级及以上TAD边界有40% ～ 50%。说明高层级TAD边界对辐射更敏感，辐射条件下细胞的低层级的TAD结构更加稳定。此外，GM12878细胞与BJ细胞辐射后30 min和24 h的TAD边界层级变化表现出了不同的模式，GM12878细胞辐射后24 h较辐射后30 min，随着TAD边界层级升高，未变化TAD边界和获得层级的TAD边界比例增加，缺失层级的TAD边界的比例减少，BJ细胞辐射后24 h较辐射后30 min，不同层级的未变化TAD边界比例均下降，而层级增加的TAD边界比例变化较小，缺失层级边界的TAD边界比例上升(Fig.3A)。GM12878细胞与BJ细胞辐射后变化模式的不同，可能与淋巴细胞和成纤维细胞的辐射敏感性差异相关。为了观察高层级TAD对辐射更敏感的规律是否在基因组范围内存在，对GM12878细胞5Gy X射线照射后24 h各染色体1层级TAD边界和3层级TAD边界变化情况进行分析，发现照射后GM12878细胞各染色体的低层级TAD边界均表现出更高的保守性，说明辐射诱导的TAD边界层级变化规律在全基因组范围内存在(Fig.3B)。

Fig.3 Changes of the hierarchical TAD boundary (A) Proportion of the hierarchical TAD boundary changed/unchanged in GM12878 and BJ cells after (30 min and 24 h) exposure to 5 Gy X-ray compared with control. N.S.: P > 0.05;*: P <= 0.05;**: P <= 0.01;***: P <= 0.001;****: P <= 0.0001). (B) Proportion of level 1 and level 3 TAD boundary changed/unchanged in GM12878 after (24 h) exposure to 5 Gy X-ray compared with control

2.4 辐射诱导的基因拓扑关联结构域嵌套水平变化

前期研究表明，TAD层级结构与基因表达密切相关[14]，为了探究辐射诱导的TAD层级结构变化是否参与细胞辐射响应，对各样本中18 000多个基因的TAD层级嵌套分数“TH score”进行计算，并对各样本中基因的TH score是否存在差异进行了统计学分析。结果显示，辐射后30 min和24 h，GM12878细胞和BJ细胞的TH score均发生显著变化，辐射敏感的淋巴细胞辐射后30 min基因的TH score变化较大，而辐射相对不敏感的成纤维细胞基因的TH score变化较小，说明辐射后TAD层级结构变化与细胞辐射敏感性相关。辐射后5 d，BJ细胞的TH score与对照组无显著差异，说明BJ细胞完成辐射损伤修复后，基因的TAD层级恢复至初始水平(Fig.4A)。为了验证细胞辐射后基因的TAD层级嵌套水平变化是否参与细胞的辐射响应，基于TH score变化较大的淋巴细胞辐射后30 min的Hi-C数据，筛选出TAD嵌套结构变化较大的基因(Fig.4B)，进行基因功能富集分析，发现TH score变化较大的基因富集了与淋巴细胞的功能紊乱及细胞辐射损伤修复相关的多种生物学功能，例如淋巴细胞黏附、B细胞激活、DNA双链断裂修复和细胞凋亡等，说明辐射诱导的基因TAD嵌套水平改变会影响细胞的生物学功能，并参与细胞的辐射响应(Fig.4C)。

2.5 毛细血管扩张性共济失调突变基因在辐射诱导的拓扑关联结构域层级结构变化中的作用

为了评估ATM在辐射诱导的TAD边界层级变化中的作用，对X射线照射后24 h的BJ细胞与ATM缺陷的BJ细胞各层级TAD边界变化进行比较，发现ATM突变后，BJ细胞辐射后24 h的TAD层级结构变化减弱，未变化TAD边界变化比例显著增加，说明ATM在辐射诱导TAD边界层级变化中具有重要作用(Fig.5A)。对ATM抑制和突变的BJ细胞5Gy X射线照射后24 h和对照组的TH score进行比较，发现正常BJ细胞辐射后5 d，基因的TH score恢复正常(Fig.4A)，但ATM缺失的基因TH score并未完全恢复(Fig.5B)，说明ATM对辐射后的基因TAD层级嵌套水平的恢复具有重要作用。

Fig.4 Radiation-induced changes of TH score (A) The TH score of lymphoblasts and fibroblasts. Statistical significance between two independent samples was calculated by the wilcoxon test (N.S.: P > 0.05;**: P <= 0.01;***: P <= 0.001;****: P <= 0.0001). (B) TH score heatmap of GM12878 before (control) and after (30 min) exposure to 5 Gy X-ray (gene names are partially displayed). (C) Gene ontology analysis of genes whose TH scores were significantly changed (ΔTH > 1.5, bule: TH score increased, orange: TH score decreased in GM12878 cells after (30 min) exposure to 5Gy X-ray

Fig.5 Changes of the TAD hierarchy in ATM-deficient cells (A) Proportion of hierarchical TAD boundary changed/unchanged in BJ cells and ATM mutant BJ cells after exposure (24 h) to 5 Gy X-ray compared with control. (B) TH scores of ATM-deficient BJ cells. Statistical significance between two independent samples was calculated by the wilcoxon test (N.S.: P > 0.05;**: P <= 0.01;***: P <= 0.001;****: P <= 0.0001)

2.6 拓扑关联结构域层级结构变化参与细胞辐射响应模型

前期研究表明，TAD层级结构与基因表观遗传状态及表达活性密切相关，高层级TAD结构具有更活跃的表观遗传状态，随着TAD层级增加，GM12878细胞的基因表达密度明显上升，且TAD内部嵌套的基因表达活性比外层更加活跃[12]。辐射后淋巴细胞和成纤维细胞的TAD边界层级及基因TAD嵌套水平发生明显改变，辐射后30 min淋巴细胞TH score变化较大的基因主要参与了细胞生物学功能变化及辐射响应，辐射诱导的淋巴细胞TAD层级结构变化可能通过改变基因的TAD嵌套水平调控基因表达参与细胞辐射响应。辐射诱导TAD层级结构变化，及其参与细胞辐射响应机制可总结为以下模型：首先细胞受到射线照射后，细胞的TAD层级结构发生显著改变，基因的TAD嵌套水平随之变化(Fig.6A)；导致基因局部的染色质相互作用改变，调控基因表达(Fig.6B)；基因表达变化一方面导致细胞的生物学功能紊乱，另一方面通过DNA损伤修复、细胞周期调控等生物学过程参与细胞的辐射响应(Fig.6C)。

Fig.6 Models of radiation-induced TAD hierarchy change regulating cellular responses to radiation (A) Radiation-induced TAD hierarchy change. (B) Regulation of gene expression. (C) Radiation-induced cell dysfunction and cellular responses to radiation

3 讨论

基因组三维结构是DNA复制、基因表达调控、变异和进化的基础，阐明基因组功能是理解遗传物质如何决定细胞命运的关键。TAD层级结构分析工具近年才开发，目前对于辐射诱导TAD层级变化的机制及其在细胞辐射响应中的作用尚未有深刻探讨。本文使用OnTAD算法分析了成纤维细胞、淋巴细胞和ATM缺陷的成纤维细胞5Gy X射线照射后30 min、24 h和5 d的TAD多层级结构，总结了不同细胞系、辐射后不同时间、不同层级TAD及边界的变化规律，并观察了ATM缺失对TAD层级变化的影响，发现TAD多层级结构变化是细胞辐射响应的潜在机制，并提出了辐射诱导的TAD层级结构变化参与细胞辐射响应的模型，为理解辐射条件下，细胞如何通过TAD层级结构变化维持基因组三维结构稳定性提供了基础，同时从TAD层级变化角度解释了ATM诱导细胞辐射响应的潜在机制。

前期研究表明，辐射后TAD边界通过富集CTCF和黏连蛋白得到加强，以维持基因组三维结构稳定[5, 6]。TAD可以调控DNA损伤修复关键分子在辐射诱导的DNA断裂位点的募集和扩散，一方面促进TAD内部DNA损伤位点的快速修复，另一方面限制γH2AX在DNA损伤位点局部扩散，以免影响其他未受损位点基因的表达[5, 6]。本文从一个新的角度——TAD多层级结构出发，分析了辐射条件下TAD多层级结构的变化规律，发现辐射对高层级TAD结构影响更大，低层级TAD相对保守。辐射诱导的TAD层级结构变化可能通过调节细胞辐射响应相关基因的表达，参与DNA损伤修复、细胞周期调控、凋亡等生物学过程，参与辐射后细胞的命运决定。

本研究从TAD边界层级和基因的TAD嵌套水平角度分析了辐射对TAD多层级结构的影响。5Gy X射线照射后5 d，成纤维细胞的TAD边界层级未完全恢复，而基因的TAD层级结构已经恢复至对照水平，这可能是染色质空间结构的高度动态变化导致的。辐射诱导的染色质的拓扑性质改变，伴随着TAD边界的不断变化，这种变化既包括TAD层级结构改变导致的边界层级丢失或获得，也包含TAD边界位置的动态变化。成纤维细胞基因的TAD层级在辐射5 d后恢复至对照水平，说明尽管辐射后TAD层级边界发生了复杂的动态变化，但在TAD多层级结构的总体调控上，仍能保证基因的TAD层级恢复正常，保证细胞基因表达的稳定性。此外，ATM缺失的成纤维细胞辐射后5 d，其基因嵌套水平与对照组仍存在明显差异，说明ATM作为辐射诱导的DNA损伤修复的关键激酶，不仅可以参与TAD内DNA损伤修复中心的形成，还能通过调控辐射诱导的TAD多层级结构变化和恢复，参与细胞辐射响应。辐射会导致TAD层级结构发生改变，伴随基因TAD的嵌套水平显著变化，5Gy X射线照射后30 min，GM12878细胞发生TAD嵌套水平显著变化的基因与细胞辐射损伤修复功能密切相关。我们的前期研究表明，TAD层级结构与基因表达水平密切相关[14]，据此我们提出了辐射诱导TAD层级结构变化，及其参与细胞辐射响应机制的模型：辐射诱导TAD层级结构变化介导基因表达改变，使细胞的生物学功能紊乱，并通过DNA损伤修复、细胞周期调控等生物学过程参与细胞的辐射响应。然而，Sanders等人提供的Hi-C数据缺少与之匹配的转录组测序数据，辐射诱导TAD层级结构变化与细胞辐射响应基因表达之间的相关性还需进一步验证。

TAD层级结构可能是细胞辐射敏感性的决定因素之一。尽管TAD结构在不同细胞和物种中高度保守，但不同细胞系的TAD层级结构仍存在一定异质性，例如本研究中的淋巴细胞和成纤维细胞，成纤维细胞中低层级TAD比例相较于淋巴细胞细胞高，而低层级的TAD结构可能使成纤维细胞的基因组三维结构更加稳定，辐射造成其TAD层级结构发生变化需要的能量更大，使其TAD多层级结构在辐射条件下更加保守。由于Hi-C数据的测序结果来自于数百万个细胞，所以分析得到的TAD多层级结构通常是许多细胞TAD结构叠加的结果[17]。细胞的TAD层级结构具有一定的异质性，即同一细胞系中存在不同TAD层级分布模式，基因组上的同一位置，部分细胞可能表现为高层级TAD结构，部分细胞表现为低层级TAD结构，这些细胞在辐射过程中受到的辐射损伤及响应机制可能存在差异。随着单细胞测序技术逐渐成熟，目前已经可以获取单个细胞的染色质高级结构，基因表达和表观修饰的数据，单细胞数据分析算法的开发也为细胞辐射响应异质性机制研究提供了条件[18-21]。随着单细胞测序技术的成熟和相关算法的进一步开发，染色质高级结构的异质性及其对细胞辐射敏感性的影响将得到更深刻的理解。