灵长类动物中G-四联体形成序列保守性的模式研究

《Journal of Molecular Evolution》：Mapping Patterns of G-Quadruplex-Forming Sequence Conservation in Primates

　　摘要
虽然非编码G-四链体（G4s）在基因启动子和剪接位点处作为保守的调控元件发挥作用，但蛋白质编码区域中G4s的进化保守性仍鲜有研究。为了解码G4s的进化动态，我们绘制并分析了22种灵长类动物同源编码DNA序列（CDS）中的潜在G4形成序列（PGQS）。研究发现，PGQS的数量与可用的同源基因数量相关，而在一定程度上也与与人类的系统发育距离相关。PGQS基序在亲缘关系密切的物种中表现出高度共定位。从最小折叠自由能推断出的热力学稳定性成为与进化模式相关的重要因素：低稳定性的PGQS（最小折叠自由能≥–10 kcal/mol）更具保守性，而高稳定性的PGQS（最小折叠自由能≤–30 kcal/mol）的保守性较低，尽管这两类序列的保守性仍高于CDS基线。一致地，插入/缺失（indel）评分与PGQS的最小折叠自由能呈负相关，这表明稳定的基序与插入或缺失事件之间存在关联。与此一致的是，富含G的串联重复序列表现出较高的插入/缺失突变率，这与它们倾向于折叠成高度稳定的G4s相符。总体而言，这些发现表明PGQS同时作为保守元件和结构不稳定性的来源，反映了在保持序列功能的同时通过结构形成产生不稳定性的拮抗性选择压力。

