比利时蓝牛的生长特点：在欧洲牛和非洲牛的全基因组对比中

文|爱叨叨的小胖

编辑 | 爱叨叨的小胖

非洲牛栖息在超过5个不同的农业生态区，总的来说，瘤牛在干旱和半干旱的北部萨赫勒-苏丹地区以及包括高地在内的大陆东部很常见；而牛磺酸牛构成了西非半潮湿地区的大多数牛群，那里严重感染了非洲锥虫的媒介：采采蝇。

桑伽牛主要分布在美国西部地区、非洲中部围绕五大湖地区和大陆的南部（图1a）。非洲牛群受到强烈的环境压力，包括炎热、干燥或潮湿的热带气候条件以及沉重和多样的疾病挑战。

因此，它们有望表现出独特的适应特征。居住在采采动物感染地区的恩达马品种和其他西非牛磺酸品种的锥虫耐受性特征就证明了这一点，非洲牛也受到人类对毛色和角大小等特征的影响。

图1 非洲牛种群的地理位置。b核苷酸基因组多样性

尽管它们的生产力远低于集约化生产系统下的商业品种，但本土牛往往是非洲农牧系统中数百万农民的唯一选择，而外来改良品种在传统管理系统中表现不佳。

非洲牛种群，在农业生态区域和生产系统中的地理分布不断扩展，为理解自然选择和人工选择在塑造反刍动物物种功能多样性中的作用，提供了一种独特的遗传资源，此外揭示它们的基因组多样性可能，为了解它们适应各种农业生态系统的遗传机制提供新的见解。

在这篇文章中，我们首次报道了5个非洲本土牛品种的基因组特征，它们代表了大陆牛的多样性：恩达马，属于西非牛磺酸群，对多种传染病具有耐受性；非洲瘤牛和牛磺酸牛种群之间的中间杂交，有大而独特的角和毛色由人类选择；博兰和基纳纳，两个东非瘤牛，分别具有牛和奶制品特性；欧加登，一种生活在炎热干燥环境中的东非瘤牛。

对三种欧洲和1种非洲商业牛品种的全基因组比较分析，使我们能够确定非洲牛品种对热带挑战的独特基因组反应。现在我们准备展开本次试验程序。

我们从10只安科尔、10只博兰、9只基纳纳、10只恩达马和9头欧加登牛中，采集了全血样本（10毫升），我们使用将这些生成了对端读取从全血中分离DNA，结合我们之前发表的53个商业品种样本的数据，进行了每个碱基的序列的质量检查。

最后的读取量对参考序列的总体比对率为98.84%，平均读取深度为10.8×，在整个样本中，reads平均覆盖了参考UMD3.1基因组的98.56%。

获得了样本和DNA重测序数据后，我们将利用SnpSift来鉴定品种特异性富集snp。为此我们使用SnpSift 进行了富集分析，以鉴定杂交特异性snp，在SnpSift中，实施了一些统计检验，如Fisher精确检验和科克伦-阿米蒂奇趋势检验，用于分析由两个因素组成的基因型计数数据。

一般来说，其中一个因素是固定的遗传模型，可以是显性、隐性或共显性，另一种是品种信息，本研究用于鉴定品种特异性富集snp，我们有9个不同的牛品种，用于育种和其他品种，根据这两个因素，我们构建了2个×2或2个×3的列联表。

因此，为每个品种和snp生成了三个不同的列联表，我们分别对2×2和2×3列联表进行了Fisher精确检验和科克伦-阿米蒂奇趋势检验，该测试共使用了37,460,739个SNPs；这就导致了多个测试问题。

为了纠正多重测试误差，我们采用了Bonferroni校正，这是最保守的方法，在发现显著的品种特异性富集SNP后，我们使用snpEff对每个SNP进行了注释，并且还重点关注了这个注释的非同义SNPs。

接下来开始探索非洲牛的选择性扫描区域的统计数据。为了揭示每个品种的本地适应的遗传变异，我们在种群之间进行了比较： (1)恩达马与所有其他非洲，(2) Ankole与所有其他非洲；(3)所有非洲与所有商业；(4)所有B.indgus与所有商业。

然后，我们将基因组分割成50 kb的非重叠片段，使用一个窗口内所有snp的最大XP-EHH得分作为每个窗口的汇总统计量，为了考虑SNP密度的变化，我们根据以500个SNPs为增量的SNPs数量将基因组窗口分类。

在每个bin中，对于每个窗口i，统计值大于在i中的窗口的分数被定义为经验P值，遵循之前报告的方法，P值小于0.005（0.5%）的区域被认为是感兴趣的品种组中的显著信号，这种方法是合适的，特别是对于不可靠的人口统计模型和参数，如牛的情况。

然而使用这种方法检测到的选择位点可能不足以代表站立变异的选择。

