人多能幹細胞（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
🏷️ 人多能幹細胞 神經分化 視黃酸響應 轉錄組分析 表觀遺傳調控