引言
G-四链体（G4）是由富含鸟嘌呤的序列形成的非典型二级结构，可以在DNA和RNA中折叠（Lombardi和Londo?o-Vallejo 2020）。G4的典型基序由共识序列G3+N1?7G3+N1?7G3+N1?7G3+描述，其中G3+代表三个或更多连续的鸟嘌呤，N1?7代表1到7个任意类型的核苷酸环（Lombardi和Londo?o-Vallejo 2020）。四个鸟嘌呤通过Hoogsteen键相互作用形成平面G-四联体，这些G-四联体堆叠起来构成四链G4结构（Burge等人2006；Gellert等人1962；Sen和Gilbert 1988）。G4结构还通过单价阳离子（主要是K+和Na+）进一步稳定，这些阳离子位于结构的中心空腔中（Bhattacharyya等人2016；Tateishi-Karimata和Sugimoto 2014）。然而，也存在偏离典型基序的结构变异，包括长环、凸起和错配（Jana等人2021；Mohanty等人2025；Mukundan和Phan 2013；Palumbo等人2009；Papp等人2023；Vannutelli等人2022, 2023；Varizhuk等人2017）。在DNA中，G4结构与转录调控、DNA复制、染色质重塑和端粒维持等基本过程相关，从而影响基因组稳定性和基因表达的控制。在RNA中，G4结构主要参与转录后调控机制，调节翻译调控、miRNA和piRNA处理、前miRNA剪接、前mRNA多聚腺苷酸化、mRNA运输和线粒体转录等过程，并有助于维持端粒稳态（Bhattacharyya等人2016）。DNA中的G4形成通常需要局部解开双螺旋，而RNA分子不受互补链的限制（Bhattacharyya等人2016；Cheong和Moore 1992；Kim等人1991）。因此，RNA G4结构可以更容易形成，并且通常比DNA G4结构具有更高的热力学稳定性（Bhattacharyya等人2016；Cheong和Moore 1992；Kumari等人2007；Sacca等人2005）。计算分析表明，G4在人类基因组的重要区域（如复制起点、端粒区域、核糖体DNA、免疫球蛋白重链类别转换重组区域以及多个基因和癌基因的转录调控区域）富集（Maizels和Gray 2013）。此外，已经在体内证明了G4的存在（Biffi等人2013），高通量测序方法的使用使研究人员能够在人类（Chambers等人2015）和其他生物的基因组中实验性地绘制G4（Marsico等人2019）。使用基于抗体的G4染色质免疫沉淀技术研究了人类染色质中G4的分布（H?nsel-Hertsch等人2016, 2018）。这些研究表明，人类基因组中约有10,000个区域可以有效地折叠成G4；此外，这些G4结构在基因启动子处高度富集，并与增强的转录活性密切相关。RNA G4已在人类基因中通过计算方法绘制出来，并在转录起始位点（Huppert和Balasubramanian 2007）、5′-UTR（Huppert等人2008）和第一个内含子的5′端（Eddy和Maizels 2008）富集，在编码区域则减少（Maizels和Gray 2013）。最近的研究表明，由于密码子偏倚选择，人类编码基因组中的稳定G4显著减少（Mirihana Arachchilage等人2019）。编码G4的删除可能是由于它们对翻译的负面影响（Benhalevy等人2017；Endoh等人2013；Endoh和Sugimoto 2016）。特定稳定G4基序（G3N1）4的基因组范围代表性不足似乎取决于所分析的系统发育群体：虽然热力学上非常稳定的G4基序（具有相同的单核苷酸环组成G、C或T）倾向于被抑制，但(G3A1)4基序在哺乳动物中尤其是灵长类动物中得到强烈保留（Lombardi等人2019）。此外，作者还表明，编码区域中的遗传不稳定性增加与G4的高热力学稳定性直接相关。稳定G4基序的保守模式表明，积极的G4相关生物学作用得到了维持，而潜在的有害效应在进化过程中被消除。多项研究表明，稳定G4结构在调控非编码区域（主要是基因启动子）的出现在不同生物中是保守的（Capra等人2010；Marsico等人2019；Rawal等人2006）。所使用的方法基于使用计算和/或DNA测序方法在不同基因组中识别潜在的G4形成序列（PGQS），然后比较多个物种在特定基因组位置的PGQS共现模式。结果是得到了跨物种的详细非编码PGQS地图。在编码区域内，进化力量既反对（Eyre-walker和Bulmer 1993；Kudla等人2009；Mirihana Arachchilage等人2019）也支持（Guiblet等人2021；Katz和Burge 2003）稳定二级结构。研究表明，从编码DNA序列（CDS）衍生的mRNA转录本中形成的稳定G4结构可能导致核糖体停滞，而影响G4稳定性的沉默突变可以增强蛋白质表达（Agarwala等人2015；Endoh等人2013）。这表明G4折叠/解折叠动态在编码区域具有功能意义。在系统发育框架内对灵长类动物进行比较基因组分析对于阐明人类遗传结构和灵长类多样性至关重要（Shao等人2023）。对灵长类动物CDS的进化分析使我们能够理解反映作用于该群体的进化力量的一般保守和分化模式（Hellmann等人2003）。尽管G4对突变、细胞功能和临床相关性很重要，但它们的进化仍研究不足（Mohanty等人2025）。特别是，对CDS区域中G4的进化保守性的研究很少，可能是因为已知这些基序在外显子中减少（Maizels和Gray 2013）。在这里，我们报告了对蛋白质编码区域中G4进化趋势的全面基因组范围分析。基于灵长类基因组中最高的保守性，我们选择了这些物种进行研究，重点关注PGQS的位置和精确序列保守性分析。尽管基于G4基序序列构成的计算机预测PGQS数量可能会高估基因组中PGQS的数量（Zhang等人2023a, b b），但这些算法提供了表征潜在G4结构的宝贵工具（Hon等人2017；Huppert和Balasubramanian 2007；Kikin等人2006）。此外，使用基于基序的PGQS识别可以与G4折叠能量的估计相结合（Lorenz等人2011），以便更真实地了解PGQS的折叠潜力。对PGQS出现的实验验证以及这些结果与先前计算机预测PGQS映射的相关性突显了它们的实用性（H?nsel-Hertsch等人2016, 2018）。我们的分析不仅限于绘制和比较人类特有的编码G4；所有PGQS都被识别和分析，不论其来源的灵长类基因组如何。我们发现，PGQS的数量与可用的同源基因数量相关，在一定程度上也与与人类的系统发育距离相关。PGQS基序在亲缘关系密切的物种中表现出高度共定位。从最小折叠自由能推断出的热力学稳定性成为与进化模式相关的重要因素：低稳定性的PGQS（最小折叠自由能≥–10 kcal/mol）更具保守性，而高稳定性的PGQS（最小折叠自由能≤–30 kcal/mol）的保守性较低，尽管这两类序列的保守性仍高于CDS基线。一致地，插入/缺失评分与PGQS的最小折叠自由能呈负相关，这表明稳定的基序与插入或缺失事件之间存在关联。与此一致的是，富含G的串联重复序列表现出较高的插入/缺失突变率，这与它们倾向于折叠成高度稳定的G4s相符。总体而言，这些发现表明PGQS同时作为保守元件和结构不稳定性的来源，反映了在保持序列功能的同时通过结构形成产生不稳定性的拮抗性选择压力。

方法
**同源编码序列（CDS）的获取与分析**
我们从Ensembl数据库（http://www.ensembl.org/info/data/ftp/index.html）下载了所有可用灵长类物种（n=22）的CDS fasta序列。物种科学名称、基因组数据和PGQS数量的列表见补充信息表S1。对于所有物种，原始CDS fasta文件被分割并过滤，以选择每个基因的最长完整转录本（包含起始密码子ATG、终止密码子TAA/TGA/TAG以及长度为三个核苷酸的倍数）。然后，使用Biomart工具（https://m.ensembl.org/info/data/biomart/index.html）以过滤后的转录本标识符作为输入查询，以获取跨灵长类动物的同源基因。分别为每个物种获取同源基因列表，并过滤以包括与人类基因具有80%或更高身份匹配的非人类灵长类基因。最后，对于每个在多个物种中至少有五个同源基因的基因，将CDS序列合并成一个单独的fasta文件（n=18,346）。

