雅鲁藏布江岩蜥（Laudakia）沿海拔–经度复合梯度的连续基因组分化与谷首种群基因组脆弱性

《Biology》：Genomic Differentiation Along an Altitudinal–Longitudinal Gradient in Yarlung Zangbo Rock Agamas (Laudakia)

　　

极端海拔梯度是检验环境梯度如何塑造基因组变异的自然实验，但在高海拔爬行动物中仍缺乏表征。研究人员分析了来自青藏高原雅鲁藏布江谷五个种群、跨越2311–3758 m的44只岩蜥（Laudakia sacra和L. wui）的全基因组重测序数据（~17×）。L. sacra内变异形成连续海拔渐变而非四个离散单元，第一主成分（PC1）解释41.29%的方差；多样性随海拔单调下降（π下降约9倍），连锁不平衡（LD）半衰减延长27倍，近交系数（F_ROH）达到0.155，有效种群大小（N_e）在谷首最低。结合置换的冗余分析（RDA）识别出593个海拔相关单核苷酸多态性（SNP），并得到基于秩的渐变筛选（相对独立性的6.1倍富集）共同支持；选择清除分析在十个两两比较中检测到4573个候选基因，极端种间对比（C1 vs. C5）涵盖细胞骨架、表观遗传、糖酵解和氧感知功能。低多样性、高近交、小有效规模和较大预测遗传偏移在谷首重合，使其成为遗传上最贫乏和气候最脆弱的种群。这些结果揭示沿海拔–经度梯度的基因组分化，符合推定的初期局域适应，并确定高海拔爬行动物的保护优先事项。

