CN111243667B - 华西牛基因组选择方法 - Google Patents

Abstract

本发明公开了华西牛基因组选择方法，要解决的是现有肉牛育种群制种供种能力不足的问题。本发明的具体步骤如下：步骤一，构建参考群体并且测定经济性状；步骤二，对参考群体进行基因分型并且对数据进行处理和质量控制；步骤三，进行基因型数据填充，得到测序数据，然后进行QTL位点的筛选；步骤四，计算参考群体中各个性状的全部SNP标记的效应值；步骤五，得到候选群体基因组育种值，依据每个性状的基因组估计育种值计算个体的综合选择指数。本发明的目的是建立基于基因组选择技术的高效优质华西牛选育方法，全面提高我国肉牛种群选择效率，大幅节约育种成本，为优质华西牛的培育提供分子育种方法，推动肉牛种业的快速发展。

Description

华西牛基因组选择方法

技术领域

本发明涉及肉牛育种领域，具体是华西牛基因组选择方法。

背景技术

品种及种源问题是限制我国肉牛产业发展的首要因素。我国肉牛业现有肉牛育种群制种供种能力不足，种公牛主要依赖国外进口。整体肉牛存栏中，华西牛杂交牛为2788万头，占肉牛存栏的65％。国内现有涉及肉用公牛饲养的种公牛站38家，存栏采精公牛1900头，其中华西牛1100头，育种群每年仅能提供100头左右的合格公牛，更新的公牛年需进口250头左右，进口率高达70％。因此人们也在进行相关方面的研究。

发明内容

本发明实施例的目的在于提供华西牛基因组选择方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明实施例提供如下技术方案：

华西牛基因组选择方法，具体步骤如下：

步骤一，构建参考群体并且测定87个经济性状；

步骤二，对参考群体的每头牛采血保存，提取DNA并且利用生物芯片进行基因分型，对基因分型后的数据进行处理和质量控制；

步骤三，对处理后的基因分型数据进行基因型数据填充，得到测序数据，然后进行QTL(quantitative trait locus，数量性状基因座)位点的筛选；

步骤四，将筛选出的QTL位点与测序数据整合，得到标准数据，计算各个性状的全部SNP标记的效应值；

步骤五，对候选群体进行采血和提取DNA，利用高密度770K芯片进行基因分型，质量控制后填充至标准数据中，利用参考群体估计的全部SNP标记的效应值和QTL位点效应值，将候选群体的基因型向量与位点效应值向量相乘，即可得到候选群体的基因组估计育种值(GEBV)，依据每个性状的基因组估计育种值计算个体的综合选择指数(GCBI)。

作为本发明实施例进一步的方案：步骤二中的质量控制为剔除不合格的个体和SNP，生物芯片为Illumina BovineHD高密度芯片。

作为本发明实施例进一步的方案：步骤一中经济性状包括生长发育、育肥、屠宰、胴体、肉质、繁殖六类指标，共87个性状。

作为本发明实施例进一步的方案：步骤三中QTL位点包括与断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度显著相关的SNPs。。

作为本发明实施例进一步的方案：步骤四中采用优化贝叶斯算法(BayesBT)计算各个性状的全部SNP标记的效应值。

作为本发明实施例进一步的方案：优化贝叶斯算法中假设有效应的SNPs概率为(1-π)，没有效应的SNPs概率为π，效应方差服从逆卡方分布，其模型如下：

其中，y是表型观察值向量，b是固定效应向量；W为筛选的性状QTL的基因型关联矩阵，V为筛选的性状QTL的位点效应值，为固定效应向量；g_i是770K高密度芯片第i个标记的效应值，方差为

n是总的标记数；e是随机残差向量，方差为

是残差方差。

作为本发明实施例进一步的方案：步骤五中个体的综合选择指数(GCBI)公式如下：

其中Gebv_CE为产犊难易度基因组估计育种值，Gebv_WWT为断奶重基因组估计育种值，Gebv_{DG_F}为育肥期日增重基因组估计育种值，Gebv_CW为胴体重基因组估计育种值，Gebv_DP为屠宰率基因组估计育种值。

与现有技术相比，本发明实施例的有益效果是：

本发明的目的是建立基于基因组选择技术的高效优质华西牛选育方法，全面提高我国肉牛种群选择效率，大幅节约育种成本，为优质华西牛的培育提供分子育种方法，推动肉牛种业的快速发展。

附图说明

图1为华西牛基因组选择方法的流程图。