**参考系统发育树**
物种的系统发育拓扑结构是从使用TimeTree Version 5构建的主树中获得的（https://timetree.org/）。随后使用ETE 3 v3.1.3库（Huerta-Cepas等人2016）编辑和操作该树，该库支持系统发育树的编辑、修剪和根定。对于每组同源CDS，提取了一个仅包含该组中代表的物种的子树。如果给定基因集中没有这些物种，则应用ETE3中也提供的中点根定方法。

**多序列比对（MSA）**
同源CDS被聚合成fasta文件，并使用MACSE v2.07中实现的密码子敏感算法进行比对（Ranwez等人2011）。该程序同时生成核苷酸和氨基酸比对结果，同时保留了密码子边界，从而确保了在蛋白质水平上体现选择性约束的生物学上一致的比对结果。进化与祖先模型构建：进化推断和祖先序列重建是使用BASEML（PAML v4.10.7；Yang 2007）在HKY85核苷酸替换模型（Hasegawa等人1985年）下进行的，该模型考虑了离散伽马率异质性（四个类别；Yang 1994年）。祖先状态是通过经验贝叶斯重建方法推断出来的，为每个核苷酸位点提供了后验概率（Yang等人1995年）。为了纳入祖先的不确定性，我们通过从这些后验分布中抽取核苷酸来生成100个随机实现。

潜在G-四联体形成序列（PGQS）的识别：PGQS的检测使用了R包pqsfinder（Hon等人2017年），没有应用评分阈值，并对PGQS的膨胀和错配进行了搜索容忍度设置。MSA中的每个CDS都单独进行了PGQS的筛查。这种方法使我们能够无偏地识别出跨物种的所有PGQS，即PGQS可以存在于任何一种灵长类物种中。当PGQS的起始和结束MSA位置重叠时，这些基序被视为独特的（PGQS基序）。此外，我们应用了一个过滤步骤，只保留平均得分≥40的PGQS基序。

每个PGQS基序的5′和3′端相邻的序列被定义为上游和下游侧翼区域。每个侧翼区域的长度被设置为与相应的PGQS基序长度相匹配。然而，当PGQS基序位于CDS的非常开始或结束时，侧翼区域可能比PGQS基序本身短。

串联重复序列的搜索：采用与PGQS识别相同的方法，我们使用了一个内部脚本来寻找符合正则表达式5′-NxL1NxL1NxL1Nx-3′的串联重复区域，其中N可以是A、T、G或C，x≥3个核苷酸，L代表单核苷酸环。与PGQS类似，当它们的MSA起始和结束位置重叠时，这些基序被视为独特的（串联基序）。

最小自由能计算：使用ViennaRNA包中的RNAfold工具（Lorenz等人2011年）计算每个PGQS的最小自由能。使用了–g选项来考虑G4序列的折叠能量。PGQS基序的能量被确定为所有与比对坐标重叠的PGQS的平均能量。根据其热力学稳定性，每个PGQS基序被分为两组：（i）高稳定性PGQS基序，其最小折叠自由能值≤?30 kcal/mol；（ii）低稳定性PGQS基序，其最小折叠自由能值≥?10 kcal/mol。

对于每个PGQS基序，从MSA中选取了一个等长的随机区域，确保不与原始PGQS基序坐标重叠。随后，从该随机区域中提取与特定物种的PGQS长度相匹配的随机子序列。每个随机基序的能量被确定为组成它的所有特定物种子序列的平均能量。然后根据应用于PGQS基序的相同标准，将这些随机基序分类为高热力学稳定性和低热力学稳定性组。这些基序在我们的分析中用作结构对照。

核苷酸水平指标：使用BASEML在HKY85模型下估计每个位点的替换率，该模型考虑了离散伽马率异质性，如上所述。对于每个比对位置，BASEML提供了该位置属于每个伽马率类别的后验概率，预期替换率是按类别加权后的平均替换率。替换率不仅提供了进化变化的估计，还间接反映了保守性，较低的替换率表示更强的进化稳定性，而较高的替换率反映了更快的分化。

对于相同的MSA，插入/缺失（indel）得分被定义为现代序列在每个位置上显示间隙的比例。

然后，通过提取基序范围内的每个位点的值，并取中位数作为基序的值，计算每个PGQS基序、随机基序和串联重复序列的基序水平替换和插入/缺失得分。为了实现跨基因比较，所有基序指标都在每个CDS比对中使用z分数进行了标准化，z分数定义为：（基序指标中位数 – CDS指标平均值）/ CDS指标标准差。

