衰老細胞染色質開放變化的初步分析

宋 喬,王亞琦,侯玉麗,劉 靜,張曉敏,曹 敏,王培昌

(首都醫科大學宣武醫院檢驗科,中國北京100053)

復制性衰老又稱細胞衰老,當細胞周期停滯、凋亡抵抗以及部分基因如P16INK4a表達改變時,細胞出現衰老表型[1～2]。復制性衰老是整體衰老的基礎,近期研究發現其具有拮抗癌癥、使癌細胞得以清除的能力[3～4]。近年來,隨著人口老齡化的發展及癌癥等老年病發病率的增加,細胞衰老作為老年疾病的主要病理基礎,已成為分子生物學研究的熱點之一[5～8]。

在真核生物中,DNA圍繞組蛋白分級折疊,形成核小體,而核小體則通過進一步緊密折疊形成高級結構——染色質[9]。此時,根據細胞所處的狀態,其基因組可劃分為活躍表達的基因與相對沉默的基因。在DNA分級折疊的這一過程中,非活躍基因表達相關的結構區域相對封閉、核小體結構致密,呈轉錄抑制狀態；而與活躍表達基因相關的區域則呈開放狀態,核小體相對松散,有利于轉錄因子、組蛋白等調節元件與這些基因的啟動子、增強子相互作用,促進轉錄的進行[9～10]。通常情況下,染色質的開放程度可代表轉錄因子、組蛋白等調節元件與DNA的結合狀態,因此根據染色質的開放程度明確生理學狀態、疾病發生發展中表觀遺傳學狀態的改變,也已成為表觀遺傳學技術的主要切入點[9,11～13]。

迄今為止,細胞衰老過程中表觀遺傳學狀態的改變尚不明確,從染色質空間結構、基因組學角度對染色質開放性的變化進行分析探討的研究較少,在細胞衰老過程中染色質開放程度的具體變化仍未知[14～15]。故本研究采用基于轉座酶Tn5和高通量測序的染色質開放性分析技術ATAC-seq(assay for transposase-accessible chromatin with highthroughput sequencing),研究在復制性衰老過程中全基因組染色質開放程度的變化,為探索復制性衰老過程中表觀遺傳以及轉錄調控水平的變化提供初步的研究基礎。

1 材料和方法

1.1 材料

年輕代齡(PD26)和年老代齡(PD55)的人胚肺二倍體成纖維細胞(2BS)購自武漢大學中國典型培養物保藏中心。

1640培養基、胎牛血清(fetal bovine serum,FBS)、青霉素鏈霉素混合液和磷酸鹽緩沖液(phosphate buffer saline,PBS)購自美國Gibco公司；Tris、NaCl以及MgCl2購自BBI生命科學有限公司；Tn5轉座酶源自美國Illumina公司；NP-40和nuclease-free H2O購自上海碧云天生物技術有限公司；DNA純化試劑盒源自美國Qiagen公司；PCR擴增試劑盒源自美國New England BioLabs公司。

二代測序儀(150PE,Illumina公司,美國)；qPCR儀(480,Roche公司,瑞士)；低溫高速離心機(CR3I,Thermo公司,美國)。

1.2 Tn5酶切反應

參照文獻[16]的方法,對2BS細胞進行Tn5酶切反應,主要步驟如下：取凍存于液氮中的PD26、PD55的2BS細胞進行復蘇。復蘇后的細胞置于T75培養瓶,采用添加10%胎牛血清、1%青霉素鏈霉素混合液的1640培養基培養,培養箱條件為37℃、5%CO2。當細胞長滿且生長狀態良好時,用胰酶將細胞消化,并使用臺盼藍染色,過濾除去死細胞后進行細胞計數,計數50 000個細胞。取計數后的細胞,在4℃條件下500 g離心5 min后棄上清；再加入50 μL預冷PBS吹打洗滌,4℃條件下500 g離心5 min,棄上清；加入50 μL預冷的裂解液(由10 mmol/L Tris pH 7.5、10 mmol/L NaCl、3 mmol/L MgCl2以及 0.1%NP-40 配置而成)懸浮細胞,在4℃條件下500 g離心10 min后棄上清。此時向細胞沉淀中加入50 μL轉座反應體系(表1),吹打混勻后,37℃孵育30 min。待孵育結束后使用DNA純化試劑盒純化DNA,純化完成后使用10 μL試劑盒中的洗脫液將其洗脫,隨后進行PCR擴增。PCR反應體系如表2所示,循環條件見表3。其中,PD26對應的primer 1序列為AGGCAGAA,PD55對應的primer 1序列為TCCTGAGC；barcoded primer 2序列為AATGATACGGCGACCACCGAGATCTACACTCGTCGGCAGCGTCAGATGTG。PCR結束后,再次使用DNA純化試劑盒純化DNA[17]。