图2为华西牛基因组选择方法中加快选择效率的流程图。

图3为华西牛基因组选择方法中770K芯片数据质量控制前后SNPs在各常染色体上对比图。

图4为华西牛基因组选择方法中770K芯片标记检出率图。

图5为华西牛基因组选择方法中770K芯片SNP标记最小等位基因频率分布图。

图6为华西牛基因组选择方法中770K芯片标记间隔的概率密度图。

图7为华西牛基因组选择方法群体的连锁不平衡图谱。

具体实施方式

下面结合具体实施方式对本专利的技术方案作进一步详细地说明。

实施例1。

华西牛基因组选择方法，具体步骤如下：

步骤一，构建参考群体并且测定经济性状：某团队从2008年开始在内蒙古锡林郭勒盟乌拉盖管理区构建华西牛基础母牛群2320头，经过逐年扩群，截止2019年，华西牛基础母牛群体数已超过4000头。该基础母牛群的后代用于构建华西牛的参考群体，测定的表型数据包含生长发育、育肥、屠宰、胴体、肉质、繁殖6类87个重要经济性状；

步骤二，对参考群体的每头牛采血保存，提取DNA，用Illumina BovineHD(770k)高密度芯片进行基因分型，并采用plink 1.0软件对分型后的数据进行处理和质量控制；

步骤三，利用44头测序数据，利用Beagle软件将参考群770K芯片数据填充至测序数据(16165263SNPs)，进行断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度候选QTL位点筛选；

步骤四，通过对参考群中的个体进行SNP基因分型，将筛选出的QTL(quantitativetrait locus，数量性状基因座)位点与质量控制后的770K高密度芯片位点整合，利用Beagle软件将770K芯片数据填充至770K+QTL位点整合的数据。结合各表型数据，采用基于断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度5个性状QTL位点的优化贝叶斯算法(BayesBT)估计各个性状的全部SNP标记的效应值。BayesBT假设有效应的SNPs概率为(1-π)，没有效应的SNPs概率为π，效应方差服从逆卡方分布。其估计模型如下：

其中，y是表型观察值向量，b是固定效应向量；W为筛选的性状QTL的基因型关联矩阵，V为筛选的性状QTL的位点效应值，为固定效应向量；g_i是770K高密度芯片第i个标记的效应值，方差为

n是总的标记数；e是随机残差向量，方差为

是残差方差；

步骤五，候选群体牛只采血，提取DNA，利用Illumina BovineHD(770k)高密度芯片进行基因分型，将原始芯片数据导入GenomeStudio软件(参考基因组UMD3.1)进行质量控制以获得高质量SNP基因分型数据，并利用Beagle软件将候选个体基因型填充至770K+QTL位点整合的数据。利用参考群体估计的SNP位点效应值和QTL位点效应值，将候选个体的基因型向量与位点效应值向量相乘

即可得到候选群体的基因组估计育种值(GEBV)，依据每个性状的GEBV对计算个体的综合选择指数(GCBI)。

根据国内肉用种牛育种数据的实际情况，选取断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度共5个主要性状进行基因组遗传评估，各性状基因组估计育种值经标准化后，按经济加权比例进行加权，得到中国肉牛基因组选择指数GCBI。各性状的基本统计量如表1所示，基因组育种值估计准确性的评估是通过采用肉牛数量性状基因组选择BayesB计算软件V1.0进行5倍交叉验证获得。

表1各性状基本统计量及基因组育种值估计的准确性

中国肉牛基因组选择综合指数(GCBI)的公式如下：

其中Gebv_CE为产犊难易度基因组估计育种值，Gebv_WWT为断奶重基因组估计育种值，Gebv_{DG_F}为育肥期日增重基因组估计育种值，Gebv_CW为胴体重基因组估计育种值，Gebv_DP为屠宰率基因组估计育种值。根据候选群体计算的每个个体的GCBI，依据该参数的大小对候选个体进行选留(选择GCBI>120的个体)。

华西牛种牛培育方案如图2所示，利用常规育种技术进行华西牛种公牛培育，平均世代间隔要3.5-4年，若利用100头以上的后裔信息，育种值估计的可靠性可以达到0.75，至少需要5年时间才能培育一头优秀的华西牛种公牛。而利用基因组选择技术培育华西牛种公牛，2年就可以培育一头优秀的华西牛种公牛，大大缩短了世代间隔。由于候选种牛出生后就可以采血获得基因型，因此该华西牛基因组选择方法可以进行早期选择，大大节约了饲养成本，提高选择效率和准确性。