氨基酸水平指标：使用MACSE编码器意识到的MSA在氨基酸水平上计算保守性、替换和插入/缺失得分。在核苷酸水平生成的100个随机祖先重建被翻译成氨基酸序列。对于每个比对位置，通过将现代序列中的每个氨基酸与这些重建中的每个氨基酸进行比较来计算保守性得分。匹配的情况被赋予值1，不匹配的情况被赋予值0。最终得分是所有现代-祖先比较的平均值。此外，替换得分被定义为保守性得分的补数（1 - 保守性得分）。较高的保守性值表示现代和重建的祖先氨基酸之间有更强的一致性，而较高的替换值表示更大的进化分化。最后，插入/缺失得分被定义为现代序列在比对的每个位置上显示间隙的比例。

然后通过提取基序范围内的每个位点的值，并取中位数作为基序的值，计算基序水平替换、保守性和插入/缺失得分。之前应用于核苷酸水平指标的相同z分数标准化程序也应用于所有基于氨基酸的指标。

所有统计分析都是使用Python版本3.11.11进行的。突变率的变化是通过应用配对Wilcoxon秩和检验并调整P值来分析的。

为了评估潜在G4序列在CDS中的全基因组发生情况，我们获取、处理并比对了来自二十二种灵长类的同源CDS序列（补充信息：表S1）。PGQS被映射到18,346个多物种比对中。尽管我们使用了基于人类的同源性，但使用先前比对的序列使我们能够在灵长类同源CDS中进行非人类偏见的PGQS搜索。不同物种之间的PGQS数量有所不同（补充信息：表S1）；人类（n=276,464）、黑猩猩（n=257,366）和猩猩（n=252,812）的PGQS发生率最高，而丛猴（n=130,101）和跗猴（n=135,853）的PGQS数量最少。总体而言，18,285个灵长类CDS至少包含一个PGQS。PGQS的数量与灵长类之间可用的同源CDS数量有很强的相关性（Spearman’s rho=0.89，p=6.89e-08），并与与人类的系统发育距离有中等程度的相关性（Spearman’s rho = ?0.65，p=2.12e-03）。一些灵长类基因组组装，特别是鼠狐猴和丛猴的基因组组装，大部分是不完整的。因此，PGQS数量的变异性可能是组装质量的结果。

我们通过估计PGQS基序共现的比例（在≥2个物种中存在的PGQS基序）来分析系统发育关系对PGQS基序出现的影响。使用以人类为中心的（9,142个PGQS基序）物种×物种比较，我们发现PGQS基序共现的比例在一定程度上反映了人类与其他物种之间的系统发育距离（Spearman’s rho = ?0.56，p=7.21e-03；图1b）。人类×黑猩猩、人类×绿猴和人类×跗猴的PGQS基序共现比例分别为90%、74%和41%。总的来说，PGQS基序共现分析结果表明，密切相关的灵长类物种之间存在高度的PGQS基序共定位。

图1：完整尺寸图像。a 分析的22种灵长类物种的系统发育表示。每种物种识别的PGQS基序总数用括号表示；b 9,142个同源基因中共享的PGQS基序的以人类为中心的成对比较。

为了评估潜在G4序列对同源基因保守性的影响，我们基于其热力学稳定性分析了DNA PGQS的保守性。较低的最小折叠自由能值表示更稳定的结构，较高的最小折叠自由能值表示较不稳定的结构。我们发现CDS替换得分与PGQS基序的最小折叠自由能呈负相关，这反过来反映了CDS保守性与PGQS基序最小折叠自由能之间的正相关（图2和3）。

图2：完整尺寸图像。a PGQS基序及其上游和下游侧翼区域的替换z分数中位数；b 按热力学稳定性类别分层的侧翼区域中的替换z分数分布；c 根据热力学稳定性类别比较PGQS基序及其相邻区域的替换z分数。

图3：完整尺寸图像。a 三个基序类别中的替换和插入/缺失动态：上面板显示串联重复序列，中间面板显示高稳定性PGQS和低稳定性PGQS，下面板显示高稳定性随机和低稳定性随机；b 三个基序类别中的插入/缺失z分数分布：上面板显示串联重复序列，中间面板显示高稳定性PGQS和低稳定性PGQS，下面板显示高稳定性随机和低稳定性随机。

在比较不同热力学稳定性类别之间的替换值分布时，我们观察到低稳定性PGQS（最小折叠自由能≥?10 kcal/mol，n=5,920）和高稳定性PGQS（最小折叠自由能≤?30 kcal/mol，n=101）之间存在重叠（图3a，中间面板）。尽管如此，低稳定性PGQS的替换分数分布倾向于更强的保守性。为了进一步将这些模式置于上下文中，我们使用全长CDS替换率作为基线来比较PGQS的z分数（见方法部分）。PGQS的替换水平根据预测的热力学稳定性显著不同（配对Wilcoxon秩和检验并调整P值，p=6.36e-15；补充信息：表S2）：低稳定性PGQS倾向于比其宿主CDS更保守（中位数=?0.44），而高稳定性PGQS的替换率更接近CDS本身（中位数=0.0002）（图3a，中间面板）。此外，紧邻PGQS基序的区域（无论是上游还是下游）通常表现出与PGQS本身相似的替换模式，无论是在总体分析中还是按预测的热力学稳定性类别分层时（图2）。

