人多能干细胞（human pluripotent stem cell, hPSC）系之间的变异性仍然是生物医学中的一项挑战，同时也是一种机遇。我们鉴定出不同 hPSC 系在早期谱系出现过程中的差异，这些差异可映射到胚胎发育的前后轴上。RNA-seq 分析揭示了动态的转录组模式，这些模式界定了在不同 hPSC 系间保守的中胚内胚层与神经外胚层谱系的出现，以及在整个分化过程中保持不变的细胞系特异性转录特征。
稳定的细胞系特异性转录组模式可预测各细胞系对视黄酸（retinoic acid, RA）的响应，从而导致其在前脑与后脑命运之间表现出不同的偏向。重复建立的 hPSC 系以及与之配对的成人供体组织表明，来源于不同个体的人类细胞表达独特且持久的转录组特征，这些特征与进化上较近出现的基因相关。除这一遗传贡献之外，我们还发现，来自同一供体的重复细胞系表现出不同的脑区命运，这与不同的染色质状态相关，提示表观遗传机制也会促成神经命运差异。
在大量 hPSC 系中，也观察到了这种谱系偏向的变异及其与 RA 响应性基因表达的相关性。这些结果界定了 hPSC 中可启动决定前部或后部神经命运这一关键早期步骤的转录组差异。
在哺乳动物胚胎中，上胚层中数百个多能细胞通过在空间上协调一致的细胞状态限制，生成不同的组织（
满足这一需求的紧迫性可由近期研究加以说明：这些研究报道了不同 hPSC 系在产生区域特异性神经前体方面存在发育差异，以及这些差异在神经发育障碍病因学中的潜在意义（
近期研究已经界定了多种具有谱系形成能力的状态，这些状态在无约束条件下可产生不同的胚内与胚外命运（
此前尚未以系统性方式检验这种空间组织在不同 hPSC 系之间的情况。为测量这些自组织过程在不同 hPSC 系之间的变异，研究将 hESC 系 SA01 与 hiPSC 系 i04 进行了比较（
多条细胞系之间早期神经命运指定的变异性，通过三条 hESC 系（H9、UC06、SA01）和三条 hiPSC 系（i04、i07、i13）在自我更新（self-renewal, SR）和神经外胚层（NSB）条件下 SOX21 蛋白表达的持续性差异得到了证明（
为更全面地分析这些细胞系及条件下低维转录组变化，我们采用了 Genome-Wide Coordinated Gene Activity in Pattern Sets（GWCoGAPS）非负矩阵分解算法（
为了将这些在体外界定的转录组动态与体内发育过程中谱系出现时的基因表达变化联系起来，我们对发育中小鼠原肠胚的单细胞 RNA 测序（single-cell RNA-seq, scRNA-seq）数据进行了投射（
GWCoGAPS分解还识别出各细胞系中特异性的表达模式，这些模式在时间和处理条件变化下均保持不变（
在神经外胚层起始模式P15中，SOX21具有较高排序（
为进一步检验SOX21在调控中内胚层与神经命运差异性出现中的作用，在第2天（D2T）对SOX21-KO细胞进行BMP4处理，此时自我更新细胞中的核心区域已经形成。与在SR最初几天对BMP反应性的早期限制一致（
SOX21在前部神经外胚层早期特化中的这一作用，还得到如下观察的支持：在NSB条件下SOX21 KO细胞中上调的基因，在发育中的小鼠上胚层和前原条中表现出更高的表达（
为进一步探究细胞系特异性转录组特征如何与谱系出现相关，我们采用多维尺度分析（MDS）显示动态GWCoGAPS模式与细胞系特异性GWCoGAPS模式之间的相关性（
为进一步探索这些细胞系中前后轴模式化的差异，我们检测了已知受视黄酸（RA）上调的基因表达（
细胞系特异性表达模式中RAR基因的差异性富集提示，这些模式可用于预测各个细胞系在前后轴方向上的分化行为。我们通过监测已知驱动后脑细胞命运的HOXB1表达，考察了多个细胞系在神经分化过程中的RA反应（
为进一步检验这些人类细胞变异测量指标的有效性，我们从3位供体（2053、2063和2075）的头皮成纤维细胞中新生成了6株hiPSC细胞系；这些供体的死后脑组织RNA-seq数据此前已由我们获得（
最初的基因型-组织表达（GTEx）联盟数据调查了来自多位供体的多种组织中的基因表达和遗传变异，并鉴定出跨组织eQTL（
这种对个体人类间基因表达变异的遗传控制可能是进化上较近时期才出现的。为检验这一假设，我们通过系统发育分层法估算了基因的进化年龄（
此外，还估算了基因剂量敏感性（
除3种供体特异性特征外（
为界定其神经分化轨迹中的变异，我们将这些新细胞系的RNA-seq数据投射到GWCoGAPS-I神经外胚层模式P15和P12中
为排除大规模体细胞突变或重编程相关基因组突变导致2053-2与2053-6细胞系之间观察到的分化偏倚的可能性，我们对这些重复细胞系以及供体2053和2075的成熟死后脑组织进行了基因组测序。我们相对于参考基因组鉴定了这些基因组DNA样本中各自获得和丢失的变异（
为了寻找这两组细胞系中这种假定表观遗传谱系偏倚背后的共同转录组要素，我们比较了 i04 和 SA01 细胞系特异性特征中的基因权重差异，以及 2053 供体特异性与 2053-6 细胞系特异性特征之间的差异。GBX2 在 i04 和 2053-6 细胞系中显示出较高的相对表达，这两者均偏向于后部神经命运（
与可能影响稳定基因表达表型的表观遗传机制相关的是，我们发现，在我们鉴定出的所有细胞系特异性和供体特异性转录组特征中，KRAB-ZNF 基因均显著富集（补充表 S5）。已有研究表明，KRAB-ZNF 基因在早期发育过程中抑制转座元件的表达，并构建由 H3K9me3 介导、持久存在的异染色质，以调控成熟组织中的基因表达（
为进一步探究这一点，我们在 SR 条件下对 hiPSC 细胞系 2053-2 和 2053-6 进行了 H3K9me3 ChIP-seq 数据采集。将这些数据投射到 2053-6 细胞系特异性转录组特征中显示，在 2053-6 细胞系中过表达基因的启动子区域（TSS 上下游 ±3 kb）中，2053-2 细胞系富集了这种抑制性表观遗传标记（
为在更大规模的人群中进一步探索谱系偏倚中的遗传和表观遗传差异，我们对来自 NextGen Consortium 的 101 名供体所生成的 317 个未分化 hiPSC 细胞系的 RNA-seq 数据进行了主成分分析（PCA）（
将 NextGen PC1 投射到小鼠早期发育图谱上表明，hPSC 细胞系中的这种变异包含显著的谱系偏倚（
除高 PC1 中具有中胚内胚层偏倚的 hPSC 细胞系中 RAR 基因高表达外，在中等 PC1 水平的一簇细胞系中也观察到 RAR 基因高表达（图中星号所示
类似地，为检验细胞系特异性特征是否也与更广泛的 hPSC 变异相关，我们将 NextGen PC1 的基因载荷与显著偏倚的 i04 和 SA01 细胞系稳定转录特征之间的差异进行了比较（
在干细胞研究领域，预测单个 hPSC 细胞系的分化能力一直是长期关注的重点。越来越多的研究正在揭示 hPSC 在转录异质性和分化潜能方面异质性的遗传起源（
此外，我们观察到来自同一供体的重复细胞系在前脑与后脑神经命运之间存在偏倚，这提示表观遗传机制也对早期神经命运选择的变异作出重要贡献。我们的数据支持这样一种模型：在多能性状态内部，存在若干替代性的、具有表观遗传预倾向的状态，它们先于实施这些早期区域性神经命运选择的高水平转录组输出而存在。这一机制与近期在小鼠原肠胚中的发现相呼应，即在由转录组定义的细胞身份形成之前，已存在表观遗传层面的结构化过程（
我们在此界定的自组织转录组架构，使得能够在经典组织者细胞类型完全表达之前，对hPSC内区域性胚胎特征进行映射。我们此前已报道，沿端脑背腹轴存在一种区域性神经命运偏倚；这一偏倚在来自同一供体的hiPSC细胞系中也存在差异（
与保守的动态表达模式相反，我们观察到，那些在不同hPSC细胞系之间系统性不同的模式，同样富集于进化上较新的基因。我们的结果显示，在细胞系特异性特征中存在KRAB-ZNF基因富集以及H3K9me3介导的调控，这提示了控制这些稳定表达表型的可能机制。同样受H3K9me3机制调控的、进化上较新的转座元件及其对应的KRAB-ZNF抑制基因，已被直接证明参与人类合子基因组激活以及着床前状态，并在发育后期指定基因表达（
神经发育障碍的遗传学（
近期在为干细胞治疗选择特定hPSC细胞系方面取得的进展（
使用Accutase（A11105，Life Technologies）将hPSC解离为单细胞，并以1 X 10的密度接种
hiPSC细胞系i04、i07和i13（NIH-i4、NIH-i7、NIH-i13）此前已有报道（
SOX21-KO hESC细胞系通过CRISPR/Cas9介导的基因组缺失系统构建。使用CRISPR Design Tool，即MIT优化的CRISPR设计工具，为Sox21NHEJ4设计SOX21特异性gRNA（
细胞以4%多聚甲醛固定10 min，并使用PBS中0.1% Triton X-100（Sigma-Aldrich）透化40 min。随后，细胞用10%驴血清（Sigma-Aldrich）封闭，并与一抗孵育过夜。
所使用的一抗及其稀释比例如下：BRACHYURY（AF2085，R&D，1:500）、CDX2（AM392，Biogenex）、EOMES（ab23345，Abcam，1:400）、GATA3（MAB6330，R&D，1:200）、GATA4（AF2606，R&D，1:400）、GBX2（AF4638，R&D，1:200）、HOXB1（AF6318，R&D，1:200）、HOXB4（ab133621，Abcam，1:400）、HOXB9（ab66765，Abcam，1:400）、ID1（AF4377，R&D，1:200）、ISLET1（AF1837，R&D，1:200）、NANOG（AF1997，R&D，1:200；Reprocell 1:200）、OCT4A（MAB17591，R&D，1:200）、OLIG2（AF2418，R&D，1:200）、OTX2（AF1979，R&D，1:200）、PAX6（PRB-278P，BioLegend，1:500；AF8150，R&D，1:200）、PHOX2B（AF4940，R&D，1:200）、p-SMAD1/5（9516，Cell Signaling Technology，1:200）、p-SMAD2/3（8828，Cell Signaling Technology，1:200）、TUJ1（PRB-435P，BioLegend，1:1000）、SOX1（AF3369，R&D，1:400）、SOX17（AF1924，R&D，1:500）、SOX2（AF2018、MAB2018，R&D，1:200）、SOX21（AF3538，R&D，1:200）、SOX3（GT15119，Neuromics，1:200）以及TUJ1（MAB1195，R&D，1:400）。
二抗孵育采用Alexa Fluor偶联抗体，稀释比为1:400（Life Technologies）。对于直接免疫染色，使用Alexa Fluor单克隆抗体标记试剂盒（A20181、A20184、A20186，Life Technologies）对一抗进行偶联。细胞核采用DAPI（Life Technologies）复染。
图像使用Operetta（Perkin Elmer）采集，在Columbus服务器（Perkin Elmer）上利用自定义模块进行批处理分析，并使用Spotfire（Perkin Elmer）进行可视化。
hPSC上皮的空间分析（“距边缘距离”测量）通过在Columbus中运行自定义Acapella脚本（Perkin Elmer）实现，具体命令如下：1）使用20×物镜，将用户定义的3×3相邻重叠视野拼接为蒙太奇图像；2）对单个细胞核的DAPI信号进行分割和二值化，以生成核对象；3）对每个细胞核周围细胞质中的DAPI信号进行分割和二值化，以生成细胞质对象。hPSC表现出强蓝色荧光，这是由于视黄酯在细胞质脂质小体中的隔离所致。
按照制造商说明书，使用 mirVana 试剂盒（Ambion）提取总 RNA。采用 Agilent 2100 Bioanalyzer System 进行 RNA 质量控制。按照制造商说明书，使用 Illumina mRNA sequencing sample Prep Kit（用于 Poly-A 文库）或 TruSeq Stranded Total RNA RiboZero sample Prep Kit（用于链特异性文库）构建 RNA-seq 文库。
简而言之，从约 800 ng 经 DNase 处理的总 RNA 中纯化含 poly-A 的 mRNA 分子，或使用 RiboZero 磁珠去除核糖体 RNA。纯化后，在 94°C 高温条件下利用二价阳离子将所得 RNA 片段化 2 min，形成小片段。片段长度范围为 130–290 bp，中位长度为 185 bp。采用逆转录酶和随机引物将切割后的 RNA 片段复制为第一链 cDNA。随后使用 DNA Polymerase I 和 RNase H 合成第二链 cDNA。
这些 cDNA 片段随后依次经过末端修复处理，所用酶包括 T4 DNA polymerase、T4 PNK 和 Klenow DNA polymerase；再使用 Klenow exo（3’→5’ 外切酶活性缺失）加上单个“A”碱基；最后利用 T4 DNA Ligase 连接 Illumina PE 接头。如有必要，可在 Illumina 接头中插入索引，以便在一个 8 通道流动槽的单个通道中对多个样本进行测序。
RNA 浓度通过 Qubit（Life Technologies）测定。RNA-seq 文库质量采用 LabChipGX（Caliper）并结合 HT DNA 1K/12K/HiSens Labchip 进行检测。最终的 cDNA 文库使用 HiSeq 2000（用于 Poly-A 文库制备）或 HiSeq 3000（用于 RiboZero 文库制备）进行高通量 DNA 测序。
测序完成后，采用 Illumina Real Time Analysis（RTA）模块进行图像分析和碱基识别，并使用 BCL Converter（CASAVA v1.8.2）生成包含序列读段的 FASTQ 文件。测序深度超过 8000 万条（4000 万对双端）可比对测序读段（补充表 S1）。采用 FastQC（v0.10.1）进行读段水平的质量控制（Q/C）。
cDNA 序列的双端读段基于 Ensembl Build GRCh37.75 的已知转录本，使用剪接读段比对软件 TopHat（v2.0.4），在默认参数并加入“--mate-innder-dist 160”选项的条件下，比对回人类基因组（Illumina iGenome 中的 UCSC hg19）。对于链特异性的 RiboZero 样本，TopHat 使用“--library-type fr-firststrand”选项。
采用 samtools（v0.1.18）和 RSeQC（v2.3.5）完成比对统计和质量控制，以对所得比对读段计算质量控制指标，分别提供关于可比对性、基因体覆盖均一性、插入片段长度分布以及剪接位点注释的有用信息。
为了获得基因水平的表达谱，使用“samtools sort –n”选项获得正确成对且成功比对的读段，并根据基因注释（Illumina iGenome）使用 htseq-count v0.5.3 对这些读段进行计数（采用 intersection-strict 模式，且对 RiboZero 样本采用 stranded 选项），随后计算 RPKM。由此获得了 23,368 个基因水平的表达谱。
为确定来源细胞系对核内 SOX21 蛋白水平的影响
主成分分析使用 R 中的 prcomp() 函数完成。基于 RNA-seq 数据中基因水平 RPKM 的基因凝聚型层次聚类采用 R 中的 hclust() 和 cutree() 函数，并使用相关距离（dist=1-r）进行。GWCoGAPS 采用默认参数运行，如前所述（
基因年龄采用系统发育分层法进行估算，该方法利用 BLASTP 评分的蛋白质序列相似性来确定蛋白编码基因的最小进化年龄（
细胞在室温下于持续振荡条件中用1%甲醛交联10 min，随后用125 mM甘氨酸淬灭5 min。在自我更新培养的第2、4和6天收集细胞核。使用微球菌核酸酶（MNase）对染色质进行片段化，直至大多数DNA片段长度处于200–700个碱基对范围内。染色质与H3K9me3抗体（13969，Cell Signaling）在4°C下持续振荡孵育过夜。抗体结合的染色质在4°C下于持续振荡条件中用Protein G Dynabeads进行1 h免疫沉淀。
将蛋白-DNA复合物从Dynabeads上洗脱，并在0.3 M NaCl条件下于65°C孵育4 h进行DNA逆交联。使用QIAquick Gel Extraction Kit（28706X4，Qiagen）纯化游离DNA。采用TruSeq ChIP Library Preparation Kit（Illumina，IP-202-1012）构建ChIP-seq文库，并在Illumina HiSeq 3000平台上进行测序。
使用Burrows-Wheeler Aligner（BWA）将ChIP-seq原始reads比对至人类基因组组装hg19，并使用MACS2.1.1进行峰值识别。为了解H3K9me3在启动子区域的占据情况，使用bedtools统计了转录起始位点（TSS）上下游±3 kb范围内的ChIP-seq reads。
S.K.、S.S.、G.SO.、A.J.、J.G.C.、D.J.H.、C.C.和R.D.M.构思了本研究。S.K.、S.S.和Y.W.进行了细胞培养和分化。Y.W.构建了iPSC细胞系。K.O.构建了SOX21-KO细胞系。S.K.、A.J.、Y.W.和K.O.生成了RNA-seq数据。S.S.生成了ChIP-seq数据。G.SO.、S.S.、E.J.F.、C.C.和J.S.开发并应用了信息学方法以分析RNA-seq数据。V.L.和A.K.进行了进化基因年龄分析。
S.S.和C.C.分析了ChIP-seq数据。T.M.H.、J.K.和D.R.W.提供了用于iPSC细胞系构建的成纤维细胞和人脑组织数据。S.K.、S.S.和N.M.进行了免疫细胞化学实验。S.K.、S.S.、A.J.、T.V.和D.J.H.进行了高内涵图像分析。S.K.、S.S.、G.SO.、N.M、J.G.C.、D.J.H.、N.S.、C.C.和R.D.M.解释了数据。N.S.、R.B.、A.J.C.、N.J.B.、D.R.W.和R.D.M.指导了本研究。
S.K.、S.S.、C.C.和R.D.M.撰写了手稿。所有作者均参与了结果讨论和手稿修改。
我们感谢Lieber家族和Maltz家族对Lieber Institute for Brain Development（LIBD）该项工作的慷慨支持。本研究还获得了授予E.J.F.的NCI/NIH项目资助R01CA177669、P30CA006973、U01CA212007和U01CA253403以及约翰斯·霍普金斯大学Catalyst Award的支持，并获得了授予N.S.的R01NS116418、R01HG010898、MH116488和U01MH124619资助。
通过NeMO Analytics进行的数据共享与可视化获得了R24MH114815和R01DC019370项目资助。我们感谢J. Jessee和MTI-GlobalStem在CRISPR/Cas9质粒转染方面提供的技术支持。我们感谢LIBD和Sestan实验室的众多成员对本工作的宝贵意见。
手稿已根据新的发现和补充数据进行了修订。
📄 原文链接：https://www.biorxiv.org/content/10.1101/2021.03.17.435870
🏷️ 人多能干细胞 神经分化 视黄酸响应 转录组分析 表观遗传调控