此外，我们还使用XP-CLR 来寻找基因组中，在该位点上的等位基因频率的变化发生得太快的区域，这是通过受影响区域的大小来评估的。

得知了等位基因频率的变化影响后，现在我们将开始研究种群分化与结构。对于本次研究最重要的PCA，我们使用全基因组复杂性状分析来估计特征向量，它与在特征链中实现的PCA的特征向量渐近等价，包含了来自所有样本的基因型数据，对于混合分析，我们使用plink将基因型数据限制在~ 20,000个总snp的随机子集上，以运行结构版本2.3 中的“混合”模型。

在经过5万次迭代之后，我们选择了2万次迭代并对每个K值重复分析10次，然后在结构收割机中分析结构的输出，该收割机实现了Evanno方法来推断最有可能的聚集种群数量，我们使用VCFtools 4.0 来估计核苷酸多样性和Fst差异统计量，使用VCFtools实现的Fst和加权Fst估计量。

表1 在每个品种的比较中，从XP-EHH和XP-CLR中确定的主要候选区域的汇总

基于通过PLINK计算出的所有样本的IBS距离矩阵数据，利用FigTree v1.4.0构建了邻域连接树。我们还利用TreeMix中实现的最大似然（ML）方法推断了种群水平的系统发育，使用窗口大小为1000来解释连锁不平衡（-k）和“-global”来生成ML树。

