太阳花马基因流模式与适应性渗入特征分析

Gene flow patterns and adaptive introgression characterization in Taiyanghua horses

摘要

为系统鉴定太阳花马在杂交育种过程中的基因流模式与适应性渗入区域，基于全基因组重测序SNP数据，对太阳花马（TYH，n=31）、蒙古马（MON，n=31）和纯血马（TB，n=32）进行群体遗传结构、基因流（f₃、D统计量、窗口化D）及选择信号（Fst、π ratio、XP-CLR、iHS、XP-EHH）整合分析。结果表明：1）太阳花马与纯血马遗传分化（Fst=0.040 7）显著低于其与蒙古马（Fst=0.068 5），ADMIXTURE分析显示纯血马祖源成分均值为76.5%；D统计量证实存在由纯血马向太阳花马方向的显著基因流（D=0.221 6，Z=29.1），渗入信号分散分布于全基因组。2）渗入-选择共定位整合分析识别出推定受选择的纯血马渗入区域168个窗口（126个基因），功能富集于免疫调节（IL-17信号通路）和神经信号通路。3）推定受选择的蒙古马祖源保留区域841个窗口（441个基因），富集运动步态与神经肌肉功能，其中包含马步态调控关键基因DMRT3。本研究在太阳花马中实现基因流扫描与多方法选择信号的系统整合，揭示杂交育种中纯血马渗入与蒙古马祖源保留并存的分子适应特征。

Abstract

To systematically identify the gene flow patterns and adaptive introgression regions of Taiyanghua Horses during crossbreeding， an integrated analysis of population genetic structure， gene flow （f₃-statistic， D-statistic， windowed D-statistic） and selection signatures （Fst， π ratio， XP-CLR， iHS， XP-EHH） was performed based on genome-wide SNP data from whole-genome resequencing in Taiyanghua Horses （TYH， n=31）， Mongolian Horses （MON， n=31） and Thoroughbred Horses （TB， n=32）. The results showed that： 1） The genetic differentiation （Fst=0.040 7） between TYH and TB was significantly lower than that between TYH and MON （Fst=0.068 5）， and ADMIXTURE analysis revealed a mean Thoroughbred ancestry proportion of 76.5%. The D-statistic confirmed significant gene flow from TB to TYH （D=0.221 6， Z=29.1）， with introgression signals dispersedly distributed across the whole genome. 2） The introgression-selection co-localization integrated analysis identified 168 windows （covering 126 genes） in putative selected TB introgression regions， which were functionally enriched in immune regulation （IL-17 signaling pathway） and neural signaling pathways. 3） It was estimated that 841 windows （covering 441 genes） were detected in putative selected regions with retained Mongolian horses ancestry， which were enriched in locomotor gait and neuromuscular function， including DMRT3， the key gene regulating equine gait. This study achieves the systematic integration of gene flow scanning and multi-method selection signature analysis in TYH， and reveals the coexisting molecular adaptive characteristics of Thoroughbred introgression and Mongolian horses ancestry retention during crossbreeding.

Graphical abstract

关键词

太阳花马 / 种群遗传结构 / 基因流 / 选择扫描 / 适应性渗入

Key words

Taiyanghua horses / population genetic structure / gene flow / selective sweep scanning / adaptive introgression

马在人类社会中长期承担交通运输、军事乘骑、竞技比赛和休闲骑乘等功能，其遗传改良过程具有显著的人工选择特征。纯血马是全球最具经济价值和影响力的赛马品种之一，经过长期定向选育，形成了鲜明的速度型遗传特征^［1］。在分子层面，MSTN位点变异与不同比赛距离偏好密切相关，其中g.66493737T>C多态性对应的不同基因型与短途冲刺能力和耐力表现显著相关^［2］。与之相对，蒙古马及其衍生类群在草原环境中形成了较强的生境适应能力。已有研究表明，蒙古马在天然草场放牧条件下能够维持较稳定的泌乳与代谢状态，并对季节性草场波动和低投入饲养环境表现出良好适应^［3］；高负荷运动训练前后转录组分析显示，其骨骼肌差异表达基因显著富集于肌肉结构发育、心血管系统发育、肌肉收缩和能量代谢等通路^［4］。群体基因组研究进一步在中国蒙古马不同类群中检测到与区域适应相关的选择信号，并在乌审马和百岔铁蹄马中分别发现与步态、持续运动能力及蹄部健康相关的候选区域^［5］。基于全基因组测序的正向选择分析亦鉴定出DMRT3、MAPK13、CCN2和FSCN3等与运动能力相关的候选基因^［6］，蒙古马遗传多样性的群体特征也已得到系统评估^［7］。在马品种改良实践中，如何在引入外来品种速度优势的同时保留地方品种环境适应性，是杂交育种面临的核心遗传问题。