研究背景与问题：青藏高原是生物多样性热点，雅鲁藏布江自西向东形成从高山荒漠到亚热带峡谷的连续生态梯度，为研究海拔适应分化和物种形成提供天然实验室。Laudakia Gray, 1845（鬣蜥科）约20种，中国已知6种，其中L. sacra、L. wui、L. tuberculata和L. papenfussi为西藏特有；L. sacra主要分布于雅鲁藏布江中上游较高海拔干旱河谷，L. wui分布于下游帕隆藏布和雅鲁藏布峡谷热带亚热带季风雨林，约700–2300 m，二者分布、海拔生态位和表型差异明显。Yan et al.组装了L. sacra首个参考基因组（1.80 Gb，N50=20.19 Mb），鉴定HIF1A、TIE2、ERCC6和MSH2等高海拔适应基因，并检测到L. sacra与Thermophis baileyi在MLIP基因上的趋同信号；Wei et al.和Tan et al.先后发表L. wui染色体级基因组。但Laudakia种群基因组研究仍落后，其种群结构、高海拔适应机制和物种形成过程在全基因组尺度尚未解决。本研究以全基因组重测序检验两个假设：基因组变异沿谷海拔–经度梯度呈连续体而非离散单元；沿梯度分化相对种界较浅，种内海拔分化弱于L. sacra与L. wui的种级分化。预测高海拔种群遗传多样性降低、基因组近交升高。该研究为高海拔爬行动物保护优先排序提供基线。
关键技术方法：研究人员于2025年9月沿雅鲁藏布江自西向东采集5个种群44只岩蜥（C1–C4为L. sacra，C5为L. wui；C1 Nyemo 3758 m，n=9；C2 Zetang Bridge 3582 m，n=9；C3 Gyaca 3265 m，n=10；C4 Mainling 3017 m，n=8；C5 Layue Hot Spring 2311 m，n=8），肌肉组织保存于95%乙醇，凭证标本存西藏自治区高原生物研究所。采用全基因组重测序（whole-genome resequencing, WGR）（~17×），比对L. sacra参考基因组，GATK调用SNP和InDel，ANNOVAR注释。以IQ-TREE、主成分分析（principal component analysis, PCA）、ADMIXTURE和F_ST分析种群结构；VCFtools计算π、H_o、H_e、F_IS和Tajima’s D，PopLDdecay分析LD，Mantel/偏Mantel检验分化与海拔、地理距离关系；TreeMix推断基因流；F_ST与π比选择清除及GO/KEGG富集；RDA与Spearman基因型–环境关联分析并置换种群海拔标签；候选基因区间渐变分析；d_XY、固定差异和Hudson F_ST量化种间分化；RDA投影SSP585估计遗传偏移并计算非适应风险（risk of non-adaptedness, RONA）；SMC++（基于顺序马尔可夫合并的种群历史推断方法）和LD法重建N_e；PLINK检测纯合片段（runs of homozygosity, ROH）并计算F_ROH；ANNOVAR分类SNP计算非同义/同义负荷比。
主要研究结果：3.1 测序数据与变异检测：获得1357.18 Gb clean data（平均30.85 Gb，Q30>96.0%），平均比对率99.87%、深度~17×、覆盖98.7%；保留21,911,886个SNP和2,651,703个InDel，其中外显子SNP 224,115个（非同义89,100个）；LD修剪后594,884个SNP用于结构分析。3.2 种群遗传结构与物种分化：IQ-TREE树支持C1、C2各自单系，C5形成长枝；PCA的PC1（41.29%）呈C1→C2→C3→C4→C5连续梯度，PC2（12.82%）分C3，PC3（7.37%）分C4与C5；ADMIXTURE最优K=5，K=2分高海拔C1/C2与中低海拔C3/C4/C5，K=5时各种群近纯合（Q>0.99）。C1–C5 F_ST最高0.898，C5相关比较均>0.80，种内为0.374–0.675；种间d_XY为0.646、种内0.136（4.8倍），固定差异35.0% vs. 约2%，Hudson F_ST 0.754 vs. 0.398；海拔轴种内0.142 vs. 种间0.639（4.5倍）。Mantel检验中海拔差异r=0.943（p=0.0166），地理距离r=0.771（p=0.104）；偏Mantel控制地理后海拔仍显著，控制海拔后地理不显著，提示环境梯度成分而非纯距离隔离，但海拔为复合代理。3.3 遗传多样性与演化动态：L. sacra内π从C4 0.0013降至C1 0.00015（约9倍；Spearman ρ=?1.00，线性r=?0.99，p<0.05，R²=0.97）；LD半衰减C3为2.4 kb、C4为2.9 kb、C1为65.3 kb、C2为42.4 kb，C5为11.5 kb。3.4 基因流与种群演化关系：TreeMix拓扑为（C2，（（C4，（C5，C3）），C1）），一条权重0.444的迁移边连接C1谱系与C3祖先；方向不可识别，读作历史基因流，近期基因流很有限。3.5 选择清除与候选基因：十个两两比较共4573个候选基因命中，C3 vs. C4最多（877个），C1 vs. C5为112个，C1 vs. C2为100个；C1 vs. C5的GO富集仅细胞骨架相关（肌动蛋白丝组织、肌动蛋白丝结合、钙调蛋白结合），其他比较富集嗅觉受体、趋化因子和微丝马达、伤口愈合和单糖转运；KEGG仅嗅觉转导和尼古丁成瘾显著（C1 vs. C3）。C1 vs. C5候选包括SHROOM4、FLNC、HDAC8、DGKD、BMP15、NR2E1和RPS4，FLNC和SHROOM4在C1 vs. C5与C2 vs. C5重复，提示细胞骨架重塑作用；因C1与C5海拔和物种均不同，不能仅归因于海拔。3.6 海拔相关基因组变异：四个L. sacra种群平均等位基因频率从C4 0.177降至C1 0.036；|ρ|≥0.9初筛18,402个（9.3%），其中3791个（1.9%）完全单调；RDA最大载荷（>99.5百分位）1049个与|ρ|≥0.9重叠593个，富集6.1倍；RDA海拔解释39.4%基因型变异，与PC1 41.29%接近，种群海拔标签置换精确p=1/24=0.042；两方法高度一致（r=0.958；符号一致99.1%，RDA选中SNP为100%），593个为稳健海拔集合。完全单调SNP平均跳跃指数（jump index）0.60，线性零模型为0.48（Kolmogorov–Smirnov p<10^?300），但最大增量分布于三个区间（32.5%、37.8%、29.6%），未识别单一海拔阈值。3.7 候选基因锚定海拔渐变：13个候选基因平均|Spearman ρ|为0.70（0.51–0.76），全基因组为0.29，全部位于上四分位（Mann–Whitney p<0.001）；SHROOM4、HDAC8、BMP15、DGKD在scaffold 26约115 kb区间聚集，与可能共适应区域一致但未证明；Eno3、PKM、HK2、CYBA、Piezo1、ARRB2、MTHFD2、SLC4A5也有强渐变（|mean ρ|=0.51–0.76）。3.8 气候变化下的基因组脆弱性：SSP585（2071–2100）下遗传偏移从C4 1651 RDA单位升至C1 2102，C1同时具最低多样性、最高近交和最大预测偏移；RONA同序（C1>C2>C3>C4，Spearman r=+1.000），但二者同源投影，属一致性检查而非独立验证。3.9 人口历史与基因组近交：SMC++显示C1、C2持续小有效规模（C1 N_e≈164；C2≈60），C3保持较大（原文~104），C4轨迹最不一致，C5古代扩张后中全新世收缩；LD法当代N_e排序为C1ROH随海拔升高：C1最高均值0.155（范围0.072–0.272；约278 Mb ROH），C2为0.031，C3/C4为0.006，C5为0.001；C1含最多0.5–1 Mb及最长2–4 Mb ROH，提示近期近交和小近期N_e。F_IS在C1为?0.181（95% CI ?0.246至?0.115，9个体均负），C2和C5亦负，C3/C4近零；该符号与近交预期相反，未过滤重算变化≤0.015，故近交陈述基于F_ROH。非同义/同义次要等位基因负荷比从C1 1.34降至C4 1.15，与C1小种群纯化选择效率降低一致。
讨论部分总结：讨论认为L. sacra四个种群沿连续海拔–经度梯度组织而非四个离散单元，PC1排序、K=2分高海拔与中低海拔；因L. wui仅一个种群，物种轴与种群轴不可分，不能声称海拔优先于物种身份；海拔与经度强相关，偏Mantel支持环境成分，但六个种群对无法识别具体气候变量。多样性梯度、LD半衰减和F_ROH支持高海拔小有效规模；C3/C4高多样性符合低峡谷冰川避难所。Serén et al.在Tenerife蜥蜴Gallotia galloti中发现阶梯式等位基因频率变化，本研究位点水平跳跃指数提示渐变较线性更陡且集中，但未识别阈值，且仅对线性零模型比较，不能证明适应性。L. wui种地位由F_ST、长枝和ADMIXTURE独立成分支持，但F_ST对比最弱且单种群限制；拓扑提示L. wui可能经边域物种形成（peripatric speciation）源自与东部L. sacra共享的中低海拔祖先，TreeMix提示C1与C3历史基因流但方向不可确定。候选基因涵盖细胞骨架重塑、表观遗传、脂/能量信号、糖酵解和氧感知，scaffold 26约115 kb聚集需单倍型和LD证据验证。保护上C1多样性最低、N_e最小、F_ROH=0.155、预测遗传偏移最大且非同义负荷升高，应优先保护；建议各种群作为独立管理单元并优先C1和未采样低海拔L. wui种群。局限包括单参考比对、SMC++近期和最深epoch弱约束、偏Mantel不能量化相对贡献、VQSR样本较小、F_IS符号异常、Tajima’s D对频率过滤敏感、N_e绝对值依赖突变率，以及Pailong–Mêdog未采样。
研究结论翻译：深切割的雅鲁藏布江谷作为一维地理模板，使种群水平基因组分化呈连续排序：在L. sacra内，海拔–经度轴是结构主导轴，多样性、LD、有效规模和近交梯度均沿其变化，符合沿梯度推定的初期分化，其中种内海拔分化（d_XY=0.142）仍远低于物种轴（0.639）。L. sacra和L. wui代表两个独立演化物种，L. wui可能经边域物种形成源自与L. sacra共享的中低海拔祖先。局域适应是否已进展到降低基因流的程度，本研究未检验，需同园或交互移植数据；高海拔适应候选基因（FLNC、SHROOM4、HDAC8、Eno3、PKM、HK2、CYBA、Piezo1）亦需功能验证。需在其他深切河谷开展类似样带，以检验该巧合是否普遍。