简介

转录因子(TF)与独联体-调控元件对基因调控至关重要。所有这些调控相互作用的全基因组图谱通常被称为囊区。它的动态使具有相似基因的物种和个体产生不同的转录程序，从而产生不同的表型，极大地促进了物种的适应性和表型可塑性1，2，3.．借助SELEX-seq、ATAC-seq、DAP-seq、ChIP-seq等多种高通量测序技术4，5，6在美国，现在可以在全基因组范围内绘制TF结合或识别其足迹，这在酵母和动物模型中引发了一系列跨物种TF比较。这些研究表明，具有保守生物学功能的同源tf在不同物种中共享的结合位点很少。例如，一项开创性的研究发现，与假菌丝发育相关的TF仅结合了三个相同靶基因的20%狭义敏感酵母菌物种7．对五种脊椎动物(人类、小鼠、狗、负鼠和鸡)肝脏组织中两个保守的肝脏bZIP和HB tf的ChIP-seq也发现，<10%的结合位点是保守的8．ENCODE项目最近进行的大规模TF分析也发现人类和小鼠之间TF足迹的保守程度很低9，10这表明动物池在进化过程中是高度动态的。

为了测量TF结合对转录的影响，在酵母中进行了详尽的ChIP-chip实验，以确定所有已知TF的结合位点11，12，13．令人惊讶的是，许多大多数tf结合位点没有明显的转录作用。在具有大基因组的高等真核生物中进一步的研究证实，tf可以结合到数量出乎意料的位点上，并且大多数tf似乎对附近基因的转录几乎没有影响，这表明转录调控网络具有较高的冗余性和系统鲁棒性9，14，15，16，17，18．有趣的是，在跨物种比较中发现的保守的tf结合位点通常对附近的基因表达有最强的影响8，19．此外，具有较强调控潜力的tf结合位点通常位于超增强子区，是多个共结合tf靶向的基因组热点区域15．因此，有人认为多个tf以定量的方式共同参与转录输出。单个的tf结合位点可能不足以解释转录，而增强子中的一组结合位点是在哺乳动物基因组中实现精确而强大的转录调控的关键20.．

尽管从动物和酵母模型中获得了丰富的数据，但对池的进化模式及其在植物中的相对重要性仍未充分研究。为了研究这一点，人们可以比较在多种植物物种中一个已知的TF与保守生物学功能的结合。光合作用可以说是植物中最重要和最保守的生物过程之一，GOLDEN2-LIKE (GLK)转录激活因子是控制叶绿体生物发生和发育的知名转录激活因子21．GLKs在水稻和玉米等单子叶植物中的亚功能化已被报道22．它们也在果实发育中发挥作用，其中GLK1被关闭，GLK2适应纬度梯度表达模式，导致果实组织颜色不均匀23，24．两份重复的GLK存在于迄今已测序的大多数二倍体被子植物基因组中25．他们的双功能丧失突变体拟南芥番茄和水稻叶片呈淡绿色表型，光合作用基因下调21，23，26．除了少数已知的靶基因外，GLKs的全基因组结合谱尚未确定，这使其成为跨物种湿地比较的理想候选基因。

在本研究中，我们利用ChIP-seq技术绘制了5种代表性植物中glk结合位点的图谱，发现大多数glk结合基因是种特异性的。保守的GLK结合基因通常与光合功能相关，其表达更容易受到GLK突变的影响。我们的研究结果揭示了植物进化过程中广泛存在的囊区分化和植物tf结合位点的冗余性。

结果

GLK结合的全基因组鉴定

为了研究植物贮水池是如何进化的，我们使用ChIP-seq技术确定了植物叶片组织样本中GLK1和GLK2的结合位点拟南芥（答:芥)、烟草(烟草benthamiana)、米(栽培稻)和玉米(玉蜀黍属可能)，以及番茄的叶片和未成熟的绿色果实组织(茄属植物lycopersicum)(图。1及补充表1)．之所以选择这五个物种，是因为它们可以转化为表达与表位标签融合的TF，并有参考基因组用于数据分析。我们使用了ENCODE2 ChIP-seq处理管道和MACS2-IDR算法27将通过QC临界值的库(NSC > 1.05, RSC > 0.9, FRiP > 1%，相关性> 0.8)用于后续分析(补充数据1)．

图1:植物进化过程中GLK TF结合的分化。
图1

ChIP-seq用于检测GLK1和GLK2的结合拟南芥烟草、番茄、水稻和玉米的基因组。鉴定了GLK的保守和特异的ChIP-seq靶基因，并对其特征进行了研究。

全尺寸图像

GLK被认为是一种控制光合作用和叶绿体发育的保守转录激活因子。我们首先检查了一些已知的光合作用相关基因拟南芥glk，比如叶绿素结合蛋白光收集系统a / b（LHCA / B)和光系统I/II亚单位的编码。我们发现，这些基因及其在其他四个物种中的同源基因都有5 '端近端GLK1和GLK2 ChIP-seq峰(图2)。2)．它们也与ATAC-seq峰重叠(图。2)，表明GLKs结合到其基因启动子中的开放染色质区域。总的来说，我们鉴定了960、1286、956、332和1089个GLK在近端启动子结合的基因拟南芥分别为番茄、烟草、水稻和玉米(图;1、补充数据2- - - - - -6)．