迁移事件（-m）已按顺序添加到该树中，利用PopSizeABC，设置牛种群(突变率和重组为1e-8，MAF >为20%，每个片段突变大小=为2000000数据集，推断最近的牛种群统计历史。

为了等到最优结果，我们生成48头非洲土著牛的个体基因组，分别覆盖~11 X，并与公开的商业牛品种基因组进行基因分型，这些品种包括印度牛、非洲牛、欧洲-亚洲牛和桑，共生成了65亿条reads或~644 Gbp的序列。

使用Bowtie 2，将reads与牛磺酸参考基因组序列UMD 3.1进行比对，平均比对率为98.84%，覆盖了参考基因组的98.56%，与之前对瘤牛的分析一致，非洲印度样本与参考基因组UMD 3.1的总体比对率与非洲牛磺酸样本的比对率相当。

在过滤了潜在的PCR重复序列并纠正了由于INDELs的存在而导致的错位后，我们使用GATK 3.1检测了单核苷酸多态性（SNPs），在进一步的分析中使用候选snp之前，我们应用了几个过滤步骤来最小化假阳性调用的数量。

特别是，根据以下标准去除snp：phred标度质量评分、映射质量、质量深度和phred标度P值。最终总共保留了3700万个~的SNPs，并使用SnpSift鉴定出了品种特异性的SNPs（详情请见上文图1b）。

来自45个非洲样本的基因组DNA此外，使用BovineSNP50基因分型珠芯片进行基因分型，以评估重测序数据中SNP调用的准确性，我们观察到BovineSNP50基因分型珠芯片SNP与整个样本的重测序结果之间的~95%的总体基因型一致性，为SNP调用的准确性提供了信心。

上文中的图1b显示了每个品种中存在的snp的数量，包括特定品种的snp，从不同的牛谱系来看，在瘤牛，那里绝大多数的SNPs是纯合的三个品种代表候选非洲瘤牛谱系特异性变异纯合的，大部分（65.13%）存在于基因间区。

其余的SNPs分别位于开放阅读框的上游（3.90%）和下游（3.96%）、内含子（26.0%）和非翻译区（UTRs，0.240%），外显子包含总snp的0.69%，有115,439个错义突变和1336个无义突变。

核苷酸多样性衡量，是一个群体内的多态性程度，它被定义为从样本群体中随机选择的任意两个DNA序列之间每个位点的核苷酸差异的平均数量，在10 Mb的全基因组窗口尺度上，与所有非洲本土品种相比，欧洲商业品种的核苷酸多样性水平都有所降低（图2d）。

图2 非洲牛与商业牛之间的人口结构和关系

在这里，全基因组水平上的核苷酸多样性水平的降低是预期的，这可能是世代间的密集人工选择和/或遗传漂变的结果，随后是以低有效种群规模为特征的人口统计历史。有趣的是，恩达马也显示出相对较低的遗传多样性，这可能是最初的低有效种群规模和/或疾病挑战后的种群瓶颈的遗留问题。

非洲瘤牛和安科尔桑的核苷酸多样性最高，这些是混合牛磺酸×瘤牛品种，具有相对较大的有效种群规模，与其他商业品种相比，商业汉宇相对较高的核苷酸多样性可能反映了较弱、有针对性和较短的选择历史。

我们使用特征链对常染色体snp基因型数据进行了主成分分析（详情请见上文图2a），该分析忽略了品种成员，但仍然揭示了明确的品种结构，因为样本来自同一品种集群在一起。前两个pc，分别解释了总变异的16.0%和3.4%，将非洲品种与非非洲品种分开，其中Ankole牛处于中间位置。

分别基于非洲、商业和牛磺酸样本的PCA显示，没有品种之间的混合或品种内存在异常动物的证据。

为了进一步了解群体中的混合程度，我们对一个随机抽样的SNPs子集（~20,000个SNPs）使用了结构。我们将K从1增加到9，其中K是假定的祖先种群数量。

分析表明，K = 2是我们的样本中最可能的遗传不同群体数量，反映了牛磺酸和瘤牛在牛种群中的差异，在K = 3中，Ankole显示出明显的遗传异质性，与非洲、亚洲瘤牛和商业牛磺酸遗传背景具有共同的基因组祖先。

K值的增加表明，商业种群的品种同质性水平高于非洲瘤牛品种，此外一个邻居连接树（详情请看上列图2c）将每个品种分离成它自己的独立的分支，欧洲品种聚集在一起，然后是汉宇和恩达马，同样，所有的非洲瘤牛品种聚集在一起，安科尔动物位于瘤牛和恩达马之间的中间位置。

与其他非洲牛相比，恩达马牛似乎遭受了更严重的人口下降，这一观察结果与在热带半湿润和潮湿的西非环境中，祖先种群到达和适应后的初始种群瓶颈相一致，这些西非牛种群最近受到了新的环境压力，施加了强大的适应性限制。

此外，奥加登和基纳纳的估计显示，大约1000年前的牛口规模略有增加，这与第一波瘤牛通过大陆之角到达的时间相对应，它们都有一个共同的种群下降，从大约10,000 BP开始，这可能是新石器时代驯化事件的结果。

然后，我们利用Treemix重建了9个品种的最大似然树（图3b）和残差矩阵，以解决种群历史关系，并识别出相互相关的种群对被这棵树捕获，将顺序迁移事件添加到树中，我们发现一个推断的迁移边产生一个残差最小的树，因此它是最适合数据的。

图3 非洲牛的有效种群规模和历史

我们观察到一个具有统计学意义的迁移边缘（P < 2.2E-308），估计权重为11.4%，近年来，Ankole牛越来越多地与牛磺酸品种杂交，包括50年前首次引入乌干达的荷斯坦牛。

我们比较了非洲牛品种的基因组与SNP芯片数据，多样性高估了牛磺酸血统和低估指示血统全基因组测序，可以克服这个限制确定偏差，适当使人口分析人口和识别目标在非洲b的选择。

特别是，我们使用跨群体扩展单倍型纯合度（XPEHH）和跨群体复合似然比（XP-CLR），检测了扩展连锁区域上的极端单倍型纯合度和等位基因频率分化，考虑到非洲籼牛之间的遗传距离较近，将恩达马牛和安科尔牛品种分别与所有其他非洲品种进行比较，以鉴定非洲品种特异性特征。

XP-EHH在小样本量（低至10个样本）下保持功率，此外当群体间之间的遗传距离（FST）估计大于或接近0.05时，如我们的分析，每个群体少于20个个体就足以用于群体分化分析。

为了能够比较群体间的基因组区域，我们将基因组划分为50 Kb的不重叠片段。离群区域被认为是繁殖特定的候选区域，用于进一步分析（单倍型和多态性）。

我们首先研究了，对锥虫挑战的耐受性对非洲牛的基因组的影响。非洲锥虫是一种细胞外原生动物寄生虫，可导致人类（昏睡病）和家畜的严重疾病；大约有6000万人和5000万头牛生活在锥虫感染的风险中。

在少数“耐锥虫”的本土非洲牛品种中，西非的恩达马是最典型的，而“新来的”印度通常对锥虫病高度敏感。因此，我们将恩达马基因组，与所有其他非洲牛品种进行了比较。XP-EHH和XP-CLR分析的离群窗口分别包括124和106个基因，其中28个是两种分析中共同的。

这种相对适度的重叠，可能是由于设计用于检测受完全或不完全选择性扫描影响的区域测试之间的功率差异造成的。

其中，我们发现HCRTR1（XP-CLR = 597.3）编码下下缀素受体A（图4），属于偶联受体超家族中的I类亚家族，与Ca2 动员偶联。下丘脑分泌素是由下丘脑外侧和穹窿周围区域的一小群神经元产生的，它们参与控制哺乳动物的进食行为。

图4 在N‘Dama HCRTR1、SLC40A1、EPB42和STOM基因区域的选择性扫描特征

与其他非洲牛相比，N‘Dama在HCRTR1区域表现出几乎纯的单倍型纯合度，我们还在该基因中检测到7个非同义变异（详情请见上列图4b）。大量研究表明，下丘脑泌素基因的多态性与进食和饮酒行为的改变有关。

特别是，食欲素-A、G蛋白偶联受体的内源性配体、刺激的食物消耗和食欲素信使RNA通过禁食而上调，这些独立的研究表明，下丘脑泌素在进食调节中起着重要作用。这可能解释了恩达马在锥虫感染后，在维持体重和抵抗酗酒和消瘦方面的优越能力。

恩达马牛在耐锥虫，至少有两个额外的特征：抵抗贫血的能力和控制寄生虫增殖的能力。贫血是锥虫感染最突出和最一致的临床体征，是治疗的主要指标。我们在基因组区域内发现了5个与贫血相关的阳性选择的基因（SLC40A1、STOM、SBDS、EPB42和RPS26）。

铁出口商SLC40A1（XP-EHH = 3.32，XP-CLR = 831.1）对铁稳态是必不可少的，因此它与缺铁性贫血有关。该基因显示出核苷酸多样性的局部减少，单倍型模式有所扩展（详情请看上列图4c）。

值得注意的是，我们在N‘Dama中发现了一个固定的SLC40A1单倍型，在其他非洲牛和商业品种中发现的频率分别为24%和58%，强烈支持该基因的选择（详情请看上列图4d与4e）。

胃蛋白（STOM，XPCLR = 525.0）是一种以罕见的人类溶血性贫血命名的基因，它编码一个31-kDa的完整膜蛋白。SBDS（XP-EHH = 2.91）EPB42（XPCLR = 511.1）基因的突变分别导致低色性贫血和遗传性溶血性贫血，而在钻石-黑扇贫血患者中发现了RPS26（XP-CLR = 562.8）基因的突变。

我们进一步筛选了这些候选基因的非同义突变代表的假定的功能变异。值得注意的是，错义snp改变了STOM中的氨基酸(p。Met48Val)和EPB42 (p。Arg503His)蛋白质。与所有其他品种相比，这两种等位基因变异在恩达马牛中都是完全固定的（详情请看上列图4f和4g）。

非洲牛品种已经进化到，能够适应撒哈拉以南非洲普遍存在的恶劣环境条件，如热带牲畜疾病、高太阳辐射和温度、干旱和营养不良的，这些环境条件在撒哈拉以南非洲地区普遍存在，积极选择的信号可能在非洲品种中普遍存在。

为了调查这一点，我们将所有非洲品种合并，并与商业品种进行比较，以确定常见和独特的非洲基因组特异性选择特征，在这项比较中，XPCLR和XP-EHH分析显示了包含252个基因的异常窗口。

其中，我们发现了包含牛淋巴细胞抗原基因的区域。通过对该区域的详细检查，我们确定了6个BOLA单倍型块，其中主要的非洲牛单倍型对应于商业牛的对比或次要单倍型。

BOLA-DRB3等位基因与牛对蜱侵染的抗性有关，牛淋巴细胞抗原复合物由于其在宿主免疫中的重要性，在过去30年来得到了广泛的研究。

现在我们开始研究非洲牛的耐热性。为了鉴定非洲牛中负责热调节的基因组区域，我们利用先前鉴定的13个耐热性数量性状位点（QTL）区域和18个热休克蛋白，选择了一个先验候选基因，这些区域都没有得到我们所使用的XP-EHH和XP-CLR的共同指标的支持。

然后，我们分析了非洲牛与欧洲和亚洲牛磺酸的单倍型纯合子模式。与我们之前的结果一致，我们发现在检测随机基因组区域时，商业品种的单倍型共享更为广泛。

然而，看看非洲品种的候选地区与商业相比，牛的耐热基因（图5a）和一个热休克蛋白，热休克70 kDa蛋白4，表明该地区的耐热性的选择性扫描，细胞对热应激的耐受性是由一个热休克蛋白家族介导的。热休克蛋白70可以促进细胞抗热损伤的保护和防止蛋白变性。

在这些地方的单倍型共享的程度有两个地区，更为广泛非洲牛则在恩达马，这是与之前的一份报告是一致的。

图5 一种与非洲牛的耐热性相关的选择性扫描

我们的氧化石墨烯分析显示了Wnt信号通路（P00057）最显著的富集，以及参与调节皮肤血流的通路：内皮素信号通路（P00019）和组胺H1受体介导的信号通路（P04385）。

在热稳态挑战时，皮肤血流的体温调节控制对于维持正常体温至关重要，特别是在体温过程中皮肤血流的上升含有H1组胺受体成分，这些途径可能在非洲牛中迅速进化，这可能解释了它们在细胞和生理水平上与温带牛品种相比具有完全不同的耐热性程度。

我们的研究结果在全基因组水平上揭示了非洲牛独特的适应性多样性，同时强调了在非洲大陆可持续提高牲畜生产力的机会。