太阳花马，又称“草原纯血马”或“新锡林郭勒马”，是以蒙古马系统母本为基础、纯血马（Thoroughbred）为主要导入父本，经多代级进杂交和持续选育形成的草原改良乘用马品系^［8］。其母系来源为锡林郭勒马，该品种于1952—1985年间在蒙古马基础上导入苏高血、顿河马和卡巴金马等血统育成的乘挽兼用型地方良种^［9］。太阳花马的培育目标是在保留草原马环境适应力、耐粗饲性和耐力优势的基础上，引入纯血马的速度潜力、运动能力和竞赛性能，形成兼具速度与耐力特征的复合型改良乘用马。培育始于1995年，自1997年起引进纯血公马与锡林郭勒母马进行级进杂交，先后经历杂交导血、横交固定和选育提高3个阶段，经近30年的持续选育，截至2024年种群规模已达约2 500匹^［8］。围绕太阳花马的基因组研究已取得一定进展，但总体上仍以遗传结构描述和选择信号筛查为主。已有研究从遗传结构、纯合子片段（ROH）、群体分化（Fst）、核苷酸多样性（π）及拷贝数变异（CNV）等层面对太阳花马进行了初步解析，结果一致表明其整体遗传结构更接近纯血马，并鉴定出ATF2、IGFN1、PPARGC1A和SLC25A15等与运动性能、肌肉功能和能量代谢相关的候选基因^［10-13］；同时发现育种过程中伴随遗传多样性下降与近交积累的动态变化^［11］。

现有研究至少还存在三方面不足。第一，虽然太阳花马的杂交育种背景明确，但缺少对纯血马导入方向、导入强度及其基因组分布模式的系统量化评估。第二，既往选择扫描多依赖单一方法或有限方法组合，难以区分稳定信号与方法特异性噪声。第三，尚缺乏将基因流信号与选择信号联动起来的分析框架，因而难以识别那些既来源于外来祖源又在目标群体中受到保留或强化的候选区域。

基于此，本研究以太阳花马为核心对象，以蒙古马和纯血马为参照群体，在统一质控和分析框架下系统开展群体遗传结构分析、基因流检测、多方法选择扫描及其整合分析。鉴于锡林郭勒马与蒙古马系统之间具有明确的谱系连续性，而同批次锡林郭勒马全基因组数据尚未纳入本研究，故以蒙古马作为母系祖源的近似参照，用于表征锡林郭勒马中延续自蒙古马系统的基础适应性遗传成分。重点回答以下问题：太阳花马与双亲祖源之间的遗传关系如何；纯血马导入信号在全基因组范围内呈现何种分布格局；哪些外来渗入区域或地方祖源保留区域在杂交育种过程中可能受到选择。研究结果可为导血育种分子遗传效应的解释、太阳花马育种策略优化及地方马遗传资源保护提供依据。

1　材料与方法

1.1　样本与测序

太阳花马（TYH，n=31）横交固定，其中♂27匹，♀4匹，均来自锡林郭勒盟种马场。样本于2022年7月采集于内蒙古锡林郭勒盟乌珠穆沁旗，采样过程遵循当地动物福利法规，所有操作经中国农业大学动物伦理委员会批准（批准号：AW72203202-1-1）。基因组DNA采用华大基因DNA抽提试剂盒提取，全基因组重测序文库通过BGI Optimal DNA Library Prep Kit构建，在DNBSEQ测序平台开展双末端150 bp（PE150）测序。蒙古马（MON，n=31）和纯血马（TB，n=32）的基因组数据来源于美国国家生物技术信息中心（National center for biotechnology information， NCBI）数据库。本研究所有遗传分析仅基于常染色体（chr1-31）位点。

1.2　数据处理与变异检测

原始测序数据使用fastp（v0.23.4）进行质控，采用默认参数过滤低质量reads。质控后的reads使用GTX.CAT（v2.1.1）比对至马参考基因组EquCab3.0（GCF_002863925.1），生成BAM文件。使用GATK（v4.4.0）HaplotypeCaller模块对每个样本进行变异检测并生成gVCF文件，随后采用GLnexus进行群体联合分型，获得联合分型VCF文件。原始VCF经GATK硬过滤流程（QD<2.0、QUAL<30.0、SOR>3.0、FS>60.0、MQ<40.0、MQRankSum<-12.5、ReadPosRankSum<-8.0）进行初步质控。

1.3　样本与位点质控

使用PLINK（v1.90b7.7）^［14］对样本间亲缘关系进行评估。首先对SNP进行LD剪枝（-indep-pairwise 50 10 0.2），再基于剪枝后的SNP集合计算血缘同源（Identity by descent，IBD）矩阵（-genome）。后续分析纳入全部94个样本（TYH=31，MON=31，TB=32）。

在位点质控层面，使用bcftools（v1.17）提取双等位SNP位点，并以AC≥1保留样本中的全部多态位点。为确定质量过滤阈值，对前100 000个位点的QUAL与AQ分布进行统计，分别比较1%、5%和10%分位数阈值（QUAL/AQ≥150、≥505和≥676）对全量AC≥1数据集的过滤效果。3个阈值下保留位点数依次为20 465 395、18 304 914和16 639 843个，对应Ti/Tv比值分别为2.134、2.166和2.175，平均位点缺失率分别为0.058 6%、0.053 3%和0.053 0%。在低质量区间（QUAL<150、150~505和505~676）内，AC≤2位点的比例分别为62.6%、49.0%和44.6%，而在QUAL/AQ≥676的保留集中降至9.14%，表明1%与5%分位数阈值虽可保留更多位点，但低频噪声比例较高。综合Ti/Tv提升、低频噪声剔除效果、缺失率稳定性及位点保留规模，最终采用10%分位数对应的阈值QUAL≥676且AQ≥676进行硬过滤，获得16 639 843个高质量SNP，Ti/Tv比值为2.175，平均位点缺失率为0.053%，用于后续分析。

1.4　种群遗传结构分析