图2:ChIP-seq识别的GLK靶标。
图2

一个基因组浏览器跟踪显示GLK1和GLK2 ChIP-seq峰值，以及光合作用基因位点的开放染色质区域(ATAC-seq)。每个GLK的motif富集结果如下所示。bAtGLK1袋,k-mers机器学习模型ROC曲线。c前十名得分kAtGLK1模型中的-mers。d热图和平均信号图显示了番茄GLK结合位点附近的表观基因组特征。聚类使用叶片GLK1 ChIP-seq信号进行。从左至右依次为叶片中GLK1 ChIP-seq、叶片中GLK2 ChIP-seq、果实中GLK1 ChIP-seq、果实中GLK2 ChIP-seq、染色质可及性DNase-seq信号、H3K4me3、H3K27ac、CG、CHG和CHH位点的DNA甲基化。显示了基因转录起始位点上游和下游2 kb的区域。源数据作为源数据文件提供。

全尺寸图像

先前的基序富集和蛋白结合微阵列分析拟南芥GLK已表明其靶基因通常包含转录起始位点上游的GATTCC或rgatyyy基序28．我们提取了GLK ChIP-seq峰峰区序列并进行了motif富集分析。结果表明，RGATTYY核心基序确实在所有5个物种的glk结合位点中富集(图2)。2、补充数据7)．需要注意的是，motif的存在是TF结合的必要条件，但不是充分条件。DNA构象等特征可能影响TF在体外的结合，而DNA胞嘧啶甲基化、组蛋白修饰、染色质可及性和其他TF的共结合可能进一步使TF在体内的结合复杂化。因此，基因组中基序的数量往往远远超过实际结合位点的数量29．以GLK为例，在玉米基因组中有超过200万个GLK基序匹配(HOMER评分> 6)和40万个GATTCC完全匹配。为了克服基序搜索的局限性，建立了基于自然语言处理的机器学习模型来区分边界区域和非边界区域30.．这些模型的一个优点是，它是一个黑盒的集合，考虑了ChIP-seq峰内的多个短序列，其中可能包括TF的基序及其联合结合因子或其他影响DNA构象或表观基因组的序列。这种机器学习工具已经应用于为104种玉米tf生成精确的分类器模型6，30.．

因此，我们使用k-mer语法工具训练机器学习分类器来区分GLK结合区域和背景的DNA序列，背景是没有GLK结合的ATAC-seq区域。结果表明k-mer模型确实能够以较高的精度预测GLK在五个物种中的结合(图。2 b；补充数据8和补充图。1)．我们还提取了前10名k-mers在这些模型中的表达，并发现它们通常与GLK核心结合基序RGATTYY相匹配(图2)。2摄氏度及补充资料9)．例如，AtGLK1模型的受试者工作特征(ROC)曲线下面积为0.99，准确性评分为0.94，排名最高k-mer为GGATTTT(图;2 b, c)．如此高的准确性表明，这些模型既可以捕获识别基序，也可以捕获周围区域的隐藏序列信息，这些信息有助于tf结合的特异性。

为了了解围绕glk结合位点的染色质环境，我们将番茄GLKs ChIP-seq信号聚集在基因转录起始位点(TSS)上，并结合不同的表观遗传特征，如染色质可及性、DNA甲基化和组蛋白修饰(图2)。二维)．glk结合区域显示典型的活性染色质特征，如高染色质可及性和活性组蛋白标记H3K4me3和H3K27ac，通常标记启动子和增强子区域。他们还显示CG、CHG和CHH背景下的低DNA甲基化通常与tf结合位点有关。

先前的遗传学研究表明，这两个GLK基因在功能上是冗余的21，23，31．一致地，我们观察到GLK1和GLK2 ChIP-seq峰在光合基因启动子中重叠(图2)。2)．这两个glk有相似的主题k-mer富集(图;2及补充资料9)．他们的ChIP-seq信号热图也显示了相同的读分布模式。二维)，表明两个glk具有相同的结合位点。然后，我们比较了3472个番茄GLK1和6562个GLK2的峰值(IDR截止值0.01)，发现2815个重叠。即使在657个没有通过GLK2 IDR截止点的GLK1-only峰中，我们仍然可以检测到GLK2 ChIP-seq信号，反之亦然。2)．为了定量比较这两个slglk，我们还比较了它们在SlGLK1和SlGLK2峰并集中的ChIP-seq读覆盖，发现它们是相关的。在我们研究的所有物种中也发现了类似的结果(补充图。3.)．对仅有GLK1和GLK2峰的Motif富集分析表明，它们具有GATT核心Motif(补充图)。4)，表明不重叠的峰比重叠的峰弱，我们可以得出两个glk具有相同的结合轮廓。

GLK在光合作用相关过程中的保守靶标