表1 轉座反應體系Table 1 Transposition reaction mix

表2 PCR反應體系Table 2 PCR reaction mix

1.3 ATAC-seq測序、質量評估及原始數據過濾

使用Illumina Hiseq 150PE測序儀對Tn5酶切反應得到的DNA進行雙端測序,識別chasity＞0.6的堿基,對每個測序序列(sequenced read)進行pass filter檢驗(當某一read的前25個堿基中,僅存在一個或者不存在chasity≤0.6的堿基時,判定該read符合檢驗標準),過檢后的單克隆reads則使用FastQC進行測序質量評估。判定評估合格的標準如下：每個堿基的質量均≥5、每組堿基的中位數均≥20；tile的質量檢測中不存在暖色結果；序列質量評分均集中于高分；GC堿基含量分布均勻；GC堿基占比均勻；測序中序列長度分布較一致；測序序列中接頭(adapter)出現頻率較低；測序序列中特征性短序列重復出現的次數較低。如果評估合格,則進行數據過濾,去除低質量(QPhred≤20的堿基數占整個read長度50%以上的reads)、N(N表示無法確定堿基信息)的比例＞10%和包含adapter的reads,以保證分析質量。

表3 PCR循環條件Table 3 PCR thermal cycle

1.4 ATAC-seq的序列比對和結合位點(peak)掃描

使用Bowtie2(版本2.25)對參考基因組(hg19)建立索引回帖后,進行序列比對。然后根據泊松分布,使用MACS(版本2.1.0)對測序后各個區域基于唯一比對的read數的P-value值進行計算,當P-value＜1.0E-05時,該區域則被認為是一個結合位點(peak)。標準化后的P-value值越高,理論上該區域的染色質越致密[18]。

1.5 保守性分析

使用phastCons軟件,根據phylo-HMM模型對結合位點的保守性進行計算[19]。

1.6 回帖序列在基因附近的信號分布情況分析

將所有基因作為目標區域,使用deepTools對基因長度(轉錄起始位點TSS到轉錄終止位點TTS之間的距離)進行歸一化處理,繪制ATAC-seq reads(使用bigwig文件)在轉錄起始位點上游3 kb到轉錄終止位點下游3 kb范圍內的平均信號值,以average plot的形式將結果展示。

1.7 結合位點注釋分析

使用HOMER軟件的annotatePeaks.pl工具尋找離每個結合位點距離最近的轉錄起始位點所屬的基因(以下稱為結合位點鄰近基因),以及結合位點覆蓋的基因,并注釋其覆蓋的功能性區域。注釋的功能性區域包括啟動子到轉錄起始位點的區域(promoter-TSS,默認范圍為基因起點上游1 kb到基因起點下游100 bp)、轉錄終止位點(TTS,默認范圍為基因終點上游100 bp到基因終點下游1 kb)、外顯子(exon)、內含子(intron)和基因間區(intergenic)。如果一些功能性區域有重疊部分,則按上述順序優先選擇排在前面的功能性區域作為最終注釋。對上述注釋分析結果中結合位點在所覆蓋基因的功能性區域上的分布進行統計,并使用R語言的ggplot2包將統計結果繪制成餅狀圖和柱狀圖。

1.8 差異結合位點分析

使用bedtools工具將2BS-PD55與2BS-PD26的peak文件合并,并對合并文件進行處理。如果兩個結合位點有重疊區域,則合并成一個新的結合位點。計算PD26與PD55在每個合并的新結合位點區域上的read數目,然后使用DEGseq進行差異分析,得到樣本間的差異結合位點,即差異染色質開放區域。本次差異結合位點的篩選原則為|log2FC|＞1 且 P-value＜0.05。

