诱导多能性干细胞(Induced pluripotent stem cell,iPSC)因其具有类似胚胎干细胞的自我更新能力和多向分化潜能,在再生医学、体外疾病模拟和药物筛选等领域都有很好的应用前景
[1-2]。Oct4(Octamer-binding transcription factor 4)、Sox2(SRY-box transcription factor 2)、Klf4(Kruppel-like factor 4)和cMyc(Cellular myelocytomatosis oncogene)(简称OSKM)等体细胞重编程过程关键转录因子(Transcription factor,TF),在体细胞身份转变中发挥着十分关键的先锋调控作用
[3]。前期关于OSKM靶向调控体细胞基因组重编程的研究主要集中在蛋白编码(Protein-coding,PC)基因的启动子区域
[4-5],筛选出了一批在维持或诱导多能性状态中具有功能作用的蛋白编码基因,部分甚至可以在特定条件下替代原始的OSKM组合,例如用Nanog(Nanog homeobox)或Esrrb(Estrogen related receptor beta)替代cMyc
[6],或用Glis1(GLIS family zinc finger 1)等因子增强重编程效率
[7]。
长链非编码RNA(Long non-coding RNA,lncRNA)指长度大于200 bp的无蛋白编码功能的转录本
[8]。尽管lncRNA与蛋白编码mRNA在转录和加工机制上有诸多相似之处,但lncRNA通常缺乏有效的开放阅读框,因而不参与蛋白质的翻译过程。多项研究表明,lncRNA的表达比蛋白编码基因具有更高的细胞类型特异性
[9-10],这使lncRNA成为研究特定生物学过程中关键调控事件的重要切入点。近期,越来越多的证据表明lncRNA在体细胞重编程过程中具有同样重要的调控作用
[11-14]。然而,与蛋白编码基因相比,lncRNA的功能研究仍相对滞后,尤其是在重编程过程中,关键转录因子在 lncRNA基因组区域的靶向结合动态、结合位点(Binding site,BS)的选择偏好,以及结合是否与染色质状态和调控元件相关等方面,目前尚缺乏系统性的研究。这一空白限制了我们对非编码RNA在重编程调控网络中作用机制的深入理解。因此,亟需整合多组学数据,从转录因子结合图谱、染色质开放状态、调控元件特征等角度系统解析OSKM在lncRNA区域的调控模式。
本文围绕体细胞重编程过程中关键转录因子在lncRNA区域的结合特征与调控模式开展系统研究。作者基于公开的数据集,绘制了6种多能性相关转录因子的基因组结合图谱,并整合染色质可及性、组蛋白修饰及所结合基序的信息,揭示OSKM在重编程过程中呈时间依赖性地协同结合于特定lncRNA启动子区域,尤其偏好含内源性逆转录病毒K(Endogenous retrovirus K,ERVK)元件、较长转录本、与多能性基因邻近的长链基因间非编码RNA(Long intergenic non-coding RNA,lincRNA)。进一步构建PC-lncRNA共表达网络,筛选出可能参与iPSC多能性维持的候选lncRNA。
1 材料与方法
1.1 数据来源
本研究使用Chronis等
[15]上传至NCBI的小鼠体细胞重编程数据集(GSE90895)。包含小鼠胚胎成纤维细胞(MEF)、重编程48 h、pre-iPSC和iPSC四个阶段的RNA-seq和ATAC-seq数据,以及转录因子Oct4、Sox2、Klf4、cMyc、Nanog、Esrrb与组蛋白修饰(H3K4me3、H3K27me3、H3K79me2)及组蛋白变体H3.3的ChIP-seq数据。pre-iPSC和iPSC的划分依据Chronis等
[15]使用的Oct4-绿色荧光蛋白(Green fluorescent protein,GFP)报告系统。pre-iPSC指未激活内源性Oct4、GFP阴性的中间状态细胞;iPSC指重编程第10天GFP阳性、具有多能性特征的细胞。
1.2 转录因子ChIP-seq数据分析
ChIP-seq数据经FastQC、Trim Galore、Trimmomatic
[16]、Bowtie2
[17]和Samtools依次进行质控、剪切、比对及格式转换,使用MACS2(带宽150 bp,
q<0.005,富集倍数≥4)进行峰值探测,并借助RIdeogram
[18]可视化其基因组分布。
OSKM结合强度以log
2(RPKM+1)表示,采用
K-means(
k=7)对Oct4结合信号聚类。通过DeepTools2
[19]展示Oct4结合区域±2 kb内SKM信号分布,ChIPseeker
[20]分析其与转录起始位点(Transcription start site,TSS)的距离,启动子定义为TSS上游1.5 kb至下游0.5 kb。
1.3 染色质开放ATAC-seq数据分析
ATAC-seq数据与ChIP-seq数据一致,转录因子结合基序通过HOMER
[21]识别。
1.4 基因组lncRNA注释
小鼠lncRNA注释文件(GENCODE
[22] M21)包含位置、链别、ID与类型等信息。转座元件(Transposable element,TE)注释来自UCSC的rmsk文件,提取ERVK、ERV1、ERVL、L1等TE亚类,计算lincRNA启动子区域内的TE数量作为其TE含量指标。
1.5 基因功能富集分析
使用DAVID
[23-24]进行基因本体论(Gene Ontology,GO)分析,以气泡图展示结果,count表示富集基因数,
P.adjust为多重检验校正后的显著性水平,
P<0.05视为差异显著。
1.6 RNA-seq数据分析及共表达分析
RNA-seq数据处理同ChIP-seq数据。用Subread
[25]计数比对到基因组的短读片段。利用DEseq2
[26]分析各阶段间差异表达的lncRNA和PC基因,筛选倍数变化(Fold change,FC)≥2、错误发现率(False discovery rate,FDR)<0.05的基因。对iPSC多能性PC基因与差异lncRNA进行皮尔森相关性分析(R语言cor函数),筛选相关系数|
r|≥0.95、
FDR<0.01的基因对构建共表达网络。并用 Cytoscape
[27]进行可视化,借助CytoHubba
[28]和MCODE
[29]分别识别hub基因和关键模块。
2 结果与讨论
2.1 关键多能性转录因子的靶向结合染色体分布模式偏好
首先分析转录因子在染色体上的结合偏好,初步宏观探究其在基因组空间结构中的分布,以及重编程中染色体水平的调控差异,并进一步研究转录因子在启动子、增强子等功能区域的结合特征。为此,选取了Oct4、Sox2、Klf4和cMyc四种代表性因子,同时纳入Nanog和Esrrb(作为关键调控因子与OSKM协同作用
[15]),系统分析其结合分布及染色体偏好。
图1显示,这些转录因子主要富集于常染色体,特别是第9、10、13和17号染色体。相比之下,性染色体结合频率明显较低。值得注意的是,从图中可以观察到,转录因子ChIP-seq信号峰密度较高的位置,往往与启动子区域含有内源性逆转录病毒家族成员ERVK的lincRNA分布存在一定程度的重合。这一观察结果支持了Kelley和Rinn关于ERVK在小鼠胚胎干细胞维持中发挥重要作用的报道
[30]。
2.2 关键转录因子动态结合模式与染色质开放性关联分析
Yamanaka四因子中的Oct4是唯一不可被同家族成员替代的因子,在体细胞重编程过程中发挥核心作用
[31]。以核心因子Oct4为例,分析其全基因组ChIP-seq结合信号,发现结合位点可分为7类(
图2)。大部分位点(clusters 1~6)在重编程早期(48 h)即被占据,随后信号呈现增强或维持高水平结合(clusters 1、2、5)或减弱(cluster 6)趋势,少数位点(cluster 7)仅在iPSC阶段首次结合,显示Oct4按阶段性调控转录程序。值得注意的是,一旦Oct4的结合信号减少,几乎不会再有增强的趋势。其他SKM因子结合模式与其高度一致,表明四因子协同分批结合靶点。结合ATAC-seq数据显示,Oct4结合位点染色质开放程度与Oct4结合信号高度相关(
图2)。
2.3 lncRNA启动子区域关键转录因子的动态结合模式
启动子区域是基因表达调控的关键位点。我们对开放染色质区域的转录因子结合基序进行了富集分析,比较了lncRNA与PC基因启动子区域的差异。
表1显示,在iPSC阶段,lncRNA启动子中除常见CTCF 和Sp家族外,多能性因子(如Oct4、Sox2、Klf家族、Nanog)基序显著富集,其排序基于富集的统计学显著性(
P值)。重编程过程中,体细胞特异性转录因子如Fra1(Fos-related antigen 1)、Runx1(Runt-related transcription factor 1)、Tead4(TEA domain transcription factor 4))基序富集减少,而多能性因子的富集持续增强(
图3),反映调控网络向多能性转变。
启动子区域,lncRNA与PC基因的TF富集模式存在差异。lncRNA启动子中Oct4和Sox2基序在48 h富集后下降;PC基因中,Oct4基序从48 h起逐渐富集,Sox2则始终较低(
图3)。
依据Oct4在启动子区域的结合强度变化,将PC基因及lncRNA分为四类:1)持续结合(Persistent):48 h与iPSC的结合强度均大于等于0.66;2)结合减少(Lost):48 h的结合强度大于等于0.66,iPSC的结合强度小于0.66;3)结合增强(Gain):48 h的结合强度小于等于0.66,iPSC的结合强度大于0.66;4)无结合(None):48 h与iPSC的结合强度均小于0.66。发现SKM结合趋势与Oct4一致,且OSKM结合与H3K4me3、H3K79me2、H3.3等激活性标记呈正相关,与抑制性标记H3K27me3呈负相关(
图4),提示OSKM可能通过调控表观修饰促进基因激活。
GO分析显示,PC基因中Persistent组富集于细胞周期与RNA加工等核心生物过程,表明这些基因作为管家基因,维持着细胞的基本功能;Lost组则与分化相关,主要富集在肌肉细胞生长发育、细胞迁移等生物学过程中;Gain组主要富集在阳离子跨膜运输、减数分裂等生物学过程中;None组主要富集在应激反应和药物代谢等生物学过程中(
图5)。
2.4 体细胞重编程过程中OSKM结合启动子偏好与lncRNA特征的关系
为解析OSKM在lncRNA启动子区域的结合偏好,从几个lncRNA的特征进行了系统分析。
1) lncRNA类型
不同类型lncRNA启动子的OSKM结合频率存在显著差异(
图6)。与蛋白编码序列关系密切的lncRNA(如sense_overlapping、antisense、双向启动子型)被OSKM更频繁地结合;而如retained_intron和sense_intronic类型,它们位于内含子区域,远离外显子,其启动子对OSKM的结合频率显著较低。
2) lncRNA转录本长度
在iPSC阶段,较长的lncRNA转录本在启动子区域更倾向于显示出更多的OSKM结合位点,表现出一定的正相关趋势(
图7)。
3) 与多能性基因的距离
以lincRNA(数量最多,也最具功能潜力的一类lncRNA)为代表,分析iPSC阶段lincRNA启动子与多能性基因Oct4、Sox2、Klf4和Nanog(OSKN)TSS的距离发现:距离多能性基因越近的lincRNA,其启动子区域被Oct4结合的频次越高;特别是多能性基因TSS 10~20 kb内的lincRNA,其启动子结合Oct4水平显著高于平均水平(
图8)。这提示这些邻近的lincRNA可能通过顺式调控参与多能性建立与维持。
4) 启动子转座元件含量
转座元件(TE)是一类能改变自身在基因组中位置的重复序列,已被证实对于维持干细胞多能性具有重要作用
[32]。lincRNA中富含TE,其中ERVK对胚胎干细胞有极强的影响
[30]。据此我们假设,OSKM在重编程过程中对lincRNA启动子的结合偏好可能与其TE含量相关。为验证该假设,将小鼠lincRNA分为3类:ERVK-lincRNA,启动子区域含ERVK元件;OtherTEs-lincRNA,启动子区域含除ERVK外的其他转座元件;dTE-lincRNA,启动子区域不含任何转座元件。如
图9(a)所示,iPSC阶段ERVK-lincRNA启动子上的OSKM结合事件数量高于另外两组,表明OSK的结合偏好与TE,尤其是ERVK的存在密切相关。比较发现:ERVK-lincRNA转录本最长且拥有更长的第一内含子,见
图9(b);启动子TE数量与转录本长度呈正相关,见
图9(c)。这提示较长的转录本为TE提供了更多嵌入空间,也表明OSKM结合偏好可能更多受第一内含子结构影响。
2.5 小鼠iPSC多能性相关PC-lncRNA基因共表达网络分析
为挖掘在iPSC多能性维持中发挥关键作用的lncRNA,首先筛选出差异表达(|log
2FC|≥2,
FDR<0.05)的lncRNA和蛋白编码基因(Protein-coding gene,PCG)。结果显示在重编程48 h、pre-iPSC及iPSC阶段,上调的lncRNA数量普遍多于下调,见
表2及
图10(a)。随后将288个已知多能性相关蛋白编码基因(Pluripotency PCGs)
[33]与筛选出的4548个iPSC差异表达PCGs(iPSC-DE-PCGs)对比,发现有101个基因重合,见
图10(b)。这不仅在一定程度上验证了差异表达分析结果的可靠性,更重要的是,这些基因可能在iPSC多能性维持中发挥关键作用,我们将其定义为iPSC多能性相关蛋白编码基因(iPSC pluripotency PCGs)。
我们对iPSC多能性蛋白编码基因与iPSC中差异表达的lncRNA(iPSC-DE-lncRNAs)进行了共表达分析。发现91个iPSC多能性蛋白编码基因与210个iPSC-DE-lncRNAs存在显著相关(|
r|≥0.95,
FDR<0.01)。据此构建了由蛋白编码基因与非编码基因组成的共表达网络,共形成797条连接(
图11)。
我们进一步从共表达网络中筛选出关键枢纽(hub)基因,确定得分最高的10个hub基因,包括 Eppk1(Epiplakin 1)、Tdgf1(Teratocarcinoma-derived growth factor 1)、Dppa2(Developmental pluripotency associated 2)、Fb1(Fumonisin B1)、Phc1(Polyhomeotic homolog 1)、Pbx1(Pre B cell leukemia homeobox 1)、Nanog、Platr15(Pluripotency associated transcript 15)、AC132307.5和Stxbp2(Syntaxin binding protein 2)(
图12)。其中多个蛋白编码基因与多能性维持及重编程密切相关。如Tdgf1是Nodal信号通路的新调节分子,维持胚胎干细胞的多能性;Dppa2,这是近年来发现的在多能性细胞系中特异性表达的基因之一,可增强体细胞重编程效率并激活合子基因组激活相关基因
[34-36],并可能在小鼠胚胎干细胞状态的维持与增殖中发挥重要作用
[37];经典多能性因子Nanog也是我们鉴定出的hub基因之一,其在维持胚胎干细胞全能性中的关键功能已明确
[38-40];Phc1属于Polycomb抑制复合体,通过染色质重塑参与基因沉默,在调控多能状态及分化中扮演重要角色
[41]。
另一方面,我们发现两个高评分的lncRNA基因——Platr15和AC132307.5,也表现出典型的 hub特征。其中,Platr15属于Platr家族,该家族多个成员已被证实参与胚胎干细胞多能性的维持
[42-44]。值得注意的是,Platr家族的另一个成员Platr27亦在网络中排名第27,提示该家族成员可能共同参与多能状态的调控。
我们进一步识别出共表达网络中的关键模块,这些模块是共表达网络中基因联系最紧密的模块,其网络得分介于2.7至4.0之间(
图13)。模块中含多个多能性相关的蛋白编码基因,如Chaf1a(Chromatin assembly factor 1 subunit A),在维持多能性状态下的胚胎植入前稳定性中发挥关键作用,敲低会导致胚胎发育阻滞
[45];Mcm6(Minichromosome maintenance complex component 6)在DNA复制受阻时参与复制恢复,以维持基因组稳定性,对胚胎干细胞至关重要
[46-48]。
关键模块还富含与多能性密切相关的lncRNA,如Platr家族的Platr11、Platr27以及Panct(Pluripotency-associated noncoding transcript)家族的Panct2。已有研究表明,Panct1、Panct2、Panct3的表达与小鼠胚胎干细胞的自我更新及多能性状态密切相关,其功能涉及调控基因表达、染色质结构及细胞命运决定
[49-50]。具体而言,Panct2可能通过与核心多能性转录因子或表观遗传调控因子互作,参与维持干细胞的开放染色质状态和转录活性。
3 结论
通过整合关键转录因子的ChIP-seq与ATAC-seq数据,系统性地描绘了小鼠体细胞重编程过程中多个关键时间点上,Oct4、Sox2、Klf4、cMyc、Nanog和Esrrb等转录因子在全基因组范围内,特别是在lncRNA启动子区域的动态结合模式。研究发现,这些转录因子的结合呈现明显的时间依赖性,且与染色质开放状态高度相关,反映出它们在重编程进程中具有阶段性、协同性的调控特征。基于转录因子结合基序和组蛋白修饰数据的进一步分析显示,OSKM等转录因子可能通过与染色质修饰因子协同作用,共同调控lncRNA的转录活性。我们还发现,OSKM更倾向于结合具有特定特征的lncRNA启动子区域,如正义重叠和反义lncRNA,启动子含有ERVK、具有较长转录本或第一内含子,并且与多能性基因位于相邻区域。这些结果可能揭示了关键转录因子在lncRNA区域较为精细的调控图谱。更为重要的是,我们通过构建蛋白编码基因与差异表达lncRNA的共表达网络,鉴定出210个与iPSC多能性密切相关的lncRNA,初步勾勒出多能性调控中非编码RNA的网络架构。这些lncRNA可能通过多种机制——包括作为分子支架、调节染色质修饰酶的定位、影响转录因子活性等——充当多能性基因表达的关键调控节点。共表达网络中筛选出的候选lncRNA可能为未来深入功能研究提供线索。未来的功能验证研究,如基因敲除或过表达实验,将有助于揭示这些lncRNA在维持iPSC多能性中的具体机制。同时,这也可能为开发更高效的重编程策略和非编码RNA分子标记提供新的潜在靶点。综上所述,本研究可能为理解重编程过程中转录因子与非编码RNA协同调控机制提供理论参考,也为后续功能验证实验及重编程相关lncRNA分子标记的筛选奠定初步基础。