实施例2

步骤一，参考群体的建立及表型性状的测定：某团队组建的华西牛资源群体，该团队从2008年开始，在内蒙古锡林郭勒盟乌拉盖管理区构建了我国首个华西牛牛资源群体，经过逐年扩群，到2017年基础母牛群体数已超过4000头。在每年的7月或8月份，在乌拉盖管理区测量犊牛和基础母牛的体尺体重数据，并收集相应的繁殖记录数据。为了便于资源群体后期的生长发育性状、屠宰性状、胴体性状和肉质性状数据的测量和收集，每年10月份将乌拉盖地区5～9月龄的华西牛运送到某育肥屠宰一体化公司，所有华西牛在进场时，都要测量进场重和进场的体尺数据。在集中育肥期间，每经过三个月测量一次体重和体尺数据。该资源群体所有华西牛都依据统一的饲养管理方法进行集中育肥，同时饲喂的饲料也都是一样的。当华西牛集中育肥6个月时，静脉采血20ml并保存，并用2ml血液提取DNA(Deoxyribo Nucleic Acid，脱氧核糖核酸)，所有个体均用Illumina BovineHD(770k)高密度SNP(single nucleotide polymorphism，单核苷酸多态性以单个碱基变异为标志确定的多态性作为遗传标志)芯片获取SNP数据，同时又选择44头具有代表性的华西牛进行全基因组测序。在牛屠宰前一天，所有华西牛空腹24h，并测量屠宰前的活重和体尺。当资源群体的华西牛集中育肥10～12个月后，也就是在第二年的8-10月份进行分批屠宰，获取屠宰和胴体数据。屠宰数据包括屠宰牛在屠宰过程中的所有数据，便如：牛头重、牛舌重、牛蹄重、牛皮重、牛肾油脂肪、牛的红脏数据和牛白脏数据、睾丸等屠宰数据。屠宰后所有牛的胴体排酸48h后进行胴体分割，胴体的分割数据包括牛肉块的重量，例如：腹肉、脖肉、前腱子、金钱腱、上脑、眼肉、外脊、里脊、后腱子、米龙、黄瓜条等肉块。同时在胴体分割过程中，采集肉样标本(12～13肋眼肌一块约1Kg)，用于肉质性状的表型数据测量，肉质性状包括剪切力、系水力、ph、大理石花纹评分、胴体眼肌面积、肌内脂肪含量、肉样的脂肪酸性状等。由于断奶重作为一个综合指标，不仅反映犊牛断奶前一阶段的生长情况，更对犊牛断奶后的生长趋势起着决定性作用。产犊难易度反映牛只后代女儿的繁殖能力，育肥期日增重是评价华西牛生长发育的重要指标，在某种意义上育肥期日增重直接反应牛只对饲料的利用情况。胴体重和屠宰率是直接反映了肉牛相对胴体质量及肉牛生产性能的重要指标。因此，选取该五个性状作为重要的育种目标性状进行选种选配对于提高肉牛的生产性能具有重要经济学意义。

步骤二，基因分型与芯片数据的处理：将1407头华西牛全部用Illumina BovineHD(777962SNPs)高密度芯片进行基因分型，去除性染色体和质粒上的SNP(42669SNPs)，1407头华西牛高密度芯片在各常染色体上的分布如表2所示，共735293个SNPs。

表2芯片质量控制前各条常染色体的SNPs个数及比例

在进行数据分析之前，需要对基因型数据进行质量控制，利用PLINK V1.0剔除不合格的个体和SNP。本研究质量控制标准和代码为：plink--cow--file filename--geno0.1--maf 0.05--hwe 0.000001--mind 0,1--recode12--out filename。经过质控后，还剩1233头华西牛和597031个SNPs。质控后利用44头华西牛重测序数据，对华西牛群体的芯片数据进行填充，使用Beagle软件先对质控后缺失的SNPs进行填充，Beagle软件运行命令为：java–Xmx1000m-jar beagle.jar unphased＝file.bgl out＝output niterations＝100，然后利用Beagle软件分染色体将其填充至重测序数据(16165263SNPs)，运行命令为：java-d64-Xss100m-Xmx100g-jar beagle.08Jun17.d8b.jar gt＝chr1_vcf.gz ref＝chr1.gt.vcf.gz ne＝74map＝chr1_qc.map out＝impute_chr1.gt nthreads＝64，填充后各常染色体的数量及所占比例如表3所示。利用填充后的重测序数据，结合断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度性状进行候选QTL位点筛选，将筛选出的QTL位点与质控后的高密度芯片数据组合成标准数据，用于华西牛基因组选择。