1.9 差異結合位點注釋分析

使用HOMER軟件的annotatePeaks.pl工具尋找離PD26與PD55差異結合位點距離最近的轉錄起始位點所屬的基因,以及差異結合位點覆蓋的基因,并注釋其覆蓋的功能性區域。對差異結合位點注釋分析結果中結合位點在所覆蓋基因的功能性區域上的分布進行統計,并使用R語言ggplot2包將統計結果繪制成餅狀圖。

1.10 GO富集分析

基因本體論(Gene Ontology,GO)數據庫可對基因功能進行富集分析。我們使用GO Term Finder,將PD26與PD55差異結合位點臨近的基因與GO term數據庫比對,計算每個term的基因數量,然后應用超幾何檢驗,找出整個基因組背景下結合位點相關基因中顯著性富集的GO term,從而識別出差異結合位點相關基因行使的主要生物學功能。

2 結果

2.1 測序質量

2.1.1 測序的錯誤分布率

對Tn5酶切后的年輕與年老細胞中的DNA進行雙端測序,排除酶切后堿基質量因素以及測序儀和測序試劑等的影響后,所得reads的測序錯誤率如圖1所示,均低于0.06%。本實驗測序的錯誤分布率檢查合格、測序質量較高,提示此次測序結果可用于后續分析。

圖1 測序錯誤率分布圖(代表性結果)Fig.1 Error rate distribution along reads(a representative result)

2.1.2 堿基含量分布

使用FastQC對Tn5酶切后的年輕與年老細胞中的DNA進行堿基組成檢測,結果如圖2所示。各嘌呤堿基和嘧啶堿基的組成比例都較均衡,低質量堿基(指在某一read的某一位置中,含量小于20%的堿基)比例較低且僅出現于隨機引物擴增的開始階段,表明本次測序結果較好,可信度較高。

2.2 酶切片段長度的分布

酶切片段(insert size,也稱作fragment length)是指經轉座酶Tn5切割后獲得的DNA片段,常用于監測Tn5酶切反應中酶的用量是否合適。本實驗的酶切片段長度分布如圖3所示,從左至右依次為無核小體片段(nucleosome free,即開放染色質區域)、單核小體片段(包括一個核小體的酶切片段,長度一般在147 bp到147*2 bp之間)和雙核小體片段(包括兩個核小體的酶切片段)的長度分布,表明本次實驗酶切效果較好。

2.3 結合位點統計

對測序產生的reads的P-value值進行計算,通常標準化后的P-value值越小,理論上該區域的染色質越疏松[18],越可能成為染色質的開放區域從而結合轉錄因子。對本次測序樣本的peaks進行統計分析后發現,PD26的2BS細胞有45 681個peaks,而PD55的2BS細胞有50 075個peaks。

2.4 結合位點保守性分析

圖4展示了年輕細胞與年老細胞進化保守性的分析結果。在脊椎動物基因組學中,調節復制、轉錄和翻譯的序列,常被認為是功能性保守序列[20]。在本次實驗中,我們對結合位點進行統計分析發現,年輕細胞與年老細胞的染色質開放性變化區域較保守、相似性高,因此可以初步推斷這些結合位點的相關區域在轉錄表達中具有較高的重要性。

圖2 測序序列堿基分布示意圖Fig.2 Sequence content across all bases

圖3 酶切片段長度的分布X軸代表Tn5酶切后DNA片段的長度,Y軸代表該長度片段出現的比例。Fig.3 Fragment length distributionThe X-axis represents fragment length,and the Y-axis represents the frequency of length distribution.

圖4 結合位點保守性分析X軸代表此次檢測所得的結合位點對應基因組中該位點的相對位置信息,其中“0”表示相對位置的中點；Y軸代表保守性的評分,分值越高代表保守性越強。Fig.4 The results of conservationThe X-axis represents position information,and“0”represents relative position of peaks.The Y-axis represents conservation score(a higher score indicates a higher conservation).

2.5 回帖序列(mapped reads)在基因附近的平均信號分析