基于高质量SNP集合进行LD剪枝（PLINK，-indep-pairwise 50 10 0.2，-chr-set 31），仅使用常染色体（chr1-31）位点，获得约2.8 M独立SNP用于结构分析。主成分分析（Principal component analysis，PCA）采用PLINK（-pca 20），提取前20个主成分并计算各主成分解释方差。群体间遗传距离的置信区间通过SNP水平bootstrap估计：每次从LD剪枝后约2.8 M独立SNP中有放回地抽取等量位点子集，重新计算PCA并提取各群体的PC1-PC2质心坐标，重复1 000次获得质心距离的经验分布及95%置信区间。ADMIXTURE分析（v1.3.0）^［15］在K=1~10范围内估计种群结构，采用10折交叉验证（-cv=10）比较不同K值，选择交叉验证误差（Cross-validation error，CV error）最低的K值作为最优K。邻接树（Neighbor-joining tree，NJ树）基于LD剪枝后SNP计算个体间状态同源（Identity by state，IBS）距离（PLINK，-distance 1-ibs），使用R包ape构建。

1.5　种群分化分析

采用VCFtools（v0.1.16）^［16］计算3种群两两间的全基因组加权Weir & Cockerham Fst^［17］（–weir-fst-pop）。滑窗Fst采用100 kb窗口、50 kb步长（-fst-window-size 100000 -fst-window-step 50000）进行计算。

1.6　基因流检测

全基因组基因流检验采用f₃和D统计量2种方法：f₃统计量（3-population test）检验形式为f₃（TYH；MON，TB），使用ADMIXTOOLS（qp3Pop）^［18］计算，通过block-jackknife（block size=5 Mb）估计标准误差；D统计量（ABBA–BABA检验）以驴（Equus asinus，n=16）为外群，设置两组拓扑结构D（TYH， MON； TB， Donkey）与D（TYH，TB；MON，Donkey），采用ADMIXTOOLS（qpDstat，block-jackknife，block size=5 Mb）计算D值与Z值，以|Z|>3为显著性判定标准。

窗口化D统计量分析采用Dsuite（v0.3 r24）软件，基于D（TYH， MON； TB， Donkey）拓扑结构，窗口大小100 kb、步长50 kb，计算每个窗口内的ABBA/BABA位点数和D值（D>0表示纯血马渗入方向，D<0表示蒙古马背景保留），以Z分数阈值（|Z|>2.576对应P<0.01，|Z|>3.291对应P<0.001）标记显著窗口。

1.7　选择信号检测

本研究采用5种选择分析方法，统一使用top 5%分位数阈值进行筛选。该阈值为群体基因组学选择扫描研究的常用经验阈值，其在单方法层面保留较宽松的候选范围以避免遗漏真实信号，为降低方法特异性噪声并提高候选区域稳健性，识别结果进一步采用多方法交叉验证（n_methods≥3）进行约束，从方法一致性层面严格控制假阳性^［12］。

频率侧方法：1）滑窗Fst，以太阳花马为目标种群，分别与蒙古马、纯血马计算滑窗Fst（VCFtools，100 kb窗口、50 kb步长），以各比较的top 5%阈值筛选高分化窗口（TYH vs MON： Fst≥0.162 1；TYH vs TB： Fst≥0.113）；2）π ratio，使用VCFtools计算滑窗核苷酸多样性（-window-pi 100000-window-pi-step 50000），计算πTYH/πMON与πTYH/πTB，以bottom 5%阈值筛选低多样性窗口；3）跨群体复合似然比（Cross population composite likelihood ratio，XP-CLR），使用XP-CLR（v1.0）^［19］进行跨种群选择扫描（-size 50000 -step 25000），以top 5%阈值筛选候选窗口。

单倍型方法：首先使用Beagle 5.4^［20］（beagle.22Jul22.46e.jar）按染色体对常染色体SNP进行相位推断。4）整合单倍型得分（Integrated haplotype score，iHS），使用selscan（v1.2.0a）^［21］计算iHS并标准化（-bins 100），以|iHS_std|的top 5%阈值筛选显著位点，按100 kb窗口汇总，以窗口内显著位点占比>5%界定显著窗口；5）跨群体扩展单倍型纯合性（Cross population extended haplotype homozygosity，XP-EHH），使用selscan计算太阳花马相对于蒙古马、纯血马的XP-EHH（-max-gap 200000-max-extend 1000000），标准化后同样按100 kb窗口汇总。

1.8　多方法整合与候选基因注释

选择信号整合采用两级策略：首先在每种方法内部，将TYH vs MON与TYH vs TB两组比较的候选窗口取并集，获得该方法的综合候选窗口集；其次在5种方法间，统计每个窗口被识别为候选的方法数（n_methods，范围0~5），以n_methods≥3作为高置信候选阈值。将候选窗口与EquCab3.0基因注释（GTF）比对提取覆盖基因。

功能富集采用DAVID在线工具^［22］（马物种默认背景）进行基因本体论（Gene ontology，GO）与京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）注释，以修正Fisher精确P<0.05为初筛标准，同时进行Benjamini-Hochberg多重检验校正，对基因集规模较小的子集，多重检验校正力度较大，结果以未校正P值结合富集倍数综合解读。

1.9　适应性渗入分析