表3重测序数据各常染色体的SNPs个数及比例

图3是1407头华西牛高密度芯片质量控制前后SNPs位点在各染色体上的分布情况，从图3中可以看出不合格SNP数在染色体上所占比例相当。

图4显示的是华西牛高密度芯片各常染色体上的SNP是否被检出的分布图，表明常染色体上的大部分SNP被检出，仅有少部分未检出。图5显示的是华西牛高密度芯片所有SNP标记(774660个SNP)的最小等位基因频率(minor allele frequency，MAF)分布图，可以看出，最小等位基因频率(<0.05)的稀有位点占一定比例，去除这些位点后，余下标记的MAF是均匀分布。图6是Illumina BovineHD芯片的SNPs标记间的概率密度图谱，从图6中可以看出绝大部分的SNP标记间距离小于10kb，同时间距在4kb范围的SNP最多。图7是华西牛群体的连锁不平衡图谱，从图7中可以看出，应用高密度芯片的华西牛资源群体，其SNPs平均间距<10Kb使得该群体处于较强的连锁不平衡状态。

本发明可以弥补我国肉牛基因组选择的空白，解决肉牛生产过程中缺乏基因组选择方法的问题，为我国肉牛业提供高效优质肉牛选育的技术手段，加快我国肉牛育种的进程，推动我国肉牛业的快速发展。结果证明，通过基因组选择方法培育华西牛新品种比传统育种技术增速50％，依据我国最大的华西牛参考群体2500头，通过自主研发的新方法(BayesBT)估计的基因组育种值准确性提高10～30％，通过这种方法可以快速地估计我国华西牛的基因组育种值，对于我国肉牛育种具有重大的应用价值和推广前景。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (5)

1.华西牛基因组选择方法，其特征在于，具体步骤如下：

步骤一，构建参考群体并且测定经济性状；

步骤二，对参考群体的每头牛采血保存，提取DNA并且进行基因分型，对基因分型后的数据进行处理和质量控制；

步骤三，对处理后的基因分型数据进行基因型数据填充，得到测序数据，然后进行QTL位点的筛选，所述QTL位点包括与断奶重、育肥期日增重、胴体重、屠宰率和产犊难易度相关的显著SNPs；

步骤四，将筛选出的QTL位点与测序数据整合，得到标准数据，计算各个性状的全部SNP标记的效应值；

步骤五，对候选群体进行采血和提取DNA，进行基因分型，质量控制后填充至标准数据，利用参考群体估计的全部SNP标记的效应值和QTL位点效应值，将候选群体的基因型向量与位点效应值向量相乘，即可得到候选群体的基因组估计育种值，依据每个性状的基因组估计育种值计算个体的综合选择指数，所述个体的综合选择指数GCBI公式如下：

其中Gebv_CE为产犊难易度基因组估计育种值，Gebv_WWT为断奶重基因组估计育种值，Gebv_{DG_F}为育肥期日增重基因组估计育种值，Gebv_CW为胴体重基因组估计育种值，Gebv_DP为屠宰率基因组估计育种值。

2.根据权利要求1所述的华西牛基因组选择方法，其特征在于，所述步骤二中的质量控制为剔除不合格的个体和SNP，基因分型采用Illumina BovineHD(770K)高密度芯片进行基因分型。

3.根据权利要求1或2所述的华西牛基因组选择方法，其特征在于，所述步骤一中经济性状包括生长发育、育肥、屠宰、胴体、肉质、繁殖六类指标。

4.根据权利要求1所述的华西牛基因组选择方法，其特征在于，所述步骤四中采用优化贝叶斯算法计算各个性状的全部SNP标记的效应值。

5.根据权利要求4所述的华西牛基因组选择方法，其特征在于，所述优化贝叶斯算法中假设有效应的SNPs概率为(1-π)，没有效应的SNPs概率为π，效应方差服从逆卡方分布，其模型如下：

其中，y是表型观察值向量，b是固定效应向量；W为筛选的性状QTL的基因型关联矩阵，V为筛选的性状QTL的位点效应值，为固定效应向量；g_i是770K高密度芯片第i个标记的效应值，方差为

n是总的标记数；e是随机残差向量，方差为

是残差方差。

Images