在所有基因長度歸一化后,將reads與hg19基因組比對,計算其在基因組上的分布情況,結果如圖5所示。從結果可知,回帖序列(即染色質的開放區域)主要集中在轉錄起始位點附近,提示這些染色質的開放區域與活躍的轉錄調控相關。同時,我們發現,隨著細胞衰老程度的加深,轉錄起始位點附近的染色質開放程度明顯下降,提示隨著細胞衰老的發生發展,轉錄水平進一步下降。

圖5 回帖序列在基因附近的平均信號分布X軸表示歸一化后的基因范圍,“-3.0”表示TSS上游3 kb,“3.0”表示TTS下游3 kb；Y軸表示位點的平均信號,數值越大代表越富集。Fig.5 The relationship between mapped reads and genes in average plotThe X-axis represents the normalized gene range,“-3.0”represents 3 kb upstream the TSS,and“3.0”represents 3 kb downstream the TTS.The Y-axis represents the average signal value of the site.The enrichment degree is proportional to the Y-axis value.

2.6 結合位點在全基因組功能區上的注釋

全基因組的功能區域分為啟動子到轉錄起始位點的區域、轉錄終止位點、外顯子區、內含子區和基因間區。轉錄因子一般集中分布在轉錄起始位置附近,即TSS區域。結合位點在全基因組功能性區域上的分布如圖6所示,年輕細胞(2BSPD26)、年老細胞(2BS-PD26)的染色質開放區段主要位于：內含子區、基因間區、啟動子到轉錄起始位點的區域、外顯子區和終止子。除去在全基因組中占比極大且不參與基因表達的內含子和基因間區后,染色質開放區主要富集于啟動子到轉錄起始位點的區域。此外,圖7結果顯示年輕細胞與年老細胞中結合位點在該區域分布的比例分別為23.04%和15.4%。以上信息表明轉錄調控對復制性衰老有著重要影響,且隨著細胞的逐漸衰老,啟動子到轉錄起始位點區域的開放比例進一步下降,即轉錄的水平進一步下降。

2.7 差異結合位點統計

通過分析PD55與PD26之間的差異結合位點,可以得到衰老狀態下特異性染色質開放區域。本次實驗樣本間的差異結合位點數目為15 812。

圖6 結合位點在基因組功能性區域的分布統計圖Fig.6 The annotation of peaks in the genome

圖7 結合位點在啟動子到轉錄起始位點區域的分布統計圖Fig.7 The annotation of peaks in promoter-TSS

2.8 差異結合位點在全基因組功能區上的注釋

2BS-PD26和2BS-PD55的差異結合位點在全基因組功能區的注釋如圖8所示。結果表明,細胞衰老后染色質關閉與開放的區域主要集中在內含子區和啟動子到轉錄起始位點區域。其中,染色質關閉的區域有40.31%分布于內含子區,有26.45%分布于啟動子到轉錄起始位點區域；染色質開放的區域有52.55%分布于內含子區,有6.18%分布于啟動子到轉錄起始位點區域。

同時,在注釋中我們可以發現基因間區占比較多,推測可能因為基因間區在基因組上的占比較大,故檢測到的結合位點相對其他區域更多,但此類結合位點并非真正的轉錄調控因子的結合位點。

圖8 差異結合位點在基因組功能性區域的分布統計圖(A)細胞衰老后,染色質關閉區域在基因組功能性區域的分布；(B)細胞衰老后,染色質開放區域在基因組功能性區域的分布。Fig.8 The annotation of different peaks in the genome(A)The distribution of closing chromatin regions’annotation；(B)The distribution of open chromatin regions’annotation.

2.9 差異位點靶基因的GO分析

對2BS-PD26和2BS-PD55的每個差異結合位點的鄰近基因進行GO注釋,部分結果如表4所示。

表4 差異位點靶基因GO富集分析Table 4 GO term enrichment analysis to the target gene of different peaks

圖9 差異位點靶基因的GO富集分析Fig.9 GO term enrichment analysis to the target gene of different peaks