为了测试与PGQS基序长度相当的区域的保守性模式是否相同，我们检查了编码DNA中的替换分数是否与热力学稳定区域内在相关（图3a，底部面板）。与PGQS基序相比，低稳定性随机和高稳定性随机区域都比低稳定性PGQS和高稳定性PGQS本身更保守（配对Wilcoxon秩和检验并调整P值，p<0.01；补充信息：表S2）。与它们的宿主CDS相比，随机区域代表了所有组中最保守的组，表现出低于CDS基线水平的替换分数（中位数=?0.49和?0.54；配对Wilcoxon秩和检验并调整P值，p<0.01；补充信息：表S2）。重要的是，这表明观察到的PGQS基序的保守模式并不适用于长度相当的区域，也不适用于在同一热力学稳定性范围内的区域。因此，尽管G4基序是保守的，但它们的保守性仍然低于同等大小和热力学稳定性的随机基序。

为了评估观察到的PGQS基序的保守性模式是否仅仅归因于它们的重复组成，我们将它们与长度和碱基含量相当的短串联重复序列进行了比较（图3a，顶部面板）。具体来说，我们分析了5′-NxL1NxL1NxL1Nx-3′基序的保守性z分数，其中N可以是A、T、G或C，x≥3个核苷酸，L代表单核苷酸环。我们发现C< G < A < T的DNA保守性较低（中位数分别为?0.30、0.008和0.02、0.75）。富含A（n=442）、C（n=166）和G（n=75）的串联区域显示出与高稳定性PGQS相似的保守性水平（配对Wilcoxon秩和检验并调整P值，p>0.01；补充信息：表S2）。然而，可以观察到富含A和G的串联区域显示出与高稳定性PGQS相比更大的替换值范围，而富含C的串联区域显示出更接近PGQS基序的模式。低稳定性PGQS区域在与串联重复序列比较时仍然具有最高的保守性水平（配对Wilcoxon秩和检验并调整P值，p<0.01；补充信息：表S2）。这些结果表明，稳定性较低的PGQSs往往比串联重复序列更保守，而稳定性较高的PGQSs虽然仍然比CDS基线更保守，但其替换模式与串联重复序列相似。使用之前描述的基于核苷酸的保守性分析方法，我们评估了PGQS基序中的氨基酸保守性。总体而言，在氨基酸水平上，组间差异小于在核苷酸水平上观察到的差异（图4a和5）。我们发现，随着PGQS热稳定性的增加，相关氨基酸的保守性下降（图4a）。PGQS相关氨基酸表现出依赖于热稳定性的模式：稳定性较低的PGQSs显示出较低的替换率（即更高的保守性），而稳定性较高的PGQSs显示出较高的替换率，表明其保守性较低（配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p = 6.60e-15；图5a，中间面板；补充信息：表S3）。相比之下，无论热稳定性如何，随机基序（无论是低稳定性还是高稳定性）都始终更为保守（配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p < 0.01；图5a，底部面板；补充信息：表S3）。与串联重复序列相比，稳定性较高的PGQSs与T-、C-或G-富集的串联重复序列没有显著差异（配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p > 0.01；图5a，顶部面板；补充信息：表S3）。

图4：氨基酸PGQS基序中的替换和插入/缺失动态。a 氨基酸PGQS基序的中位替换z分数；b 氨基酸PGQS基序的中位插入/缺失z分数。

图5：氨基酸PGQS基序与对照基序之间的比较替换和插入/缺失动态。a 三种基序类别中的氨基酸替换z分数分布：顶部面板显示串联重复序列，中间面板显示高稳定性PGQS和低稳定性PGQS，底部面板显示高稳定性随机和低稳定性随机；b 三种基序类别中的氨基酸插入/缺失z分数分布：顶部面板显示串联重复序列，中间面板显示高稳定性PGQS和低稳定性PGQS，底部面板显示高稳定性随机和低稳定性随机。

我们的结果表明，灵长类动物中的序列同源性在核苷酸和氨基酸水平上都存在趋势。一般来说，PGQSs的替换率等于或低于参考CDS基线，但与相同长度的随机基序相比，它们的替换率更高。此外，保守性受热稳定性的调节：稳定性较低的PGQSs表现为更保守的元素，而稳定性较高的PGQSs与较低的保守性相关，并且其替换模式类似于某些串联重复序列的不稳定性。因此，对低稳定性PGQSs和高稳定性PGQSs、串联重复序列以及随机基序的比较分析表明，PGQSs具有更大的功能相关性，因为它们比随机DNA区域更保守，但进化速度更快，这表明它们可能作为不稳定性位点或快速基因组适应的位点。