为了了解GLK作为转录调控因子的功能，我们检测了其ChIP-seq靶基因的GO-term和MAPMAN注释，并进行了富集分析。引人注目的是，我们发现GLKs可以与叶绿体光合电子转移链中的大多数核编码基因结合，特别是光系统I和II(图2)。3及补充资料10)，以及叶绿素生物合成途径中的基因(图。3 b及补充资料11)．此外，GLK还可能在叶绿体和细胞核转录中发挥间接作用，因为它们始终与控制细胞器RNA加工的pentatricopeptide (PPR)和tetratricopeptide repeat (TPR)蛋白编码基因结合32，以及在所有五个物种中控制光周期的CONSTANS/B-BOX TF基因(图。3 c及补充资料2- - - - - -6）33，34．

图3:保守的GLK ChIP-seq靶基因。
图3

一个光合作用电子转移链示意图。在细胞核中编码的成分用紫色突出显示。方框内的数字表示该基因被glk限定的物种数量。对于番茄数据，仅使用叶片GLK靶基因。b显示叶绿素生物合成途径中的基因及其GLK ChIP-seq信号(−log(foldchange))的热图。c柱状图显示被GLK结合的保守TF基因的百分比。源数据作为源数据文件提供。

全尺寸图像

除了这些保守的GLK ChIP-seq靶基因外，我们还发现了物种特异性的靶基因。例如，在番茄和水稻中，GLKs只与编码PSII亚基S (PsbS)的基因结合，而叶绿体导入装置（CIA2)是由glk在拟南芥烟草(图;3和补充图。5)．但值得注意的是，ChIP-seq可以识别结合强度范围很广的TF结合位点。这种ChIP-seq峰的定性分析，将一个区域标记为绑定或未绑定，已知有局限性35，36．例如，在不同的ChIP-seq富集和峰值调用检测阈值下，弱结合位点在跨物种比较中容易被错误表征。为了定量比较两个物种之间的GLK结合，我们计算了它们同源基因对启动子中的ChIP-seq读取计数。保守的GLK靶基因对在两个物种中均具有较高的ChIP-seq信号，而种特异性的GLK靶基因对仅在一个物种中具有较高的ChIP-seq信号。6)．总之，我们的数据表明，尽管GLK的整体功能是保守的，但在分离这些物种的漫长进化过程中，它的结合已经出现了分歧。

GO项富集分析表明，每个物种的glk结合基因确实在光刺激响应(GO:0009416)和光收获(GO:0009768)等光合作用项上富集(图。4)．我们还使用MAPMAN注释将GLK ChIP-seq目标分配到不同的功能类别中(图2)。4 b)．光合作用类别最为丰富，尽管它们只占每个物种ChIP-seq目标的10%左右。虽然更多的GLK ChIP-seq靶点被分配到其他非光合作用类别，但没有一个被显著富集。这可能表明GLKs在祖先植物中与光合作用相关基因结合，这些结合位点处于强负选择作用下，而非保守位点在进化过程中获得或丢失。

图4:GLK ChIP-seq靶基因功能标注。
图4

一个各物种GLK ChIP-seq靶基因GO-term富集分析。bMAPMAN功能类别丰富分析(Fisher精确，双面)。右边的热图显示了每个MAPMAN类别的保护得分。源数据作为源数据文件提供。

全尺寸图像

大多数GLK ChIP-seq靶基因是种特异性的

在这五个品种中，长尾蛱蝶(拟南芥单子叶作物(玉米和水稻)的差异约为150 MYA。为了研究GLK结合是如何分化的，我们使用OrthoFinder将五个物种中的所有GLK结合基因分配到正交组中，并检查在每个物种中可以找到多少个共同的和独特的基因。结果表明，在所有5个物种中，只有极少数是保守的。例如，节点的129、100、151、205拟南芥GLK ChIP-seq靶点分别在5、4、3和2个物种中保守。其余375个基因在其他物种中没有直系GLK ChIP-seq靶点。我们还使用不同ENCODE2峰值调用管道在不同严格度下生成的基因列表重复了守恒分析。如果保守分析受到假阳性或弱结合位点的影响，当使用更严格的截止点时，保守率会增加。然而，当阈值提高时，识别的靶基因数量和保守率不受影响或下降，这表明我们的分析是稳健的(补充表2)．有趣的是，如果我们只计算MAPMAN指定的光合作用类别的基因，它们的保持率超过90%(补充表3.)，表明GLK结合的选择压力与基因功能相关，而与结合强度无关。

保守结合位点具有更大的转录调控潜能和更强的ChIP-seq信号

在动物中，研究表明最近获得的TF结合位点在调节基因表达的潜力方面往往不如那些保守的位点重要，尽管这些位点上的TF结合基序是相同的8．为了在植物中测试这一点，我们进行了RNA-seqGLK1/2双功能缺失突变体拟南芥和番茄，以及它们的野生型叶片作为对照(补充数据12，13)．我们鉴定了1105个差异表达基因(DEGs)拟南芥，其中266个为GLK ChIP-seq靶点，其中75.94%(202/266)在突变体中下调。然后我们检查了每个保存组中deg的百分比(图。5)．我们发现，在最保守的第5组中，超过一半(70/129)的基因是deg，而非保守基因的deg数量最少。对保守位点和物种特异性结合位点的motif -enrichment分析证实，它们都包含相同的RGATTYY motif(补充图)。7，8)．在番茄中也观察到这种模式，这表明保守的植物tf结合位点也具有比动物更强的转录调节潜力。