将窗口化D统计量与选择信号进行整合分析。选择信号仅纳入TYH vs MON的比较结果，该策略与人类适应性渗入研究的经典范式一致^［23］。基因流信号强度界定采用D值分布的Mean±2SD阈值（而非基于Z分数的P值阈值）。整合定义：Adaptive（+）为D>mean+2SD且n_methods≥3的窗口，代表高纯血马渗入信号与选择信号的共定位；Adaptive（-）为D<mean-2SD且n_methods≥3的窗口，代表低渗入背景下的选择信号。将候选区域与基因注释比对提取覆盖基因，使用DAVID进行功能富集分析。

2　结果与分析

2.1　太阳花马种群遗传结构与分化

PCA分析显示3种群在主成分空间中呈现明显分离（图1（a））。PC1和PC2分别解释13.51%和6.56%的遗传变异，基于1 000次bootstrap的PC1-PC2质心距离估计表明，TYH-TB距离为0.055（95% CI： 0.047~0.063），显著小于TYH-MON距离0.184（95% CI： 0.175~0.203），且两者置信区间不重叠。进一步基于前4个主成分坐标进行PERMANOVA检验，结果显示群体效应显著（R²=0.249，P=0.001），支持3种群在多维遗传空间中的稳定分离。

ADMIXTURE分析进一步揭示了太阳花马的混合型遗传背景（图1（b））。在K=2时交叉验证误差最低。太阳花马呈现典型的双祖源混合结构，其纯血马祖源成分均值为76.5%，范围为47.6%~97.1%，该跨度反映了群体中不同级进杂交代次个体的共存（F1理论值≈50%，F2≈75%，F3≈87.5%），与选育提高过程中不同个体间导血程度与杂交代次存在连续差异的育种状态相符。结合PCA结果可见，TYH个体主要分布于MON与TB之间，并沿纯血马方向形成连续梯度；但在主成分空间中，个别样本仍表现出轻微聚集趋势，且可能叠加一定的家系或育种批次效应。全基因组加权Fst分析进一步表明（图1（c）），TYH与TB的遗传分化最低（Fst=0.040 7），TYH与MON次之（Fst=0.068 5），MON与TB 最高（Fst=0.119 2）。

2.2　太阳花马种群基因流与杂交证据

基于常染色体SNP对3种群进行f₃统计量检验。结果显示f₃（TYH； MON， TB）=0.055，标准误差为0.002 37，Z值为23.18（SNPs=13 848 125）（图2（a））。f₃值为正，未检测到经典的混合信号（admixture f₃<0）^［18］，该结果可能由以下因素单独或共同导致：一是太阳花马经多代级进杂交后遗传漂变的积累削弱了f₃对混合信号的检测灵敏度；二是本研究以蒙古马而非锡林郭勒马作为母系祖源的代理群体，与真实祖源存在一定偏离。

为进一步验证基因流方向，本研究以驴（n=16）为外群进行Patterson's D统计量分析。结果显示：D（TYH， MON； TB， Donkey）=0.221 6（D值表示等位基因共享的偏离程度），Z=29.1（极显著），BABA位点（611 329）显著多于ABBA位点（389 513），表明太阳花马与纯血马共享更多衍生等位基因，存在纯血马向太阳花马方向的显著基因流；而D（TYH， TB； MON， Donkey）=-0.001 6，Z=-0.903不显著，BABA与ABBA位点数相近（388 265 vs 389 513），表明未检测到蒙古马向太阳花马的额外渗入（图2（b））。该结果与太阳花马的育种模式一致：太阳花马中的蒙古马祖源成分（ADMIXTURE估计约23.5%）源自锡林郭勒马母本的基础遗传背景，属于受体群体的共祖遗传而非独立的基因流事件，因此D统计量仅检测到纯血马导入方向的显著信号。

为探索基因流在基因组上的分布特征，本研究基于D（TYH， MON； TB， Donkey）拓扑结构进行窗口化D统计量分析（100 kb窗口，50 kb步长）。在该框架下，D>0表示纯血马方向的基因渗入（TYH与TB共享更多衍生等位基因），D<0表示蒙古马遗传背景的保留。全基因组45 099个窗口中，纯血马渗入方向（D>0）检测到139个显著窗口（P<0.01）和5个极显著窗口（P<0.001）；蒙古马背景方向（D<0）检测到365个显著窗口和25个极显著窗口（图2（c））。纯血马渗入信号分散分布在22条染色体上，未形成大规模连续热点区域，提示渗入是全基因组范围的弥散过程。

显著渗入窗口中包含多个与马运动性能相关的候选基因：脂肪酸代谢相关的ACACA（乙酰辅酶A羧化酶，曾被Gu等^［24］在纯血马选择扫描中鉴定）和PNPLA2（脂肪甘油三酯脂肪酶），cAMP信号通路的ADCY1，以及毛色表型相关的OCA2、KITLG等基因。值得注意的是，chr7：38-50 Mb区域表现为全基因组最显著的蒙古马背景保留区域（90个显著窗口，P<0.01），该区域富集了神经肌肉信号传导（NTM、OPCML、CACNA1A）和能量代谢（PRKACA）相关基因，提示这些区域在杂交育种过程中抵抗了纯血马的基因渗入，保留了蒙古马的遗传背景。

2.3　选择信号检测与候选基因

采用5种选择分析方法对太阳花马进行选择扫描，各方法在方法内整合（TYH vs MON与TYH vs TB取并集）后识别候选窗口。韦恩图（图3（b））展示5种方法的交集关系：Fst、π ratio、XP-CLR、iHS、XP-EHH各自贡献独特候选区的同时存在大量重叠。以≥3种方法支持作为高置信阈值，共识别1 760个候选窗口（1 740个基因）；其中5种方法支持的阈值共识别72个核心窗口，对应103个核心候选基因（图3（a））。