我们的结果揭示了一个趋势，即高稳定性PGQSs与较低的保守性和较高的替换率相关，而低稳定性PGQSs则更保守。这种依赖于热稳定性的模式在核苷酸和氨基酸水平上都是一致的。为了阐明与高稳定性PGQSs相关的突变模式，我们计算了PGQSs、串联重复序列和随机基序的DNA插入/缺失分数。插入/缺失z分数在不同基序中有所不同；较高的值表示富含间隙的序列，即插入或删除更受青睐。我们发现，CDS插入/缺失分数与PGQS基序的最小折叠自由能呈负相关，且紧邻PGQS基序的区域（上游和下游）通常表现出与PGQSs本身相似的插入/缺失模式，无论是在总体分析中还是在按预测的热稳定性类别分层后（低稳定性PGQSs和高稳定性PGQSs）（图6）。高稳定性PGQSs的插入/缺失率高于低稳定性PGQSs（分别为2.14和-0.20；配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p = 9.93e-17；图3b，中间面板；补充信息：表S4）。对于串联重复序列，我们观察到插入/缺失分数随A < T < C < G的顺序逐渐增加（图3b，顶部面板）。对于随机基序，无论热稳定性如何，插入/缺失分数都保持负值，反映了相对于全基因组模式的较低插入/缺失率（高稳定性随机中位数 = -0.41，低稳定性随机中位数 = -0.25；图3b，底部面板）。高稳定性PGQSs和G-富集的串联区域显示出所有分析组中最高的中间值（分别为2.14和1.84；图3b，中间和顶部面板）。一小部分映射到CDSs上的PGQS基序（13%）呈现独特的G4基序，仅存在于一个物种中。有趣的是，当考虑PGQS的热稳定性时，高稳定性PGQSs的物种特异性显著高于低稳定性PGQSs（分别为54%和13%；Fisher精确检验p = 1.27e-23，比值比 = 7.9）。这些结果表明，高度稳定的PGQSs与物种间的局部插入或删除事件相关，从而促进了物种特异性G4基序的出现，以及在其他物种中同源位点缺乏G4。

图6：核苷酸水平上PGQS基序及其侧翼区域的比较插入/缺失谱型。a PGQS基序及其上游和下游侧翼区域的中位插入/缺失z分数；b 按热稳定性类别分层的侧翼区域中的插入/缺失z分数分布；c 根据热稳定性类别比较PGQS基序及其相邻区域的插入/缺失z分数。

由于基于氨基酸的对齐通常比基于相应核苷酸的对齐更准确（Abascal等人，2010年），我们还在氨基酸水平上评估了PGQSs、串联重复序列和随机基序的突变模式（图4b和5b）。高稳定性PGQSs继续显示出明显高于低稳定性PGQSs的插入/缺失率（分别为1.92和-0.09；配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p = 2.19e-11；图5b，中间面板；补充信息：表S5）。与G-富集的串联基序（中位数 = 0.90）一起，高稳定性PGQSs在所有基序中显示出最高的插入/缺失率，它们之间没有明显差异（配对Wilcoxon秩和检验，考虑多重比较后的p值调整，p = 1；图5b，顶部面板；补充信息：表S5）。对于串联重复序列，氨基酸插入/缺失分数从A-富集到G-富集的基序逐渐增加，这与DNA插入/缺失分数的模式一致（图5b，顶部面板）。对于随机基序，氨基酸插入/缺失z分数保持负值（图5b，底部面板）。有趣的是，在氨基酸水平上，只有高稳定性PGQSs和G-富集的串联重复序列显示出z分数值的下降。所有其他基序的z分数值从DNA到氨基酸都有轻微增加。对于T-富集的串联基序，观察到相反的趋势，中位插入/缺失z分数从DNA水平的-0.12变为氨基酸水平的0.20，导致相对于氨基酸CDS基线的插入/缺失率更高。

这些结果表明，高稳定性PGQSs和G-富集的串联重复序列与编码区域内的较高插入/缺失突变率相关，这一模式在核苷酸和氨基酸水平上都是一致的，其中G-富集的串联重复序列的效应最为显著，可能反映了G-富集基序潜在地折叠成非常稳定的G4结构。