图5:保守的GLK ChIP-seq靶基因的表达。
图5

一个最保守的拟南芥GLK ChIP-seq的靶基因通常是差异表达的基因GLK突变体。根据目标基因的保守程度，将其分为五类。b保守基因组的GLK结合比非保守基因组强。的y-axis显示了峰值处GLK ChIP-seq信号的平均折叠变化(log 10) (Kruskal检验)。c火山图显示基因表达量(WT与突变体)的log2倍变化，ChIP-seq信号的log2倍变化。d-fGLK基序核苷酸多样性的累积分布函数(Ecdf)图的估计(K-S检验，单侧)。dglk绑定与未绑定的主题。e保守与非保守的主题。fdeg中的基序和非deg中的基序。源数据作为源数据文件提供。

全尺寸图像

我们假设非保守的GLK ChIP-seq靶基因最近通过启动子开放染色质区域的随机序列变化获得了GLK结合。其中一些独联体-调节元件尚未进化出有助于植物适应性的功能。如果没有选择压力，它们可能会逐渐消失。为了验证这一点，我们基于AtGLK1和AtGLK2 ChIP-seq信号在峰值处的平均折叠变化来估算TF结合强度(图2)。5 b)．结果表明，第5组中最保守的结合位点确实比种特异性结合位点具有更强的ChIP-seq信号(p-value = 1.1E−08)。接下来，我们结合结合强度和基因表达数据，发现保守基因(如光合作用功能注释的基因)通常具有更强的ChIP-seq信号，并且突变体中的基因表达比非保守基因有更大的降低(图。5度)．

非degs中的GLK结合位点也存在负选择

由于保守的GLK ChIP-seq靶点在突变体中更有可能存在差异表达，这是否表明非degs或非保守的靶基因中的结合位点是假阳性或没有生物学功能?首先，我们注意到超过40%的保守ChIP-seq靶基因不是deg，包括一些著名的光合作用基因(图2)。5)．例如，叶绿体ATP合酶亚单位基因ATPC1有很强的GLK结合位点拟南芥ChIP-seq信号分别排在第70位和299位(补充图。9)．然而,ATPC1只在番茄中被下调，而在番茄中没有拟南芥GLK突变。此外，也有具有保守和强结合位点的非deg。感光的隐花色素1在两个突变体中都没有差异表达，但在所有5个具有强ChIP-seq信号的物种中都是保守靶点(例如在拟南芥)．

我们还观察到保守位点的结合强度在不同物种之间也会发生变化。例如PSII亚基基因中的结合位点PsbTn在…方面很强拟南芥(第102位)，在番茄中较弱(第1083位)，尽管在两者中均有差异表达(补充图。9)．这表明，存在一个强的或保守的TF结合位点可能不能保证基因依赖于TF进行转录。其他tf可能以种特异性和位点特异性的方式弥补glk功能的缺失。或者，当TF被敲除时，TF靶基因可能在某些细胞类型或特定生长条件下有差异表达。因此，仅根据有限的基因表达数据往往很难推断单个结合位点的生物学功能。

为了使用不同的方法来评估tf结合位点的潜在功能，人们可以测量特定人群中该位点的核苷酸多样性。如果该位点处于负选择状态，则可能表明该位点对该物种的适应性很重要，尽管该位点附近的基因在TF突变体中没有差异表达。为了测试GLK，我们使用了1001拟南芥基因组重测序数据并计算其结合位点的核苷酸多样性。我们首先将ChIP-seq峰中的位点与背景进行比较，背景是在没有GLK ChIP-seq峰的开放染色质中发现的未结合的GLK基序(图2)。5 d)．我们还比较了保守位点和非保守位点(图2)。5 e)．正如预期的那样，它们的核苷酸多样性得分都低于背景，这表明它们处于负选择状态。然后我们比较了DEGs和非DEGs(图。5 f)，结果显示两者无显著性差异(Kolmogorov-Smirnov检验，p= 0.597)。这表明非degs中的结合位点也面临类似的选择压力，并可能以一种我们尚未理解的方式对植物的适应性做出贡献。因此，我们不应过度概括调控潜力与结合位点保护或结合强度之间的关系。

环流动力学和tf结合发散

鉴于tf结合的巨大差异，一个有趣的问题是它是否由基因组序列变异引起(独联体变异)。另外，TF蛋白的进化(反式-variation)可以改变其结合偏好。为了验证这一点，我们进行了转换拟南芥并进行ChIP-seq鉴定其在异源基因组环境中的结合位点。

