文章信息
摘要
中国地下鼢鼠的系统发育和物种形成尚不清楚，因为之前对形态和有限分子标记的研究产生了相互矛盾的结果。这项研究基于染色体水平的从头组装和23个种群的全基因组测序，揭示了中国8个鼢鼠物种的复杂进化史。我们发现核系统发育和线粒体系统发育以及不同的合并系统发育之间存在广泛的系统发育不一致，这可以通过基因渗入和不完全谱系排序（ILS）来解释。最近的青藏高原隆起（约360万年前；Mya）驱使凸颅鼢鼠属分化为进化枝A和进化枝B（〜3.22Mya），并且该节点中不一致的系统发育主要归因于基因渗入而不是ILS。由于地理隔离和冰川作用，进化枝A迅速分化为三个谱系，而冰川作用和C4植物扩Zhang促进了进化枝B的物种形成。ILS导致了两个快速辐射节点的不一致，而不是渗入。所有凸颅鼢鼠属物种的有效种群规模（Ne）均受到三次冰川作用的影响。古老的多态性和分歧搭便车促成了所有物种对的基因组岛。确定了与特定居住适应相关的积极选择的基因，例如心脏发育、神经发生、DNA修复和免疫反应。气候、地质构造和C4植被影响了中国鼢鼠的适应和物种形成。
关键词:鼢鼠亚科、群体基因组学、基因渗入、不完全谱系排序、物种形成
意义
系统发育是许多生态和进化研究的基础。然而，鼢鼠的系统发育和物种形成模式存在激烈争议。这项研究从基因组角度解开了鼢鼠的系统发育和物种形成。距今360万年前青藏高原近期隆升，将凸颅鼢鼠属中的6种分为高海拔的高原鼢鼠和斯氏鼢鼠以及其余4种低海拔的种。罗氏鼢鼠和斯氏鼢鼠在秦岭-淮河线以南形成，那里是冰川作用期间提供避难所的地方。基因渗入和不完整的谱系分类导致了鼢鼠的复杂系统发育。基因组岛是由于古老的多态性和分歧搭便车而形成的。这项研究的结论是，气候、地质和构造事件影响了中国鼢鼠的系统发育和物种形成。
引言
生态物种形成的遗传基础是进化生物学的关键主题之一。种群适应新环境的方式以及密切相关物种的谱系激发了生态学家和进化生物学家对该领域的兴趣(1)。系统发育树显示的物种之间的关系并不总是明确的，因为许多研究已经表明基因渗入、不完全谱系排序(ILS)和树构建的群体统计历史的广泛影响(2,3)。基因渗入和ILS在多样化的物种中留下了相似的基因组足迹，因此很难在基因组中区分它们。因此，面对进化过程的组合，重建代表类群之间真实关系的精确系统发育树始终具有挑战性（4）。高通量测序技术和相应分析的最新进展提高了我们推断系统发育关系的能力。进化生物学家现在可以通过整合不同的数据集和方法来阐明物种之间相对精确的关系。基因组学和遗传模型（包括合并理论）的发展促进了群体统计学、物种形成和适应方面的进步。这些进展使研究人员能够研究复杂的群体历史、揭示物种形成触发因素并确定适应性进化的潜在遗传机制(5)。
鼢鼠亚科（啮齿目，鼹形鼠科）生活在地下并在地下觅食，具有地下哺乳动物典型的形态特征，如四肢和尾巴较短，耳廓和眼睛缩小，以及突出的挖洞“工具”。鼢鼠亚科一般分布在半干旱和半湿润地区长江以北地区年降水量在200～800毫米之间（图1A）。该亚科由两个属组成：凸颅鼢鼠属和鼢鼠属，主要分布于古北界，包括中国、蒙古、俄罗斯远东地区、哈萨克斯坦和朝鲜。据推测，鼢鼠的适应性进化和生态物种形成是由快速的环境变化和气候振荡驱动的(6)。它们起源于中新世的蒙古(6)，然后向南扩散到上升的青藏高原(QTP)(7)。凸颅鼢鼠属仅分布于中国(8)，包含6个种，分别为中华鼢鼠(分布于中国东北至中西部，与其他几种有重叠)、秦岭鼢鼠(分布于整个秦岭)、甘肃鼢鼠(分布于中国东北至中西部)、甘肃鼢鼠(分布于整个秦岭)整个黄土高原、秦岭以北）、罗氏鼢鼠（发现于秦岭鼢鼠以南）、高原鼢鼠（散布于青藏高原，有两个种群被黄河隔离）和斯氏鼢鼠（主要发现于高原鼢鼠所处位置以南）（图1A）。鼢鼠属，由草原鼢鼠和东北鼢鼠组成，发现于中国东北地区。秦岭淮河线（QHL）根据气温和降水划分中国南北。线路以北为温带大陆性气候，冬季有雪；线路以南为热带、亚热带季风气候，炎热湿润。QHL将所有八个物种分为三组。斯氏鼢鼠和罗氏鼢鼠分布在QHL以南，高原鼢鼠和秦岭鼢鼠分布在整个昆仑-秦岭山脉。高原鼢鼠的分布以黄河为界（图1A），分为年降水量约400至800毫米的南部地区和年降水量约200至400毫米的北部地区。草原鼢鼠和东北鼢鼠仅分布于中国东北地区（图1A）。
图1.鼢鼠的分布和系统发育。(A)中国的采样地点和相应的地理特征。半干旱和半湿润地区分别在地图上以浅黄色和蓝色区域显示。长江和黄河用蓝线描绘。昆仑山和昆仑山分别用深棕色和绿色虚线标记。确切的采样地点已标记在地图上。(B)8个鼢鼠物种的19个气候因子的主成分分析(PCA)。(C)六种啮齿类动物和六种其他哺乳动物基于直系同源编码序列的具有分歧时间的系统发育树（CI在括号中显示）。(D)马戏图显示基本统计数据：(a)参考基因组每条染色体的长度；(b)GC含量；(c)重复的密度；(d)基因密度；(e)全基因组31条染色体之间的共线性分析。
鼢鼠的系统发育仍然存在争议。基于使用线粒体DNA数据的系统发育分析，Su等人。表明甘肃鼢鼠是一个基础进化枝，而秦岭鼢鼠、斯氏鼢鼠和高原鼢鼠则聚集在另一个进化枝中(8)。Li和Chen认为凸颅鼢鼠属起源于秦岭山脉，其祖先为中华鼢鼠。生活在低海拔亚热带地区的一个分支进化为罗氏鼢鼠，另一个分支迁移到寒冷和高海拔地区，后来形成斯氏鼢鼠(9)。Li和Wang认为，凸颅鼢鼠属起源于秦岭山脉，以中华鼢鼠为祖先，先分裂为秦岭鼢鼠，向北传播后分裂为甘肃鸸鹋。斯氏鼢鼠与罗氏鼢鼠(10)关系密切，向西扩散并与高原鼢鼠形成物种。Zhang等人。声称甘肃鼢鼠可以与凸颅鼢鼠属亚属的物种区分开来，并且斯氏鼢鼠与罗氏鼢鼠密切相关，这表明斯氏鼢鼠可能是从罗氏鼢鼠中分离出来的(11)。这些先前的研究是基于过时的形态学或有限数量的可用基因。然而，基因树和物种树是不同的（12）。一些进化历史，如趋同进化，可能会导致进化或生态生物学家仅根据形态学得出错误的分类法(13)。基于一个或多个基因(14)或仅线粒体基因组(15)的系统发育树会由于ILS、基因渗入、杂交(16)、水平基因转移、基因复制和基因删除(17)产生误导性结果，这将导致使系统发育推断更加复杂并导致基因树不一致（16）。
这项研究对八个鼢鼠物种分布范围内的代表性种群的基因组进行了测序，以阐明它们之间的关系。我们基于种群规模的基因组数据和联合模型探索了八个鼢鼠物种的分化顺序、物种形成模式、基因渗入和有效种群规模(Ne)波动。每个物种的种间分化时间尺度和Ne与气候和地质变化有关。在12个物种对中检测到分化程度较高的基因组岛，并为每个物种鉴定了正选择下的基因。基因组结果强调了中国鼢鼠的系统发育和物种形成。高通量测序技术和群体遗传分析方法为鼢鼠物种生成了比以前仅基于形态表型的系统发育更准确的系统发育。
结果
从头基因组组装、支架和注释。我们从8个鼢鼠物种中收集了90个个体，每个个体有2到5个种群（图1A和SI附录、图S17和表S1）。鼢鼠(中华鼢鼠)的染色体级基因组是通过基于99.2-Gb(∼33×)HiFi读数的从头组装生成的，并用125-Gb(∼42×)Illumina配对末端(PE)读数进行抛光（SI附录），表S3）。组装的基因组大小为2.70Gb（表1），与预测的大小大致相同（SI附录，图S1）。重叠群N50为59.2Mb，GC含量为41.5%（表1和SI附录，图S2和S3以及表S4）。我们使用Hi-Cdata（SI附录，表S5和S6）将91.6％的重叠群锚定到31条染色体（2n=62）（图1D，表1和SI附录，图S4和表S7）。中华鼢鼠的核型数与之前报道的核型一致(18)。
基因组组装覆盖了基准通用单拷贝直向同源物（BUSCO）数据库中4,104个总直向同源物中的3,967个（∼97.0%），表明基因组的完整性（SI附录，图S5和S6以及表S8）。PE读数对装配的映射率为99.9%，表明该装配的高精度（SI附录，表S9）。整个基因组中预测有21,656个蛋白质编码基因、0.08%的简单序列重复(SSR)和42.6%的转座元件(TE)（SI附录、图S7以及表S10和S11）。所有基因组组装和注释基因特征均与同时报道的高原鼢鼠(19)相似。
表1.中华鼢鼠基因组组装统计数据
鼢鼠亚科与竹鼠亚科的亲缘关系较鼢鼠亚科更为密切。基于不同标记的系统发育树显示鼢鼠亚科、鼢鼠亚科和竹鼠亚科之间存在冲突的系统发育（SI附录，图S14）。基于13个蛋白质编码线粒体基因，重建了串联最大似然(ML)树，以阐明这三个亚科的系统发育，揭示鼢鼠亚科与竹鼠亚科的关系更为密切（SI附录，图S15）。基于核基因组的合并拓扑有40.7%的窗口支持鼢鼠和竹鼠之间更密切的关系。其余两个不一致的拓扑由30.0%和29.3%的窗口支持（SI附录，图S16），归因于基因渗入（SI附录，表S12）。盲鼹鼠最初在大约2910万年前（Mya）与鼢鼠分化，而鼢鼠和竹鼠的分化发生在26.9Mya左右（图1C）。
鼢鼠物种之间的系统发育关系。19个气候因子（SI附录，表S1）的主成分分析（PCA）可以区分这些物种，除了中华鼢鼠与甘肃鼢鼠、东北鼢鼠和草原鼢鼠重叠。（图1B）。温度季节性和年降水量被确定为促进物种形成的两个主要气候因素（SI附录，图S18）。检测到19个气候因素中的11个与相应的遗传距离之间存在显着正相关（Mantel检验，R2=~0.215至0.496，**P<0.01；SI附录，表S1和S2）。
为了研究中国鼢鼠物种的系统发育关系，将PE读数与中华鼢鼠基因组进行比对（SI附录，表S9），并从8个物种中获得约736万个高质量单核苷酸多态性(SNP)（SI附录，表S9）。图S19）。在去除密切相关的个体后，90个样本中的51个保留下来用于后续分析。PCA（图2A）、系统发育树（图2B）、结构分析（图2C和SI附录，图S21）和遗传网络（SI附录，图S20）显示了八个谱系的独立遗传簇。由于其中的种群分化，高原鼢鼠被视为两个谱系，高原鼢鼠1（黄河以北）和高原鼢鼠2（黄河以南），用于下游分析（图2C）。所有物种均聚类为两个进化枝：高原鼢鼠1、高原鼢鼠2和斯氏鼢鼠聚类为进化枝A，中华鼢鼠、秦岭鼢鼠、甘肃鼢鼠和罗氏鼢鼠聚类为进化枝B。在K=5和6之间聚类甘肃鼢鼠和罗氏鼢鼠中观察到，并且包括罗氏鼢鼠聚类。高原鼢鼠1、高原鼢鼠2和斯氏鼢鼠在K=5和6之间（图2C）。观察到东北鼢鼠和草原鼢鼠分离，直到K=9，表明它们关系极其密切（图2C和SI附录，图S21）。
鼢鼠物种之间的遗传分化和多样性。通过FST测量成对遗传分化。高原鼢鼠1与所有其他谱系的差异最大（图2D）。高原鼢鼠的种内分化（FST=0.314；图2D）甚至大于许多其他物种的种间分化（SI附录，图S22）。中华鼢鼠是遗传多样性最低的群体，与其他谱系表现出高水平的遗传分化。高原鼢鼠1和高原鼢鼠2显示出最高的核苷酸多样性（以π测量），分别为4.1×10-3和3.9×10-3，其次是罗氏鼢鼠、斯氏鼢鼠、秦岭鼢鼠、甘肃鼢鼠、和中华鼢鼠（图2D）。除中华鼢鼠（SI附录，图S23）外，所有物种的连锁不平衡（LD）（以r2表示）在100kb内衰减至其最大值的一半（SI附录，图S23），与鼹形田鼠属相当。中华鼢鼠在~600kb处降至最低值（SI附录，图S23），显示出最高的LD。
图2.八种鼢鼠的系统发育和遗传多样性。(A)基于全基因组SNP的PCA。(B)基于13个线粒体蛋白编码基因的ML树（左）和核基因组SNP的邻接树（右）。(C)从K=2到10的八种鼢鼠物种的结构分析，其中K=8是最佳拟合模型（SI附录，图S22）。(D)凸颅鼢鼠属物种的分歧参数(FST)和核苷酸多样性(π)。该进化枝（左）由三个适应高海拔的物种组成，其余由四个适应低海拔的物种组成。圆内的值为π(×10-3)。连接两个物种的虚线上的值代表相应的FST。
线粒体和核基因组之间的鼢鼠拓扑不一致。基于线粒体基因组的拓扑与基于核基因组的拓扑不一致。中华鼢鼠是线粒体树中的基础进化枝，而它是核基因组树中的内部进化枝（图2B）。草原鼢鼠和东北鼢鼠聚集成线粒体树中的一个分支。相比之下，草原鼢鼠是一个基础分支，东北鼢鼠在核基因组树中与凸颅鼢鼠属聚集在一起（图2B），可能是因为分布紧密导致渗入。线粒体系统发育中的两个罗氏鼢鼠个体与中华鼢鼠聚集在一起，与基于线粒体DNA中线粒体细胞色素氧化酶亚基1(COI)基因的鼢鼠单倍型网络一致（SI附录，图S24）。
基因树和物种树之间的不一致。我们观察到73.6%的基因树（图3A）与串联物种树（SI附录，图S25）一致。不一致主要存在于两个节点，7和8（图3A和B）。在不稳定节点7中，红色拓扑(4,8j5,6)所占比例为58.1%，其余两个拓扑(4,5j6,8)为蓝色，(4,6j5,8)为紫色。分别为21.1%和20.8%。在节点8中，红色共识树(12,13j4,7)的百分比为63.0%，而其他两个拓扑(12,4j13,7)和(12,7j13,4)分别为21.3和15.7%，分别（图3B-F）。
鼢鼠物种之间的广泛渗入。为了调查基因渗入是否是系统发育不一致的主要原因，使用表示为((P1,(P2,P3)), 银星竹鼠)的有根树计算Patterson的D统计量。从50个三联体中检测到显著的基因渗入（**P<0.01；SI附录，表S13），并计算每个物种对的绝对D平均值（图4A）。通过F3统计分析，没有发现近期杂交的重要证据（SI附录，表S14）。凸颅鼢鼠属和东北鼢鼠之间的基因渗入显著大于凸颅鼢鼠属和草原鼢鼠之间的渗入（***P<0.0001；SI附录，表S13）。这一结果说明了东北鼢鼠和草原鼢鼠的线粒体和核树之间的基因渗入导致的不一致（图2B）。广泛的邻接分布有利于凸颅鼢鼠属和东北鼢鼠之间的渗入，而不是草原鼢鼠之间的渗入。更重要的是，三联体（（草原鼢鼠，东北鼢鼠），中华鼢鼠）的D值为0.098，比中华鼢鼠被始中华鼢鼠中的6种中的任何一种替换时的D值更大（D值约为0.065；SI附录，表S13），表明中华鼢鼠可能是凸颅鼢鼠属最古老的谱系。基因渗入发生在斯氏鼢鼠和除高原鼢鼠1之外的所有其他物种之间。一致地，在罗氏鼢鼠和除甘肃鼢鼠之外的所有物种之间观察到基因渗入（图4A）。这两个结果表明整个QHL的广泛迁移（图1A）。基因渗入仅限于地理上孤立的高原鼢鼠1和除罗氏鼢鼠和东北鼢鼠之外的所有其他物种之间，表明高原鼢鼠1存在异域物种形成（图4A）。
图3.树结构与系统发育树不一致。(A)主导拓扑占据100kb中性位点窗口树的73.6%。黑色块中的数字代表每个节点。进化枝A包含高原鼢鼠和斯氏鼢鼠，而进化枝B包含始鼠中的其余四个物种。(B)使用100-kb中性位点滑动窗口在ASTRAL物种树的分支7和8中诱导的三个二分的频率；标头7和8对应于A中的分支标识。x轴表示每个四重奏拓扑的精确定义。主要拓扑为红色，其他两种替代拓扑为蓝色和紫色。虚线表示1/3阈值。(C和D)节点7中的不一致树结构分别对应于B（左）中的蓝色条和紫色条。(E和F)节点8中的不一致树结构分别对应于B（右）中的蓝色条和紫色条。
渗入和ILS对不一致节点的贡献。系统发育的不一致可能归因于基因渗入或ILS。为了区分基因渗入和ILS，我们通过分支长度(QuIBL)分析对168个拓扑进行量化基因渗入，其中银星竹鼠作为外群。在所有测试的三联体中，84%显示出显着的基因渗入证据（168个三联体中的141个，ΔBIC（贝叶斯信息标准）>10；SI附录，表S15）。四个不一致的拓扑（图3C-F）平均分别归因于8.1%、7.5%、28.8%和8.3%的渗入位点（SI附录，表S15），表明种间渗入是系统发育冲突的主要原因。此外，不一致拓扑的比例相似(fo,(ba1/ba2/sm,ro/ru/ca))(∼28to30%)和(ro/ru/ca,(fo,ba1/ba2/sm)）（∼6.8%至9.1%）表明古代的基因渗入可能发生在中华鼢鼠与斯氏鼢鼠和高原鼢鼠的祖先之间（SI附录，表S15）。
(罗氏鼢鼠(甘肃鼢鼠,秦岭鼢鼠))和(高原鼢鼠2,(高原鼢鼠1,斯氏鼢鼠))的不一致拓扑完全是由ILS引起的（两者ΔBIC<10；SI附录，表S15）。基于进化枝B内部的不一致，我们获得了有根四叉树的15种拓扑，并且主要拓扑与系统发育树一致。前三个不一致的拓扑分别归因于基因渗入、ILS和基因渗入（图4C和SI附录，表S15）。
古代渗入检测。DFOIL统计数据用于使用来自凸颅鼢鼠属的一组四叉树（以银星竹鼠作为外群）来检测古代基因渗入。中华鼢鼠与高原鼢鼠和斯氏鼢鼠的共同祖先之间的古代基因渗入在凸颅鼢鼠属中是最大的，占据了约50%的窗口（图4B和SI附录，表S16），验证了我们基于古代基因渗入的推论。关于QuIBL分析。
为了解决K=2时，贝氏高原鼢鼠和斯氏鼢鼠共同祖先的属内分裂早于鼢鼠和始鼢属间分裂的问题（图2C），我们测试了中华鼢鼠和斯氏鼢鼠之间的古老基因渗入。东北鼢鼠和草原鼢鼠的最近共同祖先。大约59.5%的窗口支持中华鼢鼠和鼹鼠之间的渐渗，与PhyloNet分析一致（图4B和D），比支持中华鼢鼠和高原鼢鼠和E.的最近共同祖先的窗口多∼10%。史密斯（图4B）。因此，中华鼢鼠和鼢鼠之间较大的古代基因渗入有助于它们从K=2到5的聚类（图2C）。
图4.基因渗入检测。(A)块中的值包含每个物种对的平均成对D值。无法在数据集中测试称为姐妹物种对的空块，并且0值不显着。（A，左）基于核基因组的系统发育树。进化枝A由斯氏鼢鼠和高原鼢鼠组成，进化枝B由中华鼢鼠、秦岭鼢鼠、罗氏鼢鼠和甘肃鼢鼠组成。虚线代表DFOIL分析的谱系之间的古代基因渗入。(B)DFOIL推断的基因渗入窗口比例超过30%的Pentataxon系统发育。虚线代表古代基因渗入，如A中所示；这些数字表示支持虚线所示分支之间古代渗入的窗口的比例。(C)以竹鼠为外群的进化枝B所有可能拓扑的平均重量箱线图。Strigulas从下到上分别表示最小值、第一四分位数、中位数、第三四分位数、第四四分位数和最大值。(D)以盲鼹鼠和竹鼠为外类群的鼢鼠亚科最优系统发育网络。
群体统计史。使用成对顺序马尔可夫合并（PSMC）模型估计了八个谱系中每一个谱系的Ne波动。分布在QHL北部的5个谱系（中华鼢鼠、甘肃鼢鼠、秦岭鼢鼠、高原鼢鼠1和高原鼢鼠2）发生了3次种群数量下降（图5A），反映了三个冰川事件：Naynayxungla冰川作用（NG），〜0.5至0.78Mya），倒数第二次冰期（PG，135至〜194千年前；kya）和末次盛冰期（LGM，约20kya）。NG冰期消退后，除高原鼢鼠1和罗氏鼢鼠之外的所有物种均趋向于恢复其Ne（图5A和B）。中华鼢鼠、甘肃鼢鼠和秦岭鼢鼠开始恢复约600至900kya（图5A和SI附录，图S17）。高原鼢鼠2和甘肃鼢鼠在∼0.3至0.5Mya期间出现了小规模的种群扩Zhang，但高原鼢鼠1和罗氏鼢鼠在此期间都没有扩Zhang（图5A和B）。随着海平面上升和大气表面气温下降，所有谱系的Ne在PG和LGM期间均下降。相比之下，分布在QHL以南的斯氏鼢鼠和罗氏鼢鼠在冰期期间扩Zhang（图5A和B），并在冰河期结束后立即收缩（图5C）。
凸颅鼢鼠属分裂为进化枝A和进化枝B～3.22Mya（95%CI为3.07至3.31Mya）（图5C、SI附录、图S26和数据集S1）。后来，高原鼢鼠1和斯氏鼢鼠分别在~1.18Mya（95%CI，1.17至1.34Mya）和∼0.75Mya（95%CI，0.711至0.787Mya）处从进化枝A中分离出来。对于低海拔分支B，中华鼢鼠首先从其祖先∼2.40Mya中分离出来（95%CI为2.21至2.40Mya），证实它是始祖鸾中最古老的物种。甘肃鼢鼠与秦岭鼢鼠∼1.46Mya分叉（95%CI为1.41至1.67Mya）并定居在黄土高原。最后，罗氏鼢鼠与秦岭鼢鼠的分歧约为0.44Mya（95%CI为0.371至0.499Mya）（图5C和SI附录，图S26）。在中华鼢鼠和秦岭鼢鼠和罗氏鼢鼠的最近共同祖先之间检测到基因渗入，与D统计和DFOIL一致（图4A和B）。
图5.群体统计历史。(A)华北种高原鼢鼠、中华鼢鼠、秦岭鼢鼠、甘肃鼢鼠有效种群规模的变动。(B)亚热带种斯氏鼢鼠和罗氏鼢鼠有效种群规模的变动。全球海平面（A和B中的灰色平滑曲线）和大气表面气温（A和B中的浅蓝色平滑曲线）以浅色显示在A和B的背景中。(C)中国鼢鼠的分化和物种形成路线。高原鼢鼠旁边的灰色曲线代表黄河，黄河将高原鼢鼠1和高原鼢鼠2分开。甘肃鼢鼠旁边的工厂代表C4工厂扩建；图中标有“青藏高原隆起”的山峰表示青藏高原隆起，斯氏鼢鼠和罗氏鼢鼠旁边的山峰表示秦岭。
基因组群岛。为了进一步研究物种形成的进化力量，我们从21个物种对中鉴定了8,705个基因组岛(GI)，其FST差异最大为1%（SI附录，图S27和S28）。将相邻岛屿合并后，在整个基因组中总共检测到7,788个GI。我们将GI的DXY、π、Tajima的D、LD和重组率(ρ)与背景区域进行了比较。在地理标志区域观察到LD水平显着升高，ρ水平降低，田岛D下降，表明这些区域存在不同的搭便车现象。地理标志的数量或规模与其种群差异不相关（SI附录，图S30）。
正选择下的基因。尽管凸颅鼢鼠属中有六个物种密切相关，但由于它们对特定栖息地的适应，它们的表型是多样的。为了研究鼢鼠适应和物种形成的遗传基础，我们鉴定了每个物种的正选择基因（PSG）。在分析的20,580个基因中，总共鉴定了747个PSG（拜氏艾球菌286个、史密斯艾球菌115个、丰塔尼艾球菌53个、红毛艾球菌186个、坎苏艾球菌126个和罗斯柴尔德艾球菌91个）。对贝氏菌PSG的基因本体(GO)富集分析显示，有146个富集GO术语，主要涉及肌肉组织（GO:0060537和GO:0007517）、神经发生（GO:0051961和GO:0021766）、心脏形态发生（GO:0003007、GO:0001947、GO:0061371和GO:0007507)、DNA修复(GO:0006281)、自噬(GO:0010508、GO:0010506、GO:0016239和GO:0016241)、免疫学(GO:0030890），以及HMOX1（血红素加氧酶1）活性的调节(R-MMU-9707587)（数据集S2）。至于斯氏鼢鼠，丰富了43个GO术语，其中8个涉及DNA修饰、树突棘发育和血液循环发育（GO：0006305、GO：0003013和GO：0060999）（数据集S3）。其余四个物种，中华鼢鼠、甘肃鼢鼠、罗氏鼢鼠和秦岭鼢鼠，分别丰富了25、21、37和102个GO术语。对于中华鼢鼠，七个术语涉及神经元投射发育（GO：0010975）和细胞凋亡过程（GO：0043065）的调节，以及认知行为（GO：0050890）和学习或记忆（GO：0007611），得到丰富（数据集S4）。甘肃鼢鼠丰富了与肾脏和肾脏系统发育相关的三个术语（GO：0072001和GO：0001822）（数据集S5），而罗氏鼢鼠丰富了与肢体发育相关的六个术语（GO：0035108）以及超过七个涉及胚胎发育(GO:0048568)（数据集S6）。秦岭鼢鼠，一种高海拔物种，丰富了与神经发生(GO:0051960)、免疫学(GO:0002366)、DNA修复(GO:0006281)、心室系统发育(GO:0021591)、呼吸系统(GO:0021591)相关的术语：0044065）、淀粉和蔗糖代谢（mmu00500）、胰岛素受体信号（GO：0009141）、记忆（GO：0007613）和认知（GO：0050890）（数据集S7）。
讨论
鼢鼠在中国的分布。鼢鼠是草食性地下哺乳动物，一生的大部分时间都在地下度过。生活在此类栖息地的生物通常会应对缺氧、黑暗和高碳酸血症(20)。洞穴系统内的氧气浓度取决于透气性和气压(21)。如果土壤湿度高，通风就会受阻，导致隧道内缺氧和高碳酸血症。相反，如果土壤太干燥，有限的植被可能会存活下来，为动物提供数量和质量的食物资源。鼢鼠占据半湿润和半干旱地区（图1A），类似于另一种地下哺乳动物——以色列的鼹形田鼠属的气候生态位（22）。这些栖息地在植被丰富的湿地地区（由于土壤湿度高而导致极度缺氧和高碳酸血症）和通风较好、土壤湿度低且食物资源不足的地区之间提供了折衷方案。
鼹形鼠科三个地下亚科的系统发育。鼹形鼠亚科、鼢鼠亚科和竹鼠亚科被认为是一个科，鼹形鼠科是鼠科的姐妹谱系(23)。然而，三个亚科的系统发育关系仍存在争议。一些研究试图通过形态、核或线粒体基因(24-27)和转录组(28)来阐述系统发育，但产生了相互矛盾的结果（SI附录，图S14）。这种不一致有多种可能性。首先，系统发育是使用一个或仅几个串联分子标记构建的，这些分子标记可能正在选择中。其次，基因渗入或ILS也可能导致树拓扑不一致，从而导致这种不一致。本研究基于全基因组重建的系统发育证实鼢鼠与竹鼠的亲缘关系比盲鼹鼠更为密切，为这种不一致提供了确凿的证据。基因渗入，而不是ILS，对进化枝3中的两个不一致拓扑（鼠属、凸颅鼢鼠属/鼹形田鼠属、竹鼠属）inclade2和（鼠属、竹鼠属/鼹形田鼠属、凸颅鼢鼠属）（SI附录，表S12）负责，而不是ILS。中华鼢鼠和竹鼠属对中华鼢鼠的定位率分别为70.7和95.4%，也表明竹鼠属和凸颅鼢鼠属之间的亲缘关系较近。
八种鼢鼠的种群结构。气候，包括19个变量（SI附录，表S1），可能会影响中国鼢鼠的适应和物种形成（图1B）。降水和温度是影响鼢鼠分布和演化的最有效因素（SI附录，图S18）。第一个PC解释了57.0%的变异，将罗氏鼢鼠与其他物种区分开来（图1B）。罗氏鼢鼠是唯一分布在年降水量超过800毫米地区的物种（图1A）。第二个PC解释了25.2%的变异，并将斯氏鼢鼠和高原鼢鼠与所有其他物种分开，因为它们分布在更高的海拔（SI附录，图S17）。遗传距离与19个气候因素中的11个（SI附录，表S1和S2）和海拔呈正相关（R2=0.520；***P=1.0e-4），这与观察到的高海拔地区之间的巨大差异一致。高海拔和低海拔谱系（图2D）。
鼢鼠的遗传多样性和分化。遗传多样性变异的主要原因可能与群体历史和生态压力有关。由于其最大的Ne（9.4×107）和环境异质性，高原鼢鼠的π在所有物种中最高（SI附录，图S26）。据报道，目前的高原鼢鼠在末次盛冰期期间从四个避难所中恢复了(29)，并且环境异质性增强了多态性的维持(30)。罗氏鼢鼠和斯氏鼢鼠分布在气候温暖的QHL以南，由于温暖的庇护所，它们可能在冰期期间接收了许多迁徙的鸟类。广泛的基因渗入导致了它们第二（罗氏鼢鼠）和第三（斯氏鼢鼠）最高的π（图2D和4A）。罗氏鼢鼠中的π比斯氏鼢鼠中的更大，可能是由于Ne更大以及更多降水施加的更强的应力（图1A和2D以及SI附录，图S26）。中华鼢鼠虽然是基础分支物种，但表现出最低的π（图2D）、最慢的LD衰减率（SI附录，图S23），但不是最小的Ne，这一结果需要在未来进一步研究。
高原鼢鼠1和高原鼢鼠1之间的种内差异。高原鼢鼠2甚至比不同物种对之间的更大（图2E），表明由于黄河屏障，高原鼢鼠分裂成两个不同的物种。高原鼢鼠1遭受年降水量200至400毫米的干旱胁迫（图1A）。相比之下，所有其他物种都分布在年降水量超过400毫米的地区。先前的一项研究报告称，贝氏艾球藻由四个亚群组成(29)，这可能是由于该研究中个体较多且分布区域较大。我们建议高原鼢鼠1和E.由于遗传和表型差异，高原鼢鼠2是两个物种。高原鼢鼠2和高原鼢鼠2之间的分歧。红褐色的数量极少（图2E），可能是由于相似的地质情况（图1A）以及它们之间发生的渗入（图4A和SI附录，图S26）。
系统发育的不一致。面对基因渗入和ILS，使用串联多位点方法存在局限性，这可能会在推断系统发育时造成困难，并且可能无法提供物种的真实历史(2)。节点7和8上发生的不一致（图3B）主要是由于ILS[SI附录，表S15；拓扑（ba2，（ba1，sm））和（ca，（ro，ru））]。STRUCTURE分析中，秦岭鼢鼠、甘肃鼢鼠和罗氏鼢鼠从K=4到6的聚类与斯氏鼢鼠、高原鼢鼠1和高原鼢鼠2从K=5到6的聚类不一致（图2C)可归因于ILS。ILS在快速分化的谱系中无处不在，例如由于物种形成事件之间的时间较短而发生的适应性辐射(31-34)。因此，鼢鼠在这两个节点上相对较短的发散时间可能导致了这种不一致。较大的Ne，7.9×107（高原鼢鼠和斯氏鼢鼠的祖先）和6.3×107（甘肃鼢鼠、秦岭鼢鼠和罗氏鼢鼠的祖先），使得这两个分支具有高度的ILS（SI附录，图S26)(35)。
群体统计历史。气候塑造了鼢鼠（图2B）和盲鼹鼠（22）的适应和物种形成。虽然它们都生活在温度相对恒定的地下，但它们却受到地上温度的影响（36）。冰川作用导致鼢鼠被冻死，并导致为食草鼢鼠提供食物的植被减少。冰川作用期间植被有限，解释了QHL以北五个物种Ne的减少。然而，在QHL以南气候温暖湿润的情况下，植被预计会更好，这可以为鼢鼠提供庇护所。PG和LGM后斯氏鼢鼠和罗氏鼢鼠Ne的下降可能是由于中国南部和东部的大雨造成的(37,38)。检测到从高原鼢鼠2向斯氏鼢鼠以及从秦岭鼢鼠向罗氏鼢鼠的基因渗入（图3D），这与鼢鼠在冰期增加Ne（图5B）和π期间迁移到更温暖的庇护所的假设是一致的。(图2D)的斯氏鼢鼠和罗氏鼢鼠(图2E)。这两个物种的物种形成可能是由冰川驱动的鼢鼠向南迁徙引发的。3个物种(秦岭鼢鼠、中华鼢鼠和甘肃鼢鼠)分布区内的温度和植被因海拔不同而不同。秦岭鼢鼠分布在海拔最高的地方，相应地，New秦岭鼢鼠也是最后恢复的。另一方面，中华鼢鼠出现在所有三个物种中海拔最低的地方（SI附录，图S17），相应地，它的Ne首先扩展（图5A和SI附录，图S17）。
与高原鼢鼠1相比，NG后高原鼢鼠2的种群扩Zhang较小，因为高原鼢鼠2分布在年降水量大于400毫米的地区。NG冰期后，在甘肃鼢鼠和高原鼢鼠2中观察到种群扩Zhang（图5A和B），这与夏季降雨量增加但总降水量减少时黄土高原0.43Mya的C4植物扩Zhang一致（39）。C4植物适应炎热、干燥的环境，鼢鼠主要以C4植物(40)为食，这使得它们能够同步扩Zhang约∼0.40至0.50Mya。
积极选择的基因。在适应新环境和形成新环境的过程中，基因总是会被选择。所有生活在地下的鼢鼠都要应对黑暗、缺氧、高碳酸血症和充满病原体的环境(41)。然而，鼢鼠的环境在海拔（SI附录，图S17）、温度、降水（SI附录，表S1）和植被方面有所不同。贝雷伊蚊生活在海拔2800～4600m之间，其特点是气温低、气压低造成低氧、UV（紫外线）强。选择与免疫学相关的基因是因为潮湿地下隧道中的病原体威胁哺乳动物的健康（20）。鼢鼠的大脑和神经系统预计会在缺氧下受损，因此发达的神经系统可以促进它们对缺氧的适应，从而解释在神经发生中丰富的PSG(42)。PSG在细胞周期和DNA修复中富集，是其抵抗高海拔氧化损伤和强紫外线的基础(43)。此外，高原鼢鼠还可以通过调节HMOX1的表达和活性来很好地保护自身免受缺氧引起的超氧化物的影响(44)。与自噬相关的基因可以增加高海拔地区缺乏食物时的能量利用率(45)。正如预期的那样，多个基因在心脏发育和确定心脏左/右不对称性的功能上得到丰富，有利于氧气输送（46），有助于缺氧适应。高海拔地区的植被稀缺且质量低下，因此，为了寻找更多的食物来源，鼢鼠必须探索更大的区域（47）。相应地，与海马体发育和记忆行为相关的基因也在选择中。高原鼢鼠的PSG的肌肉组织也丰富，这可能是由于它们广泛的挖掘活动（数据集S2）。位于第二高海拔的斯氏鼢鼠丰富了许多与高原鼢鼠类似的GO术语，如神经系统发育和DNA修饰。此外，可以选择与血液循环相关、促进氧气输送的基因来应对缺氧。位于海拔较高的秦岭山区的红毛桉，富含负责氧气吸入和输送的心室系统发育和呼吸系统发育的基因，有利于其适应缺氧。淀粉和蔗糖代谢途径的转变可以促进脂质的积累，帮助鼢鼠在没有新食物供应的情况下度过漫长的冬季（48）。与胰岛素受体信号相关的基因可能会增加葡萄糖吸收以产生热量(49)，这可能有助于鼢鼠应对高海拔地区的低温。与记忆相关的PSG有利于记住复杂的地下隧道。对于甘肃鼢鼠，由于黄土高原干旱的环境，选择了与肾系统相关的基因来增加吸水性（50）。在中华鼢鼠中鉴定出与细胞凋亡相关的基因，这反映了其抗肿瘤潜力(51)。通常，罗氏鼢鼠体型较小，位于平均气温较高的QHL以南。相比之下，贝氏巨蜥体型较大，且位于青藏高原海拔较高、平均气温较低的地区。两者都丰富了生长调节和肢体发育的术语，这与伯格曼规则（52）是一致的。
基因组群岛。每个鼢鼠物种对在基因组中的差异是异质的，并且可能产生GI（SI附录，图S28）。GI经常出现在鼢鼠基因组中，特别是在密切相关的物种之间。局部适应(53)、生殖隔离(54)、选择(55–57)以及古代平衡多态性的排序(53,58)可能会产生地理标志。在本研究中，岛屿内的FST和DXY均大于背景区域（SI附录，表S17），这可以通过古代多态性的发散排序（53,58,59）或渗入的物种形成（60）来解释。如果后一种模型诱导了GI，则预计会出现具有增加的分歧时间和基因渗入的更稳健的岛屿(61,62)。然而，不同物种对的GI的分化没有显著差异，并且GI的大小或数量与渗入和分化时间不相关（SI附录，图S30），拒绝了物种形成与基因流模型。
所有物种对的GI均以LD升高、π水平降低和重组率降低为特征（SI附录，表S17）。我们观察到大多数物种对之间的GI的分化率和重组率之间存在显著相关性（SI附录，图S28），这表明可能发生分化搭便车或背景选择（53,58,60）。如果地理标志是通过连锁选择启动的，岛屿地区内的DXY不会增加(53,60)。相比之下，岛屿地区内的DXY在我们的结果中确实有所增加（SI附录，表S17），否定了连锁选择导致GI的假设。DXY升高表明，包含GI的单倍型可能比所比较的所有物种对的其余基因组更早地进行遗传分离（SI附录，表S17；曼-惠特尼U检验，***P<2.2e-16）。这些结果表明，古代多态性的发散排序产生了GI（53、58、60、63）。
鼢鼠分歧和物种形成动量。内在（遗传学、群体历史）和外在（气候、栖息地）因素可能会产生物种形成动量（64）。凸颅鼢鼠属与鼢鼠分离超过4Mya，化石资料证实[8]，凸颅鼢鼠属在青藏高原西南约3.22Mya处形成了物种（图5D），当时正值昆-黄构造运动发生，青藏高原已经隆起。（65）。高原鼢鼠1和高原鼢鼠2之间的分离发生了1.18Mya（图5C），这与导致黄河渗透青藏高原约1.2Mya的构造运动一致（66），将这两个种群分开（图1A和5C））。高原鼢鼠2向南形成了NG周围的斯氏鼢鼠0.75Mya（图5C）和斯氏鼢鼠的Ne在PG和LGM期间都有扩Zhang（图5B）。因此，冰川作用被认为是物种形成动量，驱使动物向南迁徙到QHL以南的温暖避难所。秦岭鼢鼠向北形成甘肃鼢鼠1.46Mya，向南形成罗氏鼢鼠0.44Mya（图5C）。甘肃鼢鼠的向北扩散可能是由黄土高原南部1.6Mya发生的热适应C4植物扩Zhang引发的(39)。罗氏鼢鼠在NG期间向南形成物种（图5B和C），其种群在末次盛冰期期间扩Zhang，表明冰川作用将鼢鼠驱赶到温暖的南部避难所（图5B），在那里形成了物种。从高原鼢鼠2向南迁移到斯氏鼢鼠，以及从秦岭鼢鼠向南迁移到罗氏鼢鼠，与检测到的基因渗入（图4A和D）和更高的遗传多样性（图2D）一致。尽管高原鼢鼠的分布与甘肃鼢鼠相邻（图1A）（67），但它是次要接触，最近的渗入受到限制。气候（SI附录，图S18）、构造运动、青藏高原的隆升和C4植物的扩Zhang引发了中国鼢鼠的物种形成。
材料与方法
使用PacBio Sequel平台进行长读长测序，而使用Illumina NovaSeq 6000和DNBSEQ-T7平台进行短读长测序。BWA将干净的读数映射到染色体水平参考基因组，并使用GATK调用SNP。使用ArcGis绘制采样图，并使用promp对19个气候因子进行主成分分析。使用MCMCTree估计鼢鼠亚科的分化时间。群体基因组分析（包括系统发育树、PCA和STRUCTURE构建）分别由Phylip、GCTA和ADMIXTURE进行。通过VCFtools计算遗传分化参数，例如FST和π。使用RAxML、ASTRAL和DiscoVista进行不一致树分析。通过D统计、QuIBL和DFOIL分析最近的基因渗入、不完整的谱系排序和古老的基因渗入。使用PSMC和fastsimcoal v2.7推断群体统计历史。详细信息请参阅SI附录“材料和方法”。
期刊：PNAS
文章标题：Genomic insights into zokors’ phylogeny and speciation in China
作者信息：Xi Liu, Shangzhe Zhang, Zhenyuan Cai, Zhuoran Kuang, Na Wan, Yinjia Wang, Leyan Mao, Xuan An, Fang Li, Tuo Feng, Xiaolong Liang, Zhenglei Qiao, Eviatar Nevo nevo@research.haifa.ac.il, and Kexin Li likexin@lzu.edu.cn-10Authors Info & Affiliations
原文链接：https://doi.org/10.1073/pnas.2121819119
文内图片及封面图片来源原文