在这项研究中，我们绘制并表征了所有可用灵长类动物CDSs中的PGQS基序。由于在比较G4的出现和位置时通常使用偏向人类的方法，我们实施了一种结合了多个物种同源基因对齐的方法，然后搜索PGQSs，以减少这种偏见。为了对G4进行适当的进化分析，我们专注于灵长类动物。有趣的是，我们能够高效地绘制出共享的和物种特异性的PGQSs。需要注意的是，基因组组装的质量影响了最终识别的潜在G4的数量，特别是在小鼠狐猴和丛猴中。因此，PGQSs在与人科相关的物种中得到了更好的表征。稳定性估计是使用在37°C下定义的热力学参数计算的，这与灵长类动物作为温血动物的生理特征一致。然而，G4结构的稳定性和活性会随温度变化，因为最近的研究表明RNA G4可能作为热传感器调节哺乳动物的基因表达（Zhang等人，2025年）。总的来说，这些观察结果强调了需要进一步研究以更好地表征G4编码景观，不仅在非人类灵长类动物中，还包括具有不同体温调节策略的物种。

G4在DNA复制和转录过程中引起的诱变效应已被描述（Kruisselbrink等人，2008年；Lemmens等人，2015年；Lopez等人，2017年；Wang和Vasquez，2017年；Yadav等人，2014年），这可能解释了它们在编码区域中的代表性不足。G4在不同物种中的出现已被广泛研究（Marsico等人，2019年），并且主要使用以人类为中心的保守性分析（Frees等人，2014年）。我们通过使用一种平行的、较少偏见的多物种方法在灵长类动物同源基因中识别和分析G4编码保守性来解决这个问题。我们的结果表明，PGQSs受到复杂且有时是相互矛盾的选择压力，因为它们同时构成了序列保守的区域和结构不稳定的焦点。实际上，例如RNA二级结构可以根据上下文被正面或负面选择（Gebert等人，2019年；Katz和Burge，2003年；Shabalina等人，2013年），我们的结果显示G4的热稳定性影响其在编码区域中的普遍性。这种模式表明，PGQS序列的保守性是为了功能原因而维持的，而这些序列形成的结构本身可能代表基因组不稳定的来源（Bochman等人，2012年；Mohanty等人，2025年）。

PGQS基序与系统发育接近性之间的关联可能反映了PGQS在CDSs中的进化更替（Frees等人，2014年；Mohanty等人，2025年）。PGQS基序可以通过随机突变过程出现，这些过程生成或扩展能够支持G4折叠的鸟嘌呤片段（Gong等人，2021年）。随着时间的推移，这些基序可能会根据它们的功能影响而被不同程度地保留或丢失（Guiblet等人，2021年）。因此，PGQS基序可能在更不同的谱系中逐渐衰减或消失，导致在密切相关的物种中基序的共定位程度更高。这种模式因此与PGQSs的反复出现及其随后通过突变过程和选择约束的去除之间的动态平衡一致（Frees等人，2014年；Giblet等人，2021年；Mohanty等人，2025年）。

