单细胞测序可以揭示单个细胞的微妙的变化并对细胞进行分型，在经历了动物单细胞的火热后，植物单细胞转录组的春天来了，Science, PC, PP, MP, bioRxiv各一个。本文详细解析了2019年4月，Brad Nelms和Virginia Walbot两位研究人员在Science期刊上通过单细胞转录组测序解释玉米减数分裂进程的工作。
文章能发到Science不只是因为测序的植物，更重要的是分析过程的创新和完善的思路，以及很好的描述。Science文章还是有很多值得学习的地方，看透一篇，可以用来分析很多工作，值得仔细阅读。
解读：这个小姑娘真棒
校对：生信宝典
研究背景
减数分裂使生物细胞中染色体数目减半并且通过同源重组产生新的等位基因组合，对有性生殖个体的生命过程起着非常关键的作用。植物细胞减数分裂调控对于种子育种及农作物产量具有非常重要的意义。然而植物细胞进入减数分裂的机制尚不清楚。玉米的花药细胞在缺氧条件下会诱导体细胞分化为孢子细胞，为减数分裂过程的起始细胞。经过3天的级联有丝分裂过程，孢子细胞停止有丝分裂成为花粉母细胞，进入减数分裂前期。减数分裂前期和早期的细胞形态和基因表达发生改变，但是其在减数分裂染色体减半阶段发生的时间和关系尚不明确。在前减数分裂孢子细胞分化过程中是否存在特定基因表达的中间阶段？基因表达大幅改变是否就意味着减数分裂前期的开始? 本文就是为了解决这两个问题。
文章简介
在多细胞生物体中，减数分裂是一个复杂的过程。本文利用单细胞RNA测序的方法（scRNA-seq）重塑了玉米雄性细胞进入减数分裂的发展过程。在进入减数分裂前的细线期，有一个两步的连续的基因表达转录组变化的过程，在这个过程中，有26.7%的基因转录丰度发生了至少两倍的变化。对细胞周期基因表达的研究表明几乎所有的前生殖细胞都会增殖，而非依赖干细胞模型产生减数分裂细胞。生殖细胞的转录程序是细胞自发的，即使减数分裂失败，转录过程仍然继续进行。
文章亮点
提出了根据“拟发育时间变化速度”进行定量评估的方式，根据基因表达的急剧改变来指示发育过程的转换，并且利用这个框架来定义生殖细胞发育过程的中间状态。本文为通过单细胞转录组研究植物细胞发育阶段提供了新的指引。
结果分析
1 scRNA-seq的样本取材
（1）玉米生殖细胞的发育过程
玉米雄性生殖细胞来源于未成熟的花药，一般在每一个花药叶中（图1A）。花药在早期发育中延伸，因此，花药长度可以作为发育阶段和器官年龄的指标。生殖细胞在减数分裂前期为有规律并且连续的发育变化的。
（2）取材部位：
生殖细胞一般比花药细胞的直径大2~4倍，因此，在组织裂解后，可以根据大小将生殖细胞分离出来。为了减小样本取材过程中引起的基因表达的改变，分离得到的细胞将在1~2小时冻存。在细胞裂解后，样品分装到两个试管，分别进行细胞分离，作为技术重复（图1B）。
（3）测序分析：
从24个植株中分离到144生殖细胞，横跨了孢子细胞分化到减数分裂前期的粗线期的一周的发育时间，平均每个细胞检测到101,245个转录本（图1C）。技术重复产生相似的转录本数目。64%转录本变化由组织阶段引起的。而仅有1%的变化由细胞分离时间引起的，提示组织阶段起重要作用。
图1. 前减数分裂和早期减数分裂细胞的scRNA-seq信息 (A) 玉米花药早期的发育简图，单细胞从1.3-1.5 mm的花药细胞进行分离。 (B) 细胞分离和单细胞文库构建的流程图 (C) 主成分分析表明细胞在第一主成分的分布与花药长度正相关。 AR， 孢元细胞；PMC，花粉母细胞；PC1，第一主成分
2 生殖细胞分化早期的动态过程
考虑到聚类分析忽略了发育的动态过程信息，并且不能反应内部结构的变化，为了阐明生殖细胞分化的中间过程，研究者通过拟合发育轨迹，推测细胞所处的发育时间，从而定量基因表达随时间变化的规律。
首先通过降维来拟合时间，拟合发育时间定义为在推测的发育轨迹上单个可以代表细胞之间最大基因表达变化的变量（图2A），并且两次技术重复拟合出的结果一致性为0.97，说明结果的可靠性。其次，将样品通过拟合的发育时间来排序，并且计算相邻细胞的拟合发育时间变化的斜率代表细胞之间表达谱的变化强度（图2B）。拟合时间速率曲线表明发育起始时基因表达有一个连续稳定的改变，随后有三个速率明显高于平均值的峰（图2C）。这些峰与8个单细胞聚类群的边界相一致，但是拟时间速率仍然保留了聚类分析所忽略的细节。例如，拟发育时间变化速率峰的高度提供了一种定义表达改变幅度的方法。在前半个时间进程中拟发育时间变化速率保持在平均速率的两倍以下，提示有一段连续的基因表达变化过程，而这个过程是聚类方法所不能呈现的信息。为了探索拟发育时间速率改变的原因，鉴定了3046个随时间显著变化的基因，用热图可视化它们的表达。热图重复了拟时间速率鉴定的特点，首先是平滑的连续的基因表达改变，随后有阶梯状的离散改变。
为了探索拟发育时间变化速率与已知的发育进程的关系，通过拟时间聚类细胞并且评估它们来源的花药细胞的长度（代表不同的发育阶段）。基因表达连续改变阶段（蓝色、灰色区域）包含细胞从 ≤1.1 mm花药（孢子细胞）到早期的花粉母细胞阶段。表达没有明显的改变提示孢子细胞停止有丝分裂进入花粉母细胞。基因表达在孢子发育过程中并不是静态的, 1/3的拟时间过程发生在这个阶段。然而在孢子数目上没有明显的改变反映了快速的分化事件的发生。在第二个拟时间速率峰后，所有的细胞均来源于≥1.3 mm花药。总体来说，减数分裂第一阶段开始于1.2 mm的花药，提示最后两个峰发生在减数分裂的早期。下文将称这些表达改变为前期过渡表达阶段1（Pr1）和前期过渡表达阶段2（Pr2）。估计这两个转变阶段发生在10小时之内。
图2 从渐变和离散的细胞状态转变来阐述生殖细胞分化 (A) 基于前10个主成分拟合发育轨迹获得拟发育时间曲线，技术重复的一致性很高在0.97。 (B) 细胞按拟发育时间排序，10个细胞为一个单位计算变化斜率代表拟发育时间速率。 (C) 拟发育时间速率与生物发育时间的关联。表达变化速率是相对于平均速率而言。不同的颜色区域表明不同的取材部位。 (D) 热图展示所有随拟发育时间显著变化的基因。
3 减数分裂阶段基因表达的改变
减数分裂前期两个转录水平的改变是否与减数分裂过程中细胞所处的分裂阶段有关？分裂阶段的定义基于染色体的细胞学形态。从scRNA-seq数据的Pr1和Pr2阶段选择了5个标记基因（图3A），两个在Pr1和Pr2中下调（Ago18a 和Rpl38e）, 两个在Pr1中上调（Rmf 和 C3h3）,和一个在Pr2中上调（Trps8）。然后比较从46个花中分离的1.1到1.6 mm的花药基因表达与染色体细胞学形态的关系。每朵花的同步发育的3个花药，一个固定用来观察细胞状态，另外两个用来分离减数分裂细胞用于qPCR定量。
通过qPCR定量检测标记基因的表达发生改变，这与单细胞转录组测序数据是一致的（图3B）并且与细胞学发育进程相吻合。Pr1的标记基因（Rmf 和 C3h3）在减数分裂前的间期和细线期表达丰度低，在细线期和中央核仁阶段呈现上调，在后面的细胞学阶段都呈现较强的表达水平。得出结论，Pr1转变发生在细线期和中央核仁阶段。Pr2的标记基因（Trps8）在偶线期的所有花中都上调，而在偶线期之前没有变化，提示Pr2高表达时完全进入了偶线期。这些结果表明形态学和基因表达提供了所有细胞阶段的互补的视角，它们并不完全重合（图3C）。在Top 5000表达的基因中，有26.7%的基因在Pr1和Pr2期发生两倍或更高的表达水平的变化（图3D）。
为了探索Pr1或Pr2表达变化是否受控于体细胞或减数分裂过程中染色体形态的变化，检测了Pr1和Pr2的标记基因是否在体细胞分化或减数分裂进程相关功能缺失突变体中的表达。结果发现除了Pr2的标记基因Trps8延迟表达外，其他的基因均正常表达，表明这些基因的调控并不受来自邻近的体细胞的影响。进一步检测了am1的两个等位基因突变体am1-praI和am1-489中的表达（am1-praI突变体可以进入减数分裂前期，但是停滞在偶线期；am1-489突变体在生殖细胞中进行孢子的有丝分裂而不是减数分裂）。结果发现所有的Pr1和Pr2标记基因均在am1-praI突变体中表达（图3G）。标记基因在am1-489突变体的2 mm的花药阶段（正常花药I期的后期）异常表达，而在2.5 mm花药（正常减数分裂的四分体时期）阶段标记基因的表达量与野生型相同（图3H），其染色体形态在2.5 mm是仍然保持异常。因此，这些突变体中，早期阶段标记基因的表达与减数分裂染色体的形态是相解离的。得出的结论是Pr1和Pr2状态的转变是细胞自主发生的，在正常的减数分裂受阻后，其仍能正常进行。
图3 在细线期转录组的重构 (A) Pr1和Pr2的标记基因的表达展示。发育阶段处于Pr1和Pr2时期的细胞标注为红色。 (B) 基因表达状态与细胞所处分裂阶段的比对。花药根据核中是否存在致密的染色体，核仁的定位和knob的形态分为4个细胞学形态。 (C) 在细线期早期细胞形态学和基因表达之间的关系。 (D) 在Pr1到Pr2阶段和孢子细胞到花粉母细胞的表达变化≥2倍的基因所占的比例。（E）-（H）上部为2 mm花药的生殖细胞形态展示。底部为标记基因相对于1 mm花药的表达情况。
4 核质协调
为了探索在生殖细胞分化过程中被调控的基因的功能，将所有表达发生改变的基因进行无监督聚类并探索它们富集的生物进程（图4A）。减数分裂核心功能（染色体重组和联会）基因分别富集在上调表达的聚类簇5和6，每个减数分裂基因的上调时间与已经鉴定到起计时做用的编码蛋白的功能之间关联较弱（图4B），编码在减数分裂前期发挥作用的蛋白的基因（例如内聚蛋白亚基Afd1）比在后期作用的基因（例如编码联会复合物蛋白Zyp1）提前上调。也有一些在减数分裂细胞相比于周围体细胞差异表达的基因富集在聚类簇4,5,6（图4C），减数分裂细胞在减数分裂期间的时间和空间表达并不总是一致（图3）。
随后用AgriGO对不同簇的基因进行了富集分析，持续下调的基因富集到代谢和基本细胞生物学相关的通路（图4A）。簇5和6 的基因富集到特定的细胞器，如线粒体、内质网膜和高尔基体等。同时蛋白翻译相关的基因表达下调，占到两个下调的基因簇2和3中基因的20%-30%。核糖体基因的表达降低了32.6%，在细胞分裂早期急剧降低24.3%。这些结果与已经证实的植物减数分裂前期细胞的重塑（例如，增加线粒体的数目和体积及其它的异质膜结构）一致。因此，减数分裂表达反应了染色体重塑和减数分裂胞质分化相关的转录调控。
图4 生殖细胞分化的转录调控总结 (A) 主成分分析散点图展示6个共调控基因簇的表达情况。 (B) 在减数分裂前期第一阶段已经确定功能的基因的表达情况。 (C) 在性母细胞中表达高于周围体细胞的基因在基因簇中所占的比例。 (D) 在GO数据库中注释为蛋白翻译的基因在每个基因簇中所占的比例。星号表示显著性差异。 (E) 编码胞质核糖体蛋白亚基的基因的累计表达水平。
5 有丝分裂和减数分裂之间的细胞周期调控
由于孢子细胞减数分裂的不同步，利用群体转录组很难解释减数分裂期间基因表达调控，单细胞方法则可以解决这个异质性问题。转录组数据研究发现有特定细胞周期阶段的基因簇（图5A），因此可以根据这些数据推断出每个细胞所处的细胞周期（图5B）。预测处于S阶段的细胞的比例与之前对特定花药长度范围细胞所处的细胞周期阶段的评估一致（图5C）。单细胞测序鉴定出更多的细胞周期阶段，提供了一个更加完整的细胞周期过程。(既突出结果准确性，又强调新发现。)
是否存在减数分裂和有丝分裂特异的细胞周期基因表达调控？在减数分裂前期的有丝分裂增值阶段，发现细胞周期阶段相关的基因存在周期性表达的模式（图5B）。而有丝分裂细胞周期的多个阶段的基因在减数分裂中的都有表达基（图5B）。例如，减数分裂细胞表达几种组蛋白亚基(在原始细胞S期表达)与微管结合蛋白(表达在原始细胞G2/M期)，细胞周期蛋白和细胞周期蛋白依赖的激酶的组合表达并未出现在原始细胞中。这种现象有利于区分减数分裂前期的第一阶段和有丝分裂细胞期。通过计算每个基因与(i)拟发育时间和(ii)最相似的细胞周期的相关性来鉴定原始细胞 （花药细胞长度≤1mm）发育调控基因，这一分析可以将原始细胞分化基因与细胞周期调控基因区分开（图5D），而如果使用整个花药的数据则没有办法进行这种分析。9种细胞周期蛋白选择性的表达在细胞周期的特定阶段（图5E）。但是有一小部分与细胞周期无关的基因也表现出细胞周期阶段特异的表达调控。
图5 有丝分裂和减数分裂的细胞周期调控基因的表达 (A) 细胞周期调控基因的主成分散点图，每个点代表一个基因，颜色代表细胞所处的周期。PCA主成分分析实战和可视化 附R代码和测试数据(B) 细胞周期调控基因在代表性细胞的特定阶段的表达。 (C) 不同花药长度群体中预测的细胞周期阶段的比例 (D) 基因与细胞周期和拟发育时间的关联分析 (E) 编码细胞周期蛋白和细胞周期依赖的蛋白激酶的基因的表达热图。
总结
在原始细胞增殖和有丝分裂停止的转换阶段，转录组图谱变化稳定平缓没有发生巨变。随后转录组出现两个急剧变化过程。通过比较转录组变化与染色体细胞形态学的关系，发现第一个转变发生在细线期，第二个转变发生在偶线期。研究结果表明减数分裂期间转录表达的改变不仅与核事件相关，而且与细胞形态相关。核糖体转录丢失和膜结合细胞器转录增加，减数分裂的胞质重塑与其将二倍体转变为单倍体有关。
看完了不过瘾，下面还有很多单细胞分析教程等您，更有易生信的单细胞培训班期待您的参与，自己学会分析，才能做更多特色和创新。