我们首先划分拟南芥GLK ChIP-seq靶基因与非靶基因分为三组(I组:保守，II组:拟南芥-特异性和III组:玉米特异性)(图。6)．I组保守基因，如AtLHCA1而且ZmLHCA1，在两个物种中都与GLK结合。我们在其中发现了ZmGLK1 ChIP-seq峰值拟南芥基因启动子(图;6 b, c)．II组基因没有与GLK结合的玉米同源基因。例如，AtGLKs可以与叶绿素生物合成基因结合AtHEMA1，和它的正交正交向量ZmHEMA不是GLK的目标。但是，异源表达的玉米GLK1现在可以与它的启动子结合在相同的位置拟南芥GLKs(无花果。6 b, c)．最后，我们检测了玉米特异性的III组基因，如编码叶绿体ARM重复蛋白的AT3G03440。当ZmGLK1在拟南芥，它就不能再与启动子结合了。这种异源表达的ZmGLK1的ChIP-seq信号热图显示，它再现了肿瘤的生长过程拟南芥GLK在I组和II组基因中的结合模式，但在III组中未与玉米特异性基因结合(图2)。6)．综上所述，我们的数据表明，TF结合在很大程度上是由独联体通过DNA序列的变化而不是反式．

图6:异源表达的ZmGLK1与拟南芥基因组。
图6

一个三组ZmGLK1 ChIP-seq信号的热图拟南芥基因。b基因组浏览器跟踪显示第一组到第三组的代表性基因。c显示玉米同源基因的轨迹。

全尺寸图像

利用机器学习模型预测转录结果

为了更好地理解为什么只有部分具有近端glk结合位点的基因存在差异表达，我们试图定量比较它们的基因组特征，如tf结合强度、野生型叶片中的初始表达水平、结合位点与基因TSS的距离，以及从研究中推断出的其他tf共结合信息拟南芥TF DAP-seq数据采集。不出所料，差异表达基因(DEGs)具有更强的ChIP-seq信号，叶片中更高的基因表达水平，以及更紧密的glk结合位点(图。7一个)．我们还发现TF DAP-seq数据具有与DEG状态相关的共结合模式(图2)。7 b)．然而，必须指出的是，这些特征本身都不足以预测GLK ChIP-seq靶基因是否会有差异表达。因此，我们训练随机森林分类器模型，一起考虑不同的定量基因组特征。为了避免202个deg的小数据集的异常值的影响，我们随机采样500次训练和测试数据来训练500个独立的模型。模型的平均AUC、ACC、Recall和F1-score分别为0.77、0.698、0.685和0.688(图2)。7 c)．这表明，结合不同类别的基因组特征确实可以改善我们对基因组调控密码的解释。

图7:预测目标基因表达的机器学习模型。
图7

一个差异表达和非差异表达的GLK ChIP-seq靶基因的不同定量基因组特征拟南芥GLK双突变体(Wilcoxon检验，双面，n= 328)。bDAP-seq tf共结合特征与GLK靶标DEG状态的Pearson相关系数直方图。c500个随机森林模型的ROC曲线。d箱线图显示500个随机森林模型去除不同特征后的ROC曲线下面积(AUC)。方框包含数据集的第25 - 75个百分位数，中线表示中位数，胡须表示1.5 IQR。e模型的前10个DAP-seq TF特征，根据它们在500个随机森林模型中的平均重要性得分进行排名。源数据作为源数据文件提供。

全尺寸图像

这种方法还使我们能够测试哪些特征对模型的准确性是重要的。我们发现，去除TF共结合数据会导致模型性能的最大下降，这表明共结合TF可以保存预测转录结果的重要信息(图2)。7 d)．随机森林模型还允许我们计算DAP-seq数据中每个联合绑定tf的特征重要性得分。我们发现一些众所周知的光合作用和光信号调节因子，如myb相关，C2H2和GBF tf是主要贡献者(图2)。7 e)．虽然相关性并不意味着因果关系，但有可能这些差异表达的基因在祖先植物中已经被GLKs和这些tf共同调控，并且它们已经进化成为基因调控网络中不可或缺的一部分。这与在动物基因组中观察到的情况类似，具有较强调控潜力的TF结合位点往往位于被大量TF结合的增强子中15．

追踪基因组复制后GLK的结合发散

番茄基因组经历了古代eudicots γ复制(~120 MYA)和近代的γ复制(~120 MYA)茄属植物谱系一(~70 MYA)。由于tf结合的启动子区域也在这些事件中重复，因此它提供了一个独特的机会来估计GLK与其靶标之间的相互作用是在其中一个重复之前还是之后发生的。我们发现番茄的glk与一对编码基因结合光收集装置a4（LHCA4)在第3和第6染色体的同位块中。两种基因在突变体中均下调(图2)。8)．这两个同义块具有较高的同义替换率(KS = 1.689)，并且不位于先前报道的由最近的基因组三倍性引起的亚基因组区域37，表明GLK绑定在LHCA4在经过古复制后保存了超过100亿年。然而，在基因组复制后，基因可能会丢失或重新排列，导致大多数番茄基因不再位于同步性块中。重复的基因仍然可以被识别为GLK靶向的同源基因对。我们计算了5个保守组中GLK ChIP-seq同源靶点的百分比，发现最保守的第5组确实有更多的同源靶点基因，这表明保守的结合位点是古老的，而物种特异性的结合位点是最近获得的。

