文|爱叨叨的小胖

编辑 | 爱叨叨的小胖

背景

在牛中，牛奶蛋白的组成主要是由遗传因素提供的这很有趣，因为它决定了奶酪的制作而通过分析中红外（MIR）光谱，可以大规模地预测牛乳蛋白的组成，这是一种常规的方法。

在之前的一项基于牛50 K单核苷酸多态性（SNP）阵列的全基因组关联研究（GWAS）中，法国三个主要品种(蒙特贝利亚德（MON）和诺曼德（NOR））的牛奶蛋白组成具有非常显著的影响，然而由于50 K SNP阵列只包含基因组变异总数的一小部分，我们无法直接精确定位候选突变。

在1000头牛基因组，参考群体的第4次试验中，通过分析27个不同品种1147头牛的全基因组序列（288头、28头和24头），构建了包含5600万个SNPs的数据库，但这其中包括了288头的基因缺失（InDel）。

这些数据，可以用来从实验或常规获得的50 K SNP基因型中推断WGS，这样就可以得到大量的动物数据，特别是那些具有表型的动物的推断WGS。

由于WGS包含了几乎所有的基因组变异，它们应该包含一个特定性状的因果突变，因此它们提供了更高的GWAS分辨率。

然而由于奶牛品种中存在长程连锁不平衡，近品种GWAS的解决往往受到限制对于品种间共享的因果突变，可以使用多品种模型来识别包含数量性状位点（QTL）的区域。这种方法利用了在每个品种中发生的历史重组事件，导致了较短距离上的连锁不平衡和更好的分辨率。

在此，我们研究了来自HOL、MON和NOR牛的六种主要乳蛋白含量的GWAS序列水平，即α-乳球蛋白、β-乳球蛋白以及αs1、αs2、β和κ酪蛋白的GWAS结果。

品种内、多品种和条件分析的结果显示，ft具有最显著的变异，除了其他测试变异，一起检查，以此来确定每个基因组区域的潜在候选变异。以下是我们的试验方法。

实验方法

本次试验方法简单地说从法国三个主要奶牛品种蒙特贝利亚德（MON）、诺曼德（NOR）和荷斯坦（HOL）的156,660头奶牛的848,068个牛奶样本中获得了MIR光谱。

用光谱来预测牛奶蛋白质含量（PC）和牛奶蛋白质组成，预测了六种主要乳蛋白（αs1-CN、αs2-CN、β-CN、κ-CN、α-LA和β-LG）在g/100 g蛋白中的Te含量，也分析总酪蛋白含量和总乳清蛋白含量（分别为Σ-CN和Σ-WP）。

为了调整非遗传缺陷的表型，使用基因试剂盒软件，对近种混合模型应用于测试日数据，单性状重复性模型包括遗传、永久环境和残留随机效应，以及畜群×试验日、胎次×哺乳期、产犊年×月和光谱仪×试验月×效应。

我们将该模型应用于frst两个撕裂期的数据，其中，包括在研究期间至少三个跨撕裂期的测试日记录。

10天测试日数据，对模型中包含的所有非遗传缺陷进行了校正并取每头牛的平均值，对于每个性状和每头奶牛一个单一的表型被取消，随后用于GWAS分析共分析了293,780、58,594和72,973个试验日记录，分别对应44959头MON、12428头NOR和14530头HOL牛，平均每头牛6.5、4.7和5.0个试验日记录。

在这些奶牛中，8010头用snp50珠芯片进行了基因分型，我们采用了以下质量控制方法：个体呼叫率，必须高于95%，SNP呼叫率高于90%，次要等位基因频率（MAF）高于5%，基因型频率必须与P>10-4处于哈代-温伯格平衡，Te fnal数据集包括37,332-41,028个SNPs（表1），具体取决于7907头表型奶牛（3032个MON，2659个NOR和2216个HOL）。

利用FImpute软件，将7907头奶牛的Te 50 K SNP基因型，输入到全基因组序列（WGS）中，准确、快速地分析了大数据集。

为了提高结果的准确性，我们采用了两步的方法：1、从50到777 K的高密度（HD）SNPs，然后，从推断的HD SNPs到WGS；2、对所有每个品种都分别进行推断，使用品种特殊（从50 K到HD SNPs）或多品种（从HD SNPs到WGS）参考面板取决于目标密度。

在每个MON、NOR和HOL品种中，使用品种内参考群体进行HD SNP水平的推断，分别包括522头MON、546头NOR和776头的HOL公牛进行了基因分型，在去除质量控制花瓶中失败的SNPs后，每个品种中保留了大约55万个SNPs（详情请见上列表1）。

利用1000个公牛基因组联盟中1144头公牛的WGS变异，从HD SNP基因型中推断WGS变异；这些公牛代表27个牛品种，有288个HOL、28个MON和24个NOR个体，在“1000个公牛基因组”联盟中使用的Te协议被取消。

使用Illumina合成测序技术，对所有个体的全基因组进行2×100 bp的对端测序并进一步筛选序列质量，随后与UMD3.1参考序列进行比对使用SAMtools 0.0.18 检测小的基因组变异（SNPs和InDel）。

原始变异进一步改变，产生27,754,235个常染色体变异。随后，用集成后的变异差异预测器（VEP）管道对过滤后的变异进行注释，并使用SIFT工具预测氨基酸变化的差异。

通过比较推断的基因型与使用定制芯片对同一奶牛的基因型进行重新基因分型，评估从HD到序列的推断精度，这在归责过程中没有使用附加信息。

有两个数据集：(1)168头荷斯坦牛，用EuroG10k Illumina芯片的frst版本进行基因分型，有721个额外标记；(2)2142头蒙贝利亚德奶牛，使用同一EuroG10k芯片的第四个版本，包含3082个额外的SNPs。

通过真实基因型和推断基因型之间的平方相关性，以及基因型和等位基因的一致性率来衡量推断的准确性。

为了去除推断精度最低的snp，只保留MAF大于0.02的变量进行进一步的关联分析。在每个品种分析中，约有1300万变异被纳入多品种分析（详情请见上文表1）。

全基因组序列关联分析

我们对所有多态变异与9个乳蛋白组成性状进行了单性状关联分析： PC、α-LA、β-LG、αs1-CN、αs2-CN、β-CN、κ-CN、Σ-CN和Σ-WP（表2）。

所有的关联分析都使用GCTA软件的mlma选项进行，该软件应用了一个包含候选变量的混合线性模型：

其中y是预先调整的表型向量；u是总体平均值；b是待关联检验的候选变异的加性系数；σ2是编码为0、1或2的推断基因型；u~N（0，Gσ2 u）是随机多基因缺陷的向量，G为基因组关系矩阵（GRM），使用HD SNP基因型计算，σ2为多基因方差，基于零模型（σ2=1u+ u + e）估计。

然后检验每个变异与性状之间的关联；e~N（0，Iσe 2）为随机残差缺陷的向量，I为单位矩阵，σe 2为残差方差。

在品种内测试日记录的数量在奶牛之间没有很大的差异，因此假设残差在奶牛之间是恒定的。对于多品种关联分析，模型(2)在模型(1)中添加一个fxed品种ectv，W作为与品种表型的发病率矩阵，x、b、u和e如前面所述：

TeBenferroni校正应用于阈值，经过严格矫正以后，我们认为所有1300万个测试都是独立的，因此5%的全基因组显著性阈值，对应的名义P值为3.7×10-9（-log10(P)=8.4）。

对于给定品种的给定性状，当两个连续的QTL区域有重叠的信任区间时，或当信任区间的界限之间的距离小于1 Mbp时，只保留呈现最显著结果的信任区间。

实验结果

GWAS分析中使用的SNPs，在序列水平上的推断精度的测试结果（MAF≥2%），在蒙特贝利尔德和荷斯坦品种中，验证集的推断基因型和真实基因型之间的平方相关性，分别达到76和84%，该表中还给出了一致性率的总体结果，图1为两个品种根据MAF进行的归算精度。

在1300万个测试变异中，有71,755个具有全基因组显著缺陷（-log10(P)≥8.4）。

对于每个性状，不同品种间显著相关变异的数量是相对一致的（表3）αs2-CN、β-CN、αs1-CN和PC的含量较低（从193到2394）；β-LG、κ-WP和αs2的含量较高（从8716到19952）Σ-CN；以及α-LA和κ-CN的中间体（从4110到8248）。

在这些变异中，0（PC）到2266（β-LG）在三个品种中共享，多品种分析比任何品种内分析都更强大，并且检测到更多具有显著缺陷的不同变异（34248），可是在PC、α-LA、β-LG、Σ-CN和Σ-WP的品种内分析中，每个性状的变异数大于多品种分析（表4），可能是因为品种内LD的存在。

通过合并为不同性状和品种获得的重叠QTL区域，并对相应的显著性结果进行分组，来消除QTL区域。这些区域的可信度区间，如方法部分所述。

在品种内和多品种分析中，鉴定出了34个对一个或几个乳蛋白组成性状具有显著缺陷的QTL区域，这些树位于6、11和14染色体上，几乎所有分析的蛋白质组成性状，都具有显著的多效性效应，而大多数（31 QTL）通常只表达一个性状。

Te 34 QTL分布在29个牛常染色体中的17个上，每个染色体有1~7个QTL，几乎所有（31个）在多品种分析中检测到（表5）.而在MON、NOR和HOL品种内分析中分别发现11、8和11个QTL区域。

在3个品种中检测到4个QTL区BTA6（45.8-46.9Mbp和85.2-87.4Mbp）、BTA11（103.3 Mbp）和BTA20（58.3-58.4Mbp），BTA29上另一个约9.6 Mbp的区域是MON和HOL共同的，BTA14上另一个1.7-1.8Mbp的区域，是HOL和NOR共同的（图2b），2~3个品种之间共有的6个QTL缺陷最显著，仅在BTA2上的NOR品种中，检测到一个QTL，为131.8Mbp。

多品种分析，比品种内分析检测到更多的QTL区域：在多品种分析中，检测到的31个QTL中有14个在品种内分析中没有发现。

对于在品种内，和多品种分析中发现的17个QTL区域，最显著变异的-log10(P)值，几乎总是高于品种内分析；即使是对大多数只在一个品种内分析中有显著缺陷的地区，也是如此。

对于这些QTL，最显著变异的平均-log10(P)值为64，而在MON、NOR和HOL品种内的分析中，则分别为49、46和42，此外，多品种分析产生的QTL可信区间，包含的变异数量少于品种内分析产生的变异数量。

在17个QTL区域中，在多品种分析中平均发现134个变异（2-374个），而在MON、NOR和HOL品种内分析中，分别发现189个（39-335个）、287个（61-872个）和308个（9-1236个），然而在一些QTL区域，特别是位于BTA2（131.8 Mbp）、6（38 Mbp）和19（61 Mbp）的区域，品种内分析中的显著变异数量小于多品种分析。现在我们将展开本次研究的讨论。

讨论

在实验中，我们从MIR光谱中预测牛奶蛋白组成的全基因组序列扫描结果，并利用奶牛的WGS，对来自三个法国的奶牛品种，共计7907头进行了品种内和多品种分析。

利用推断的WGS，我们能够整合之前从50 K SNP基因型中检测到的39个QTL中的22个并鉴定出12个新的QTL，除了50 K芯片的遗传参数结果和QTL检测结果外，这些结果表明，MIR预测对于遗传研究是非常准确的，重复的测试日记录，补偿了一些蛋白质适度的MIR预测准确性。

在QTL中，发现了50 K SNP基因型数据，没有发现估算WGS，可能是因为不同的方法在两个研究，也可能因为更严格的意义阈值应用，对于WGS数据上的GWAS，我们使用的非常严格的阈值可能降低了检测能力，但最小化了假阳性QTL的数量。

相反，WGS数据的更好的分辨率，结合多品种GWAS方法的力量，导致检测到了12个之前在50 K SNP研究中没有发现的QTL，评估标记密度对GWAS结果的影响，我们从WGS结果中提取了50 K和HD GWAS结果。

在几个基因组区域，例如BTA2、6和20上的区域（图3），WGS数据分辨率的提高，清楚地使识别两个或多个峰成为可能，而对50 K SNP数据的分析只检测到一个峰。

此外，WGS分辨率允许使用分析多品种的方法，这有望更好地估计罕见变异的缺陷，并减少相邻变体之间的LD，通过多品种分析，我们检测到14个未检测到的QTL。对于在品种内和多品种分析中检测到的QTL，多品种方法比无品种分析，提供了更小的QTL可信间隔。

在我们的研究中，使用的三个法国品种相关性并不强，基于50 K SNP数据，有学者认为在不同品种群体中，牛的遗传多样性具有高度一致性。

而在之前的一项研究中发现，从50 K到HD SNP密度的推断，在所有三个品种中都非常准确，这里使用的HD基因型的数量（>500）进行校准，第二步从HD SNPs到WGS，我们使用了1000头牛基因组联盟的Run 4参考群体，其中包含1147头公牛，其中288头为HOL，28头为MON，24头为NOR。

由于与其他两个品种相比，已测序的荷兰公牛数量更多，因此对HOL数据的推断比对MON数据更准确。

在NOR中，由于两个品种主要祖先的种群结构相似，全基因组序列数量相似，我们预计推断精度与在MON中获得的精度接近。

正如预期的那样，对于具有低MAF的变量计算精度下降了，为了限制归责误差对GWAS结果的影响，在分析中剔除了MAF低于2%的变异，我们在这次研究中提出的几乎所有作为候选变异的遗传变异都具有中度到高的MAF。

将品种内、多品种和条件GWAS分析与功能注释相结合，似乎是共享和品种物种QTL分化的一个很好的策略。

这种方法，还可以直接识别具有非常少量候选变异的基因，甚至在某些情况下，一种独特的变异，这似乎是解释观察到的缺陷的最佳候选变异。

平均而言根据品种的不同，我们在GWAS中检测到的QTL变异中有60-73%位于基因；这大约是全基因组规模上，基因变异百分比的两倍（35%）。

大多数显著变异位于49个不同的基因中，其中22个特别有趣，要么是因为它们在多个品种中被发现，要么与几个性状相关，要么因为它们之前被描述为注入了牛奶成分。

Tese 22基因位于11个不同的基因组区域，具有显著的蛋白-蛋白富集能力，是解释MON、NOR和HOL奶牛乳蛋白组成变化的最合理候选基因。

在四个基因组区域（BTA1、2、11和27）上，我们鉴定了一个独特的候选变异（或LD中的一些候选变异），所有三个品种（分别在SLC37A1、ALPL、PAEP和AGPAT6基因中）共享。

在其他三个基因中，我们认为存在一个品种特异性候选变异、或几个候选致病变异。最后，4个区域，包括BTA14上的DGAT1区域和BTA6上的3个区域，更为复杂，因为它们包含几个候选基因，每个候选基因都有几个候选变异。

SLC37A1（BTA1）和αs1 CN/αLA

TeSLC37A1基因编码一种葡萄糖-6-磷酸转运体在乳腺中高表达，它可能是一个很好的候选基因，来解释在MON和多品种分析中，对αs1-CN和对α-LA表型的QTL的影响。

该基因共有138个不同的变异位于QTL的合并区间内，其中133个为内含子，2个为同义子，3个位于下游。

对于αs1-CN/MON、αs1-CN/multi和α-LA/multi结果，峰中的80、81和74个最显著的变异，分别位于SLC37A1的内含子区。

在MON分析中，检测到一个下游的αs1-CN变异，在显著变异中排名第104位，而多品种分析显示有3个下游变异，排名第81、87和103位。

所有位于峰顶部的内含子变异，都处于强LD中，但只有一个变异（indel）位于144,397,274 bp，属于所有三个前10名名单；它在αs1-CN/MON排名中排名第一，在αs1-CN/多品种排名中排名第九，在α-LA/多品种排名中排名第四。

在αs1-CN/多品种分析中检测到Te前1内含子变异，在144398814bp，在αs1-CN/MON峰中排名第75位，在α-LA/多品种峰中排名第76位。

之前的两项研究，描述了SLC37A1基因变异对产奶性状的影响。在对HD SNP基因型的分析中，国外有位学者描述了，位于该区域144.325~144.525Mbp之间的6个变异；表达最显著的变异位于该基因的一个内含子区（144,414,936 bp）。

在我们的研究中，该变异，被纳入了多品种分析检测到的QTL的可信区间内（-log10(P)=10.2），但它排名第101位。SLC37A1基因中，LD的另外两个内含子变异，分别为144367474和144377960bp，先前被认为是磷浓度和产奶性状变化的最佳候选突变。

然而，在我们的研究中，尽管MAF值相对较高，但这些变异，对所有分析性状的-log10(P)值均低于6。在另一项推断WGS水平后的靶向QTL区域的研究中，缺陷最显著的变异位于144381564bp。Tis变异，与我们在分析中发现的候选变异很接近，但它可以在我们的人群中被排除为因果变异，因为它在这里分析的MON、NOR和HOL个体中是单态的。

我们进行的Te条件分析，包括两个最佳候选变量。Tese显示，在模型中包含位于144,398,814 bp的变体，完全去除了原始信号，而其他两个变体一个较小的显著峰则持久化。Tis变异在αs1-CN和α-LA表型上有不同的缺陷，但在前者上有更明显的缺陷，因此构成了我们研究中检测到的缺陷的最可能的候选变异。

ALPL（BTA2）和αs2 SubarkCN

在BTA2上发现的131.8 Mbp、QTL在几个性状（αs2-CN、β-CN和κ-CN）上具有显著差异。特别是，尽管在所有品种内和多品种分析中都检测到αs2-cn相关峰，即使在MON和HOL分析中，最大-log10(P)值，没有达到我们在本研究中应用的严格的阈值。

在所有的分析中，最显著的变异，都位于ALPL（碱性磷酸酶）基因的内含子区域。在3种品种内分析中，差异最显著的是：NOR131,806,882bp，MON131,850,456bp，HOL序列131,808,301 bp。

相反，在多品种GWAS中，检测到的排名最高的峰值变异位于131,806,882 bp。所有三个单品种条件分析，包括这些变异作为fxed缺陷缺乏峰值。结果表明，3个内含子变异在三个品种中都处于强LD，因果突变可能在品种之间共享。

在所有峰顶部的变异中，131,806,882 bp的内含子变异，是ALPL基因中αs2-CN含量最可能的候选变异；在NOR、MON、HOL和多品种峰中分别排名第1、6、26和1。

AGPAT6（BTA27）和κSubarkCN

Te多品种分析检测QTLκ-CN含量位于约36.2MbpBTA27。多品种分析，四个最重要的变异，位于基因间区域但变异排名5到17日，是位于AGPAT6基因，之前被描述为一个功能基因乳脂肪含量与多效性缺陷在其他牛奶成分，特别是蛋白质含量。

Te fve最重要的变异基因，是在完整的LD和位于上游地区（在36、209、319、36、211、252、36、211、258和36、211、708个基点）或5‘-UTR地区（36、212、352个基点）的AGPAT6基因。在fve连锁变异中，MON的MAF为0.46，NOR为0.47，HOL为0.39（多品种群体为0.44）。

当κ-CN表型以这些突变中任何一个的缺陷为条件时，关联信号在MON、NOR和多品种分析中完全消失。位于上游区域的4个变异之前，是荷斯坦牛和弗莱克维埃牛的候选因果多态性。这些作者指出，在36,211,252 bp处的多态性，是最可能的致病突变，因为它有很高的可能性出现在一个转录结合位点内。

此外，有学者描述了牛奶组成性状（脂肪、蛋白质和乳糖）与AGPAT6基因的10个变异之间的强关联。这10个变异的树，是我们研究中最显著的变异之一，分别位于36,209、319、36、211、708和36,212,352 bp。

我们鉴定了AGPAT6基因的致病变异；其中，36212352bp的变异，似乎是最可能的致病突变，因为它位于AGPAT6基因的5‘-UTR区域。然而，尽管候选变异的MAF较高，但在HOL分析中缺乏显著的缺陷，可能补充了有待解释的其他缺陷。

PICALM（BTA29）和αs1 SubarkCN

Te αs1-CN表型，由一个位于BTA29上约9.5 Mbp的基因组区域注入。在MON、HOL和多品种分析中发现了显著的相关性，在NOR中发现了一个接近显著性的峰值（-log10(P)=7.9）。在MON和HOL分析中，最显著的变异是基因间变异，同样，在多品种分析中，位于信任区间内的所有9个变异都是基因间变异。

最显著的非基因间变异，定位于MON和HOL的PICALM基因中。两个内含子变异(9,651,065，在HOL分析中排名第10位于基因上游区域（9,611,304 bp）。当进行条件GWAS分析时，包含内含子变异删除了MON中的峰值，但在HOL分析中没有，相反包含上游变异删除了HOL中的峰值，但在MON分析中没有。

在NOR中，当内含子或上游变异被排除时，相关的峰值，依旧持续存在。这些结果表明，要么是致病变异的品种是不同的，要么是，几个相关的致病变异解释了在该地区观察到的显著缺陷。Te PICALM基因编码，一种磷脂酰肌醇结合网格蛋白组装蛋白，该基因的多态性与人类阿尔茨海默病的风险相关。然而，到目前为止，还没有报道该基因的多态性与牛奶成分之间的联系

结论

我们的研究提供了证据，表明基于GWAS的方法应用于规模表型、全基因组序列和多个品种，提供了足够的分辨率来识别候选基因，并直接确定这些基因中有限数量的候选变异。

一些在品种间共有的一些变异，被确定为法国三个主要奶牛品种中，牛奶蛋白组成变化的可信候选变异。

Tey既位于先前发现的影响牛奶成分的基因，也位于没有这种关系的基因，在未来功能分析将使建立这些候选变异和牛奶蛋白表型之间的因果联系。

即使在这些研究完成之前，我们的研究结果，就已经有机会通过识别与牛奶成分变化相关的遗传变异来改善奶酪制作特性。

这些结果对实际选择的直接影响并不明显，这取决于对蛋白质组成的潜在溢价和牛奶加工业提出的激励措施，可是至少在奶酪生产中收集的牛奶中，最好更青睐酪蛋白而不是乳清蛋白，这种选择可以通过在基因组评估模型中包含影响单个蛋白质的变异来实现。