对1 740个候选基因进行功能富集分析（图3（c）和（d）），共获得351条显著结果（P<0.05）。KEGG通路中，运动蛋白通路（Motor proteins，P=9.18×10⁻⁶，校正P=0.002，富集倍数2.72）最为显著，经多重检验校正后仍高度显著，涉及动力蛋白重链DNAH2、驱动蛋白KIF家族等分子马达基因；细胞衰老通路（Cellular senescence，P=3.23×10⁻⁵，校正P=0.003）涉及乙酰胆碱受体CHRNB1、心肌肌球蛋白轻链MYL2等神经肌肉和心肌功能基因。GO-MF层面，锌离子结合（P=7.13×10⁻⁸，校正P=6.79×10⁻⁵）和C2H2型锌指结构域（P=8.38×10⁻⁷，校正P=3.02×10⁻⁴）高度富集，涉及ZNF家族转录因子。此外，核心候选区还包含神经信号传导基因（OPCML、NAV1、RBFOX3）和离子转运基因（ATP1B2）。上述结果表明，太阳花马在形成过程中，肌肉收缩、神经信号传导和能量代谢相关区域受到选择压力。

2.4　适应性渗入区域鉴定

将基因流信号与选择信号（TYH vs MON）整合分析（图4（a）），识别出168个Adaptive（+）候选窗口（126个基因）和841个Adaptive（-）候选窗口（441个基因）。对126个Adaptive（+）基因进行功能富集分析（图4（b）），免疫信号通路富集较为突出：白细胞介素-17信号通路（IL-17 signaling pathway，P=2.06×10⁻³，富集倍数9.05）、Fc ε RI信号通路（Fc epsilon RI signaling pathway，P=5.83×10⁻³，富集倍数10.69）、巨噬细胞活化（Macrophage activation，P=3.30×10⁻³，富集倍数34.39）和细胞因子活性（Cytokine activity，P=1.85×10⁻²，富集倍数4.92），涉及IL4、IL13、IL5、IL17C等白细胞介素基因。由于Adaptive（+）基因集较小（126个），上述通路经Benjamini-Hochberg校正后未达0.05阈值，但考虑到多条免疫相关通路一致性富集且富集倍数较高，提示纯血马渗入并受选择的基因片段可能与免疫调节功能相关。此外，γ-氨基丁酸能突触通路（GABAergic synapse，P=1.42×10⁻²）和生长激素合成、分泌与作用通路（Growth hormone synthesis， secretion and action，P=2.95×10⁻²）亦有富集。

对441个Adaptive（-）基因进行功能富集分析（图4（c）），运动与步态相关功能富集较为突出。GO-BP层面，微管运动通路（Microtubule-based movement，P=1.46×10⁻⁴，富集倍数5.87）涉及动力蛋白重链DNAH2、DNAH12和驱动蛋白KIF23、KIF5C、KIF18A等分子马达基因；成体运动行为（Adult locomotory behavior，P=1.54×10⁻²，富集倍数5.20）与运动行为（Motor behavior，P=1.07×10⁻²，富集倍数8.63）。DM DNA结合结构域（DM DNA-binding domain，P=1.21×10⁻²，富集倍数17.45）在功能聚类中与运动行为和性别分化构成高度关联簇（聚类富集得分=1.87），该簇中核心基因DMRT3（Doublesex and mab-3 related transcription factor 3）是已知的马步态调控关键基因^［25］。KEGG通路中，吗啡成瘾通路（Morphine addiction，P=1.97×10⁻⁴，校正P=0.052，富集倍数5.57，涉及GABA受体和cAMP信号）接近FDR校正显著阈值，运动蛋白通路（Motor proteins，P=6.90×10⁻³，富集倍数2.93）亦有富集。此外，转录调控、角蛋白/表皮形态发生和染色质调控等功能亦富集于Adaptive（-）区域。上述多个功能维度的一致性富集提示蒙古马遗传背景在运动步态调控和神经肌肉功能方面发挥重要作用。

3　讨 论

遗传结构分析一致表明太阳花马在遗传上更接近纯血马（Fst=0.040 7），ADMIXTURE分析显示其纯血马祖源成分均值达76.5%，D统计量进一步证实存在由纯血马到太阳花马方向的显著基因流，三项结果共同印证了以纯血马为父本的级进杂交育种史。Petersen等^［26］对全球马品种的分析表明纯血马因长期封闭育种而遗传多样性较低，本研究中太阳花马纯血马祖源成分均一性较高的观察亦反映了纯血马父本遗传背景的均质性。窗口化D统计量显示纯血马渗入信号呈全基因组弥散分布，未形成明显连续热点（图2（c）），符合多代回交导致供体基因组片段被逐渐打碎，分散整合至受体基因组的预期模式。在基因流方向上，D统计量仅检测到由纯血马向太阳花马方向的显著信号，而蒙古马方向不显著。ADMIXTURE估计太阳花马保留约23.5%的蒙古马祖源成分，该成分源自锡林郭勒马母本的基础遗传背景，属于共祖遗传（Shared ancestry）而非独立的基因流事件。由于D统计量检测的是相对于外群的衍生等位基因过量共享，对这类受体群体固有遗传背景的保留并不敏感，加之多代级进杂交使蒙古马来源的基因组片段逐渐碎片化，D值接近零是级进杂交育种体系下的理论预期结果。