图8:复制的GLK ChIP-seq靶基因。
图8

一个基因组浏览器跟踪显示了含有两个番茄LHCA4基因的部分同步性块。所示为开放染色质、GLK ChIP-seq和RNA-seq轨道。bGLK ChIP-seq目标基因在各保护组中重复的比例。c玉米GLK ChIP-seq靶基因在亚基因组1和亚基因组2中表达模式的Pearson相关性基因组复制后保留两个GLK结合位点的配对具有最高的相关性(t以及,片面的)。源数据作为源数据文件提供。

全尺寸图像

与番茄不同，玉米基因组最近经历了5-12 MYA的四倍体事件38．它的大部分复制基因仍然可以追溯到亚基因组1和亚基因组2。然后我们用Kendall秩相关系数检验比较了重复的ChIP-seq靶基因的GLK ChIP-seq信号。Kendall的tau为0.174 (p= 0.06322)，说明相关性很弱，重复的GLK结合位点在如此短的进化周期后已经发生了变化。我们还假设，TF结合的快速变化可能允许复制的基因在不同的组织中表达，这可能进一步导致新功能化。为了验证这一点，我们利用玉米GLK ChIP-seq的组织基因表达数据计算了两个亚基因组中玉米GLK ChIP-seq靶基因对的表达相关性。结果表明，保留GLK结合的基因对表达相关性最高，而失去GLK结合的基因对表达相关性显著降低，证实了我们的假设(图。8 c)．

讨论

GLK结合的跨物种比较表明，植物胞囊动力学也可以引起广泛的tf结合变化，就像在酵母和动物基因组中观察到的那样9，14，15，17，18．我们的观察结果也符合先前对SEPALLATA3 MADS-box TF的ChIP-seq分析答:芥而且答:lyrata他们发现，只有不到四分之一的结合位点是保守的39．尽管非保守和弱TF结合位点对基因表达的影响较小，但有人建议我们不应将其视为无功能位点14．因为真核生物的基因往往是由多个tf调控的，组合输入可能超过转录激活的阈值，导致整个基因调控网络产生冗余。此外，TF结合位点的调控潜力不能通过单个转录组分析完全证明，因为非差异表达的TF靶基因可能在另一种生长条件或组织中出现差异表达。我们的数据一致地显示，在非差异表达的GLK ChIP-seq靶基因中，GLK结合位点也与差异表达的GLK靶基因中的GLK结合位点存在相似的负选择。因此，研究tf如何共同调控基因表达对于全面理解基因组调控密码具有重要意义。

许多植物TF研究，特别是那些难以转化的作物，通常是基于突变互补或模式植物(如番茄和番茄)的过表达来推断基因功能拟南芥．我们的发现使我们有理由谨慎地解释这类研究的结果。TF结合是否受基因组影响独联体-调控动力学，人们可能错过了在外源基因组中表达的TF的物种特异性功能。例如，乙烯不敏感3类tf是保守乙烯信号通路的关键调节因子40．我们之前已经发现，一些更年期果实在成熟基因启动子中获得了EIN3结合位点，这些启动子已经进化成为控制成熟的正反馈回路中不可缺少的部分41．因此，虽然这些水果中的EIN3基因可以补充拟南芥EIN3功能丧失突变体和在过表达时产生相同的表型，它们在果实成熟中最近进化的作用无法在植物中得到检验拟南芥基因组。

我们的研究的一个局限性是，如果没有chip级抗体，我们必须过表达与表位标签融合的TF来进行ChIP-seq，这使得很难控制蛋白水平以及TF的表达时间和地点。由于ChIP-seq是一种富集分析，蛋白质丰度和细胞位置、抗体的特异性、测序深度、文库复杂性等因素，甚至ChIP实验本身的小细节，如输入抗体比、超声处理和洗涤步骤，都可能影响最终的ChIP富集结果。即使在富集程度较低的情况下，也可以很容易地检测到强tf结合位点。但是鉴定出的弱tf结合位点的数量很容易受到ChIP-seq富集的影响，使得跨物种或跨实验比较变得困难35．例如，我们之前已经证明，与传统的ChIP-seq方法相比，灵敏的chip沉着试验可以用相同的抗番茄MADS-box TF RIN抗体检测出5倍多的结合位点41，42．随着CUT&Tag和单细胞ChIP-seq技术的发展，其灵敏度有可能进一步提高。此外，基于神经网络模型或不同方法的集成的下一代峰值调用算法现在能够从嘈杂的ChIP-seq数据中找到弱tf结合位点43．灵敏度的提高意味着未来可以发现更多的弱结合位点，这对数据分析提出了重大挑战。

另一个有趣的问题是，为什么基因组会保留这么多遗传上冗余的TF结合位点，而没有强大的调节作用。据推测，获得调控复杂性是具有复杂细胞类型的多细胞生物进化的关键，通过使共同基因被多次利用来产生不同的时间和空间转录程序和生物学结果44．随着湿地的快速进化独联体-调控元件通过序列变异产生或中断，导致新的TF招募或现有TF结合减弱，从而产生不同的转录结果和表型。因此，观察到的植物物种间GLK结合的差异，以及其保守的核心调控相互作用，可能是进化过程的实时快照。