G4的保守趋势，通过我们结果中观察到的负值以及接近零的z分数得到了反映，表明了它们的功能生物学重要性。这种模式在哺乳动物中已被描述，并且在非哺乳动物生物中不太明显（Frees等人，2014年）。由于鸟嘌呤片段的中断使得G4的形成不可行，因此在这种背景下突变不太被容忍，从而维持了这些序列的结构能力（Kim，2019年；Nakken等人，2009年）。这种保守性可以通过它们在基因组的特定位置以及作为多种细胞过程中的调控元素的角色来解释，这些过程对维持这些序列在特定位置施加了选择压力（Rhodes和Lipps，2015年；Shen等人，2021年；Varshney等人，2020年）。实际上，G4在DNA和RNA中都扮演着重要的调控角色（Varshney等人，2020年）。在DNA中，它们调节复制、转录、基因组稳定性、端粒生物学和染色质组织，并参与基因表达的精细调控（Bhattacharyya等人，2016年；Varshney等人，2020年）。在RNA中，G4主要参与翻译的调节、可变剪接、microRNA成熟和细胞应激反应，作为影响RNA命运和功能效率的动态结构元素（Bhattacharyya等人，2016年；Fay等人，2017年；Varshney等人，2020年）。此外，RNA G4结构比DNA G4更容易形成，并且通常表现出更高的热稳定性，这加强了它们作为转录后生物学中调控模块的重要性（Bhattacharyya等人，2016年；Cheong和Moore，1992年；Kumari等人，2007年；Sacca等人，2005年）。有趣的是，尽管箱形图显示低稳定性PGQS和高稳定性PGQS的替换z分数在负值和接近零的值上有重叠，但低稳定性PGQS的一个范围内集中了较低的替换值，表明这一亚组具有更高的相对保守性。高度保守的G4序列虽然热力学稳定性中等或较低，但由于它们在调控具有相似表达模式的基因翻译中的潜在作用，因此得以保留（Endoh和Sugimoto 2016）。实际上，稳定性较低的G4基序往往更具保守性，这反映了正向的选择压力，这种压力减少了更稳定G4序列的负面影响，并保留了其有益的调控功能（Lombardi等人2019）。具有不稳定特征的G4序列，如较大的环结构、凸起部分和错配碱基，可能表现出更强的折叠动态性、结构灵活性和构象变异性（Meier等人2018；Tippana等人2014；Varizhuk等人2017）。这种能力使这些G4序列能够作为响应细胞环境的分子开关，从而赋予其动态和可重构的适应性（Dong等人2022；Zhang等人2023a, b a）。因此，一种能够根据细胞信号形成和分解的短暂性PGQS（PGQS），可能比永久固定的结构具有更敏感的调控作用（Gilbert和Marenduzzo 2025；Robinson等人2021）。尽管更稳定G4序列的保守性下降趋势可能表明其具有总体上的负面影响，但这些序列具有更高的物种特异性，并且其indel z分数高于稳定性较低的G4序列，这表明它们可能仅对特定物种产生有限的影响，而不会损害蛋白质的主要功能。为了进一步了解与G4相关的突变模式，我们分析了潜在非B DNA重复序列之间的indel（插入/缺失）发生率。有趣的是，G4基序与DNA水平上的indel突变增加有关，这可能导致氨基酸序列的改变。先前的研究已经报告称，人类体内的潜在G4序列中SNP（单核苷酸多态性）较少（Nakken等人2009）。最近，Du等人表明G4结构确实与较高的突变潜力相关，尽管该研究中并未区分SNP和indel的具体影响（Du等人2014）。在我们的结果中，高稳定性的PGQS以及富含G的串联重复序列在形成G4结构时表现出较高的indel z分数，这表明它们相对于识别它们的CDS（编码序列）具有更高的突变阈值。这些结果表明G4稳定性与indel发生之间存在关联，尽管因果关系的方向无法直接推断。这种模式可能反映了G4诱导的突变、indel驱动的G4基序出现，或共同的序列特性。大多数研究集中在第一种情况，即G4结构促进了插入和删除事件，因为它们对细胞复制机制构成了拓扑挑战，并且还调节基因组的稳定性、对损伤的敏感性以及修复途径的效率（Pavlova等人2021）。G4结构可能会阻碍复制叉的进展，这种现象称为复制叉停滞（Batra等人2025；Paeschke和Burkovics 2020）。确实，更稳定的G4结构往往与复制失败的风险增加相关（Bochman等人2012；Lombardi等人2019；Piazza等人2015；Williams等人2023），并且更难以被负责解开这些结构并确保复制顺利进行的专门解旋酶所识别（Budhathoki等人2014）。复制叉停滞和G4解旋效率低下带来的压力可能导致DNA双链断裂的积累（Sato等人2021；Varshney等人2020）。这些损伤的修复主要通过末端连接途径进行（Lemmens等人2015），但这些途径本身容易出错，常常导致断裂点插入和删除（Cisneros-Aguirre等人2022）。总体而言，无论这种关联背后的机制是什么，观察到的模式都与进化过程一致：高稳定性的PGQS在编码区域可能因有害效应而被纯化选择所清除，而稳定性较低的PGQS则因其调控作用大于结构成本而被容忍并可能持续存在（Endoh和Sugimoto 2016；Lombardi等人2019）。包括G4在内的DNA串联重复序列被认为是致突变的，并且已被证实与近30种人类遗传疾病有关，其中许多疾病主要通过改变蛋白质功能影响神经系统（Mirkin 2007；Paulson 2018；Usdin 2008）。这些序列之所以被认为是致突变的，是因为它们的重复结构有利于改变重复单元数量的机制。提出的机制包括减数分裂过程中的不等交叉、逆转录转座事件以及DNA复制过程中的链滑移，后者被认为是导致重复序列扩展或收缩的主要机制（Fan和Chu 2007；Tanudisastro等人2024）。我们认为，G4序列的低氨基酸变异性和高indel发生率可能反映了框架内基序的扩展事件（Castel等人2010；Iyer等人2015；Renton等人2011）。因此，稳定的G4序列可能是通过单个或密切相关的灵长类动物中富含G的串联重复序列的扩展而产生的。另一种解释G4编码基因组内容的机制是，祖先G4基序仅通过功能选择在特定物种中得以保留。虽然G4诱导的蛋白质功能改变可能对某些物种产生有害影响（如人类中的情况Conlon等人2016；Fratta等人2012），但对其他灵长类动物可能具有中性或有益的影响，这反映了遗传多样化的机制。这种潜在的G4诱导的变异性的功能后果需要进一步验证，并可能有助于理解结构易变序列对蛋白质编码基因组区域进化的影响。总之，我们的研究结果将G4序列定位为进化过程中的核心组成部分。一方面，它们对于需要序列保守性的调控功能至关重要；另一方面，它们的结构构象可能对基因组稳定性构成挑战。