本研究通过渗入-选择共定位整合框架，将基因流方向信息与多方法选择信号联动分析，揭示了太阳花马杂交育种过程中纯血马渗入与蒙古马祖源保留两类适应性区域的功能分化。Adaptive（-）区域（841个窗口）的数量约为Adaptive（+）区域（168个窗口）的5倍，这一不对称可能与太阳花马高纯血马祖源背景下D值分布向正值偏移有关：D>Mean+2SD的高渗入极端尾部窗口本身较少，而 D<Mean-2SD的蒙古马背景保留区域在基因组中分布更广泛；此外，蒙古马祖源保留区域承受的选择压力涵盖步态协调、神经肌肉功能、能量代谢等多个功能维度，受选择区域的总面积相应较大。

Adaptive（-）区域中，值得关注的发现是马步态调控基因DMRT3。Andersson等^［25］已在冰岛马中证实DMRT3的gait keeper突变（Ser301STOP）是马步态模式的决定性调控因子，该发现已在多个品种中被验证。本研究独立从适应性渗入角度识别到该基因，提示太阳花马在杂交育种中保留了蒙古马来源DMRT3单倍型并受到选择，但本研究未对DMRT3具体变异位点进行基因型分型和表型关联验证，有待后续功能试验证实。同区域内FOXP2（运动协调与学习能力^［27］）和PRKN（线粒体质量控制）等基因的共同富集，进一步提示蒙古马祖源保留区域在耐力表现和草原环境适应中可能发挥协同功能。Adaptive（+）区域中，免疫相关通路富集较为突出（IL4、IL13、IL5等），同时涵盖ACACA（脂肪酸代谢）和ADCY1（cAMP信号通路）等与运动能力相关的基因，McGivney等^［28］的运动转录组研究亦表明运动适应与免疫调节之间存在功能关联，免疫通路的富集可能反映纯血马在长期高密度饲养和竞技运动条件下免疫调控基因的适应性进化被太阳花马继承，同时也可能与其在草原环境中维持耐寒性和疾病抵抗力有关。

综合来看，Adaptive（+）与Adaptive（-）区域在功能上形成互补，前者主要贡献免疫调节和速度潜力，后者主要维持步态协调和草原适应性，共同揭示了太阳花马速度与耐力兼备这一表型特征的分子遗传基础。需要指出的是，当前整合框架主要依赖D值方向与选择信号的空间共定位推断渗入来源，尚无法完全排除部分Adaptive（+）信号可能由太阳花马自身在杂交后独立演化所产生，未来可通过local ancestry推断方法（如RFMix）进一步精确区分。此外，本研究尚缺乏体尺表型数据和运动机能测定数据，无法直接建立基因型-表型关联，上述候选基因与运动性能、耐寒性等表型性状的关系有待后续表型关联分析验证。

与Ding等^［10-13］围绕太阳花马开展的系列研究相比，本研究在遗传结构方向上与其报道一致（太阳花马整体更接近纯血马），但在以下3个方面提供了显著的研究增量： 1）首次采用f₃统计量、D统计量和窗口化D统计量对由纯血马向太阳花马方向的基因流进行系统定量评估，揭示了渗入信号在全基因组范围内的弥散分布模式；2）构建了包含5种选择分析方法的多方法交叉验证框架（n_methods≥3）较既往依赖单一或有限方法组合的选择扫描策略更为稳健；3）提出渗入-选择共定位的整合分析策略，首次在太阳花马中系统鉴定Adaptive（+）和Adaptive（-）区域，将基因流方向信息与选择信号联动，填补了太阳花马在基因流与适应性渗入交叉分析方面的研究空白。该分析框架可推广至其他导血育种物种，为杂交遗传效应解析提供可量化的方法范式。

本研究存在以下局限性：第一，D统计量的稳健性依赖于外群选择，本研究以驴为外群，虽为马属研究中的常用策略^［18］，但仍需在更广泛的外群（如普氏野马）中验证结果的稳健性。第二，f₃统计量的混合检验灵敏度受混合后遗传漂变强度和源群体代理准确性的双重制约，太阳花马经多代级进杂交后漂变积累较强，且以蒙古马代替锡林郭勒马作为母系祖源，均可能削弱f₃对混合信号的检测能力。第三，样本量相对有限（TYH=31，MON=31，TB=32），可能影响低频变异和弱选择信号的检测灵敏度。第四，缺乏体尺表型数据和个体水平的步态及运动机能测定数据，候选基因的功能推断仅基于富集分析和文献推导。第五，top 5%经验阈值的选择具有一定主观性，虽通过n_methods≥3的交叉验证加以约束，但无法完全排除因方法间统计量相关性导致的假阳性富集；此外，DAVID数据库的注释覆盖度对马基因组尚不完善，部分基因可能因缺乏功能注释而未被纳入富集分析。第六，研究未考虑环境因素（如海拔、温度、饲养管理方式等）对选择信号的潜在影响。第七，当前分析框架尚无法完全区分纯血马渗入后的适应性选择与太阳花马独立的本地选择。未来研究应重点从以下方向推进：补充太阳花马群体的体尺（体高、体长、胸围、管围等）及运动机能表型数据，开展全基因组关联分析（GWAS）直接验证候选区域与表型性状的关联；利用local ancestry推断方法（如RFMix）在个体水平解析基因组片段的祖源来源，精确区分渗入选择与本地选择；对DMRT3、IL4和ACACA等关键候选基因开展功能验证试验；结合育种记录和谱系信息，为太阳花马品系的遗传资源评价和育种策略优化提供更精准的分子依据。