在酵母和动物模型中的一些研究表明，即使是较弱的tf结合位点也可以在进化新的基因调控程序中发挥重要作用，并且是实现高tf结合特异性所必需的45，46．例如，用人工的高亲和hox结合位点取代天然的低亲和hox结合位点导致基因异位表达果蝇胚胎47．最近使用高灵敏度ChIP-seq的植物TF研究表明，每个TF可以与大约10,000个位点结合5，48．的拟南芥cistrome项目发现529个tf的270万个体外结合位点。如果一个TF可以与基因组中三分之一的基因结合，那么所有的结合位点在调控转录方面就不太可能是平等的。因此，人们可能需要定量和合作地考虑TF结合，正如我们在使用随机森林模型预测GLK ChIP-seq靶基因表达时所演示的那样。此外，最近对胚胎干细胞重编程因子Oct4、Sox2、Nanog和Klf4的分析表明，TF结合的间隔和方向也可以编码关键的调控信息49而对植物tf的类似分析在很大程度上已经落后。

方法

植物材料

的拟南芥glk1 / glk2突变体25是从欧洲获得的拟南芥库存中心(N9807)。转基因拟南芥在添加25 μg/mL湿霉素B和0.8% (w/v)琼脂的Murashige和Skoog (MS)盐上培养。植物在22°C的条件下(光照16小时)在生长室中生长。野生型选用粳稻品种日本光。水稻植株在28°C的条件下长时间在温室中生长。番茄glk1 / glk2在microTom背景下，由CRISRP/Cas9产生的双突变体由Yu Pan教授提供，我们通过Sanger测序确认了缺失(Supplementary Fig。10)．番茄和烟草(烟草benthamiana)在光照12h(25°C)和暗12h(20°C)循环的生长室中生长。