選取差異位點鄰近基因顯著富集的GO條目繪制柱狀圖,結果如圖9所示。富集程度較高的有細胞進程(cellular process)、代謝(metabolic process)、生物性調節(biological regulation)、結合功能(binding)、催化活性(catalytic activity)等,表明復制性衰老與細胞自身調節緊密相關,而這些生物學過程,均與轉錄調控有密切聯系。

3 討論

盡管目前針對復制性衰老過程中表觀遺傳學變化的研究方法有很多,如使用DNA酶1酶切基因組后進行高通量測序的DNase-seq、使用微球菌核酸酶酶切基因組后進行高通量測序的MNase-seq和利用甲醛對染色質固定后進行高通量測序的FAIRE-seq,但此類方法對實驗組細胞數量要求極高、所需測序深度較高、實驗步驟繁瑣且可重復性極低[10,16]。ATAC-seq的開發與應用,極大程度減少了樣本的使用量、增加了實驗的可重復性[21]。本實驗中的測序質量分析顯示：堿基錯配率極低、堿基含量分布均勻,這再次證明該實驗方法的準確性。同時,ATAC-seq作為檢測表觀遺傳學變化的一種新方式,可對靶基因的表觀狀態給出更為準確的推測,具有一定的可挖掘性[22～23]。

本研究通過ATAC-seq技術結合生物信息學分析方法,對年輕代齡(PD26)和年老代齡(PD55)的2BS細胞展開了分析,初步發現,2BS細胞中結合位點的基因保守性均較強；ATAC-seq的回帖序列主要富集于轉錄起始位點附近；結合位點特異性地注釋于啟動子到轉錄起始位點區域,提示復制性衰老與轉錄調控具有較強的相關性。

文中在分析回帖序列在基因附近的平均信號分布時,發現染色質的開放區域主要富集在轉錄起始位點附近,且年輕細胞與衰老細胞在轉錄起始位點的開放程度不同,衰老時染色質的開放程度明顯下降。這可能與衰老時基因的表達水平改變相關,衰老時細胞內各類蛋白質的表達大部分下降。本實驗室前期研究發現,E2F1、POLD1等基因的表達水平隨著細胞衰老出現下降的趨勢[24～26],此類基因的啟動子區在年輕細胞中呈開放狀態,當細胞衰老時則表現得較為封閉。故此類基因的回帖序列經過富集分析,可出現衰老時轉錄起始位點的富集減少。但是,轉錄起始位點區域富集減少并非意味著衰老時細胞不出現轉錄調控。此時,細胞內以E2F1、POLD1為代表的一類基因還可以出現微弱表達,而另外一些基因(如P16INK4a等)則隨著細胞衰老的出現表達上升[27～28],這些表達上調基因的啟動子區會隨衰老的發生發展進一步開放,從而使轉錄起始位點富集的回帖序列增多。由于衰老時呈表達增多的基因仍較表達減少的基因少,所以從全基因組的層面分析,衰老時染色質的開放區域在轉錄起始位點區的富集下降。

其實,結合位點在全基因組功能區上的注釋,也進一步闡明了這一問題。對于年輕與年老的2BS細胞,其開放區段均富集于啟動子到轉錄起始位點區域,表明轉錄調控對復制性衰老有著重要影響。隨著細胞的逐漸衰老,啟動子區的開放比例進一步下降,說明當細胞衰老發生時,細胞內的轉錄調控水平減弱,而在年輕細胞中轉錄調控較為活躍。當然,在年老細胞中也存在開放的啟動子區,這是因為衰老的細胞也存在需要轉錄表達的基因。

此外,我們使用GO富集分析進一步驗證了轉錄調控對復制性衰老的調節。其中,涉及細胞進程、代謝、生物性調節、結合功能、催化活性等功能的基因,在細胞衰老過程中富集,此類基因與細胞自身調節、轉錄調控均有密切聯系。同時,我們也發現,作為復制性衰老經典標志物的P53、P16INK4a[1,29～30]等,均在 GO 分析結果中富集,如 P53富集于 GO：0048522、GO：0071840、GO：0016043等；P16INK4a富集于 GO：0071840、GO：0016043、GO：0044260等。

綜上所述,從轉錄調控探究、理解復制性衰老,應為我們后續研究的主要工作之一。

致謝

感謝上海嘉因生物科技有限公司提供的技術支持。