4　结 论

太阳花马与纯血马遗传分化（Fst=0.040 7）显著低于其与蒙古马（Fst=0.068 5），ADMIXTURE分析显示纯血马祖源成分均值达76.5%。D统计量证实存在由纯血马向太阳花马方向的显著基因流（D=0.221 6，Z=29.1），渗入信号呈全基因组弥散分布，与以纯血马为父本的级进杂交育种模式一致。

通过渗入-选择共定位整合框架，识别出推定受选择的纯血马渗入区域168个窗口（126个基因），功能富集涉及免疫调节（IL-17信号通路、Fc ε RI信号通路）和神经信号通路，提示纯血马渗入并受选择的基因片段可能与免疫调控和神经信号传导功能相关。

推定受选择的蒙古马祖源保留区域841个窗口（441个基因），富集运动步态（成体运动行为、微管运动）与神经肌肉功能，其中包含马步态调控基因DMRT3，提示太阳花马保留的蒙古马遗传背景可能在运动步态协调和草原环境适应中发挥作用。

Adaptive（+）与Adaptive（-）区域在功能上形成互补，前者主要贡献免疫调节和速度潜力，后者主要维持步态协调和草原适应性，共同揭示了太阳花马速度与耐力兼备这一表型特征的分子遗传基础。

参考文献

[1] Bower M A， McGivney B A， Campana M G， Gu J， Andersson L S， Barrett E， Davis C R， Mikko S， Stock F， Voronkova V， Bradley D G， Fahey A G， Lindgren G， MacHugh D E， Sulimova G， Hill E W. The genetic origin and history of speed in the Thoroughbred racehorse［J］. Nature Communications， 2012， 3： 643

[2] Hill E W， Gu J， Eivers S S， Fonseca R G， McGivney B A， Govindarajan P， Orr N， Katz L M， MacHugh D. A sequence polymorphism in MSTN predicts sprinting ability and racing stamina in thoroughbred horses［J］. PLoS One， 2010， 5（1）： e8645

[3] Sandagdorj B， Baigalmaa T， Burenjargal S， Horiuchi M， Enkhdalai M， Bayanbat D， Janchiv D， Jamyandorj O， Ulzii-Orshikh P， Nyam-Osor P. Relationship between lactational performance and metabolic parameters of Mongolian native grazing mares［J］. Journal of Equine Science， 2021， 32（3）： 91-98

[4] 赵启南， 芒来， 白东义， 赵一萍， 李蓓， 乌尼尔夫， 乌英嘎， 苏日嘎， 敖乳嘎. 蒙古马高负荷运动训练前后转录组差异表达分析［J］. 畜牧兽医学报， 2017， 48（6）： 1007-1016

[5] Zhao Q N， Mang L， Bai D Y， Zhao Y P， Li B， Unerhu， Wu Y G， Su R G， Ao R G. Mongolian horses transcriptome differential expression analysis before and after a high load exercise training［J］. Acta Veterinaria et Zootechnica Sinica， 2017， 48（6）： 1007-1016 （in Chinese）

[6] Han H G， Randhawa I A S， MacHugh D E， McGivney B A， Katz L M， Dugarjaviin M， Hill E W. Selection signatures for local and regional adaptation in Chinese Mongolian horse breeds reveal candidate genes for hoof health［J］. BMC Genomics， 2023， 24（1）： 35

[7] 潘静， 吴凯峰， 周欢敏， 张焱如. 基于全基因组测序的蒙古马运动相关的正向选择基因筛选［J］. 畜牧兽医学报， 2019， 50（10）： 1985-1996

[8] Pan J， Wu K F， Zhou H M， Zhang Y R. Athletics-related positive selective genes screened in Mongolian horse by whole-genome sequencing ［J］. Acta Veterinaria et Zootechnica Sinica， 2019， 50（10）： 1985-1996 （in Chinese）

[9] 芒来， 杨虹. 蒙古马遗传多样性研究进展［J］. 遗传， 2008， 30（3）： 269-276

[10] Mang L， Yang H. Progress in the study of genetic diversity of Mongolian horse［J］. Hereditas， 2008， 30（3）： 269-276 （in Chinese）

[11] 丁文淇， 宫文典， 图格琴， 史晓渊， 李政， 白东义. 草原纯血马新品种培育概述［J］. 中国畜禽种业， 2025， 21（6）： 45-51，3

[12] Ding W Q， Gong W D， Tu G Q， Shi X Y， Li Z， Bai D Y. Overview of new breed cultivation of Grassland-Thoroughbred horses［J］. Chinese Journal of Animal Husbandry and Veterinary Medicine，2025， 21（6）： 45-51，3 （in Chinese）

[13] 齐伯益. 锡林郭勒盟畜牧志［M］. 呼和浩特： 内蒙古人民出版社， 2002

[14] Qi B Y. Animal Husbandry Records of Xilin Gol League［M］. Hohhot： Inner Mongolian People’s Publishing House， 2002 （in Chinese）

[15] Ding W Q， Gong W D， Bou T， Shi L， Lin Y N， Shi X Y， Li Z， Wu H Z， Dugarjaviin M， Bai D Y， Zhao Y P. Preliminary study on the genetic structure and functional candidate genes of grassland-thoroughbreds based on whole-genome resequencing ［J］. Animals， 2025， 15（10）： 1462