为拟南芥glk，我们将它们的cdna与c端GFP标签融合并将它们转化为glk1 / glk2ChIP-seq双突变体(补充表1)．番茄的转基因表达GLK2会引起共同抑制吗23．因此，我们将SlGLK2-GFP构建物转化到番茄品种Micro-Tom中，该品种在番茄中存在功能丧失突变GLK2轨迹。由于番茄glk也参与了早期果实发育，我们还使用开花后27天的未成熟果实进行了ChIP-seq。在泛素启动子的控制下，将水稻GLKs的cDNA与HA标记融合转化为野生型水稻。对于这些稳定的过表达GLKs的转基因植株，我们选择了蛋白表达水平较高的转基因株系进行ChIP-seq。烟草是一种新近的异体四倍体，每个GLK有两个拷贝，序列同一性超过90%。因此，我们选择了叶片表达量最高的两个GLK基因进行分析。我们使用农渗法瞬时表达带有HA标签的NbGLKs用于ChIP-seq。玉米(玉米GLKs ChIP-seq数据来自之前的研究6．玉米的GLK1-HA结构也转化为拟南芥glkChIP-seq双突变体。

ChIP-seq

完全展开的叶子拟南芥，米饭，西红柿，和烟草benthamiana以及27 DPA下未成熟绿番茄果皮样品，用1% (w/v)甲醛交联。然后通过过滤和离心分离细胞核，并使用Bioruptor超声检测染色质到亚kb片段。ChIP使用抗gfp (#A-11122, ThermoFisher)和抗ha (#C29F4, Cell Signaling Technology)抗体。将A/G Dynabeads蛋白(10 μL)与2 μL抗体在低盐缓冲液(20 mM Tris-HCl pH 8.0, 150 mM NaCl, 2 mM EDTA, 0.1% Triton X-100, 0.1% BSA)中孵卵1 ~ 2 h。珠子被洗了洗与低盐缓冲(20毫米Tris-HCl pH值8.0,150毫米氯化钠,2毫米EDTA,特里同x - 100 1%, 0.1% SDS)前增加了用染色质,一夜之间被孵化与Dynabeads抗体,并与低盐洗两次洗涤缓冲(20毫米Tris-HCl pH值8.0,150毫米氯化钠,2毫米EDTA,特里同x - 100 1%, 0.1% SDS),高盐缓冲(20毫米Tris-HCl pH值8.0,500毫米氯化钠,2毫米EDTA,特里同x - 100 1%, 0.1% SDS),氯化锂缓冲区(20毫米Tris-HCl pH值8.0,250毫米氯化锂,1% Triton X-100, 0.7%脱氧胆酸钠，1 mM EDTA)和TE。免疫沉淀的DNA用TS-Tn5标记在珠上45 min。然后用高盐缓冲液和TE洗涤两次。样品在65°C下反向交联过夜。用Qiagen MiniElute纯化DNA, PCR扩增。 The libraries were sequenced with a Hiseq X 150 bp paired-end read mode.

RNA-seq

总RNA用RNeasy Mini Kit (Qiagen)提取。信使rna用Oligo d(T)25磁珠(New England Biolabs)分离，用于Illumina TruSeq文库制备。文库在Hiseq X上以150 bp的端对方式测序。

ChIP-seq数据分析

ChIP-seq reads与植物参考基因组(拟南芥TAIR10,烟草benthamiana2.5.1版本,茄属植物lycopersicumSL4.0,栽培稻MSU7.0,玉米RefGen_v4)使用Bowtie 2(版本2.3.2)。对于配对端，读取150 bp，从3 '端开始的100 bp用bowtie2参数−3 100修整。SAMtools(版本1.9)使用参数“-F 4 -q 20”过滤未映射和低映射质量的读取。使用SAMtools中的子命令“rmdup”删除重复读取。使用deepTools中的multiBigwigSummary命令计算Pearson相关系数，以评估生物重复的可重复性。仅保留相关系数为> 0.8的生物重复进行进一步分析。使用带有MACS2(版本2.2.1)和IDR(版本2.0.4.2)的ENCODE2 ChIP-seq管道进行峰值调用。使用MACS2中的“- call- peaks”函数检索峰值位置。将距离峰顶位置±75 bp的窄峰格式信号折叠变化区域提供给IDR，保留超过IDR 0.01的区域，并将峰顶信号折叠变化截止点调整回150 bp。然后，它们与基于TSS峰顶距离的基因相关联。 Depending on genome size, the ChIP-seq peak summit overlaps with the putative promoter region 1.5, 2.5, 2, 1.5, and 2.5 kb upstream of the TSS were associated with genes in拟南芥分别是烟草、番茄、水稻和玉米。我们还使用了PhantomPeakQualTools(版本1.14)的ENCODE2 spp-IDR-TIP管道在两个不同的截断点(p= 0.05和p= 0.01)进行比较。

RNA-seq分析

RNA-seq reads首先被映射到rrna上，然后干净的reads被HISAT2 (version 2.1.0)映射到相应的植物基因组上。用HTSeq (version 0.11.0)计算读取计数。用DESeq2进行差异基因表达分析。

功能与母题丰富分析

基于web的农业社区工具包agriGO v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/)．使用MAPMAN进行基因功能类富集。利用OrthoFinder (version 2.2.7)比对5个物种的目的基因，根据蛋白质序列相似性确定同源关系。我们使用保守评分来量化一组基因的保守程度(例如，在一个MAPMAN类别中，eudicot共享GLK靶基因)，基于事件(观察到的)和执行的成对比较(预期的)。采用HOMER (http://homer.ucsd.edu/homer/motif/)．在核苷酸多样性分析中，我们首先确定了GLK基序的位置拟南芥GLK ChIP-seq峰值使用HOMER(截止:motif匹配得分> 6)。我们还在没有GLK ChIP-seq信号作为背景的开放染色质区域(ATAC-seq峰)中鉴定了未结合的GLK基序命中。每个基序位置的核苷酸多样性分数从1001所识别的变体集合中检索拟南芥基因组计划(https://1001genomes.org/data/GMI-MPI/releases/v3.1/)．

机器学习分析

我们使用k-mer语法工具(https://bitbucket.org/bucklerlab/k-mer_grammar/src/master/)生成机器学习模型，以区分GLK结合位点序列与没有GLK ChIP-seq峰的开放染色质区域序列。对于glk结合位点，我们测试了以ChIP-seq峰值为中心的150和300 bp序列。以ATAC-seq峰顶为中心的相同长度序列为背景，去除ChIP-seq峰弱的开放染色质重叠(IDR 0.05)。由于开放染色质区域通常比ChIP-seq峰值多10倍，因此随机选择等量的开放染色质区域以生成平衡的数据集用于训练k- m模型。平衡的输入数据集随机分为80%的训练集和20%的测试集k-mer语法工具。

为了研究glk结合位点的不同转录调节功能，我们还使用随机森林分类器来区分下调的DEGs和非DEGs。用于RF模型训练的特征是基因最近的ChIP-seq峰顶的GLK1和GLK2 ChIP-seq信号折叠变化值，峰顶到基因转录起始位点在碱基对中的距离，基因在FPKM中野生型叶片中的初始表达水平，以及拟南芥TF DAP-seq从植物池数据库(http://neomorph.salk.edu/dap_web/pages/index.php)．使用DAP-seq峰与150 bp GLK ChIP-seq峰区重叠的GLK靶基因。将DAP-seq数据转换为分类特征(1个绑定，0个未绑定)进行训练。特征-标签相关性低(<0.1)的DAP-seq TF数据集被丢弃。由于非deg比deg多，我们随机选择等量的非deg标签来匹配deg，得到一个平衡的训练数据集，80%的数据用于训练，其余20%的数据用于测试。在样本量较小的情况下，我们也随机重新选择500次训练数据和测试数据，生成500个模型来计算平均AUC和模型得分。

并列分析和K年代的计算

MCScanX在默认参数下进行共线性和共线性检测。同义替换率(Ks)由KaKs_Calculator(2.0版本)计算。

报告总结

有关研究设计的进一步资料，请参阅自然组合报告摘要链接到这篇文章。