[16] Ding W Q， Gong W D， Bou T， Shi L， Lin Y N， Shi X Y， Li Z， Wu H Z， Dugarjaviin M， Bai D Y. Analysis of ROH characteristics across generations in grassland-thoroughbred horses and identification of loci associated with athletic traits［J］. Animals， 2025， 15（14）： 2068

[17] Ding W Q， Gong W D， Bou T， Shi L， Lin Y N， Shi X Y， Li Z， Wu H Z， Dugarjaviin M， Bai D Y. Whole-genome resequencing analysis of athletic traits in grassland-thoroughbred ［J］. Animals， 2025， 15（15）： 2323

[18] Ding W Q， Gong W D， Bou T， Shi L， Lin Y N， Shi X Y， Li Z， Wu H Z， Dugarjaviin M， Bai D Y. Whole-genome resequencing analysis of copy number variations associated with athletic performance in grassland-thoroughbred ［J］. Animals， 2025， 15（10）： 1458

[19] Purcell S， Neale B， Todd-Brown K， Thomas L， Ferreira M A R， Bender D， Maller J， Sklar P， de Bakker P I W， Daly M J， Sham P C. PLINK： A tool set for whole-genome association and population-based linkage analyses ［J］. The American Journal of Human Genetics， 2007， 81（3）： 559-575

[20] Alexander D H， Novembre J， Lange K. Fast model-based estimation of ancestry in unrelated individuals［J］. Genome Research， 2009， 19（9）： 1655-1664

[21] Danecek P， Auton A， Abecasis G， Albers C A， Banks E， DePristo M A， Handsaker R E， Lunter G， Marth G T， Sherry S T， McVean G， Durbin R. The variant call format and VCFtools［J］. Bioinformatics， 2011， 27（15）： 2156-2158

[22] Weir B S， Cockerham C C. Estimating f-statistics for the analysis of population structure［J］. Evolution， 1984， 38（6）： 1358-1370

[23] Patterson N， Moorjani P， Luo Y， Mallick S， Rohland N， Zhan Y P， Genschoreck T， Webster T， Reich D. Ancient admixture in human history［J］. Genetics， 2012， 192（3）： 1065-1093

[24] Chen H， Patterson N， Reich D. Population differentiation as a test for selective sweeps［J］. Genome Research， 2010， 20（3）： 393-402

[25] Browning B L， Browning S R. Genotype imputation with millions of reference samples［J］. American Journal of Human Genetics， 2016， 98（1）： 116-126

[26] Szpiech Z A， Hernandez R D. Selscan： An efficient multithreaded program to perform EHH-based scans for positive selection ［J］. Molecular Biology and Evolution， 2014， 31（10）： 2824-2827

[27] Huang D W， Sherman B T， Lempicki R A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources［J］. Nature Protocols， 2009， 4（1）： 44-57

[28] Racimo F， Sankararaman S， Nielsen R， Huerta-Sánchez E. Evidence for archaic adaptive introgression in humans［J］. Nature Reviews Genetics， 2015， 16（6）： 359-371

[29] Gu J J， Orr N， Park S D， Katz L M， Sulimova G， MacHugh D E， Hill E W. A genome scan for positive selection in thoroughbred horses［J］. PLoS One， 2009， 4（6）： e5767

[30] Andersson L S， Larhammar M， Memic F， Wootz H， Schwochow D， Rubin C J， Patra K， Arnason T， Wellbring L， Hjälm G， Imsland F， Petersen J L， McCue M E， Mickelson J R， Cothran G， Ahituv N， Roepstorff L， Mikko S， Vallstedt A， Lindgren G， Andersson L， Kullander K. Mutations in DMRT3 affect locomotion in horses and spinal circuit function in mice［J］. Nature， 2012， 488（7413）： 642-646

[31] Petersen J L， Mickelson J R， Cothran E G， Andersson L S， Axelsson J， Bailey E， Bannasch D， Binns M M， Borges A S， Brama P， da Câmara Machado A， Distl O， Felicetti M， Fox-Clipsham L， Graves K T， Guérin G， Haase B， Hasegawa T， Hemmann K， Hill E W， Leeb T， Lindgren G， Lohi H， Lopes M S， McGivney B A， Mikko S， Orr N， Penedo M C T， Piercy R J， Raekallio M， Rieder S， Røed K H， Silvestrelli M， Swinburne J， Tozaki T， Vaudin M， Wade C M， McCue M E. Genetic diversity in the modern horse illustrated from genome-wide SNP data［J］. PLoS One， 2013， 8（1）： e54997

[32] Enard W， Przeworski M， Fisher S E， Lai C S L， Wiebe V， Kitano T， Monaco A P， Pääbo S. Molecular evolution of FOXP2 a gene involved in speech and language［J］. Nature， 2002， 418（6900）： 869-872

[33] McGivney B A， McGettigan P A， Browne J A， Evans A C O， Fonseca R G， Loftus B J， Lohan A， MacHugh D E， Murphy B A， Katz L M， Hill E W. Characterization of the equine skeletal muscle transcriptome identifies novel functional responses to exercise training［J］. BMC Genomics， 2010， 11： 398

基金资助

新疆维吾尔自治区重大科技专项（自治区级重点研发计划）(2023A02005-1)