基于GBS简化基因组测序数据重建麦洼牦牛保种群系谱

李在文,李 响,李小伟,李 飙,江明锋*

(1.西南民族大学畜牧兽医学院,成都 610041;2.四川省龙日种畜场,红原 624401)

麦洼牦牛是乳肉兼用的高原型地方品种[1],乳蛋白质含量高,具有稳定的遗传特性,对高寒牧区有着良好的适应性,是我国珍稀的地方遗传资源[2]。麦洼牦牛是全国21个牦牛地方优良品种之一,列入《四川家畜家禽品种志》和《牛品种志》[1]。改革开放后,川西北高原牧区推广草场联产责任承包制极大的促进了牧区经济发展,同时由于包产到户导致麦洼牦牛群体变小,牧户养殖牦牛近交系数上升,麦洼牦牛繁殖性能和生产性能快速下降,最终导致品种退化[3]。为了对麦洼牦牛进行提纯复壮,四川省龙日种畜场于2004年组建了纯黑(毛色纯黑)、粉嘴(毛色纯黑,鼻唇部毛色为白色)和弗洛(出生时毛色为黑,随着年龄增长逐渐变白)3个核心群进行麦洼牦牛保种选育[4]。

育种场拥有准确的系谱对保证育种工作的正常开展、实施保种选育方案具有十分重要的意义[5]。随着分子生物学技术的快速发展,已有多种方法对系谱不清的选育群体进行系谱重建。各种DNA分子标记技术应用于动植物选育群体的系谱建立,能够快速且准确地构建群体的系谱,有利于更精确的计算育种值、近交系数和亲缘系数等育种参数,使育种工作更加高效[6-8]。DNA测序技术的快速发展,使获取大量SNPs信息变得更加便宜和便捷[9]。基因组测序、基因芯片等高通量检测技术在家畜近交系数、亲缘关系及遗传结构等研究中得到了广泛应用[10-12],为家畜的系谱重建、育种和保种奠定了更坚实的基础。2011年Elshire等[13]开发了GBS技术,与其他高通量测序技术相比,该技术具有步骤简单、快速、可重复性高和成本低等特点。目前,在群体遗传学、遗传图谱构建和亲缘关系鉴定等领域应用广泛[14-16]。

受技术及经济条件限制,牦牛选育过程中往往面临系谱不全甚至缺乏的情形,严重影响了牦牛育种进程,导致牦牛缺乏高产优良品种,我国仅有大通牦牛和和阿仕坦牦牛两个培育品种[17]。在麦洼牦牛育种过程中也发现了系谱信息完整度较低,缺少关键亲属关系信息等问题,严重迟滞了牦牛育种计划的制定和执行。本课题组以龙日种畜场的3个麦洼牦牛保种群为研究对象,进行简化基因组测序[18],通过获得的SNPs分析3个保种群的亲缘关系及近交程度,并尝试重新构建完整的系谱为后续的保种和育种工作服务。

1 材料与方法

1.1 试验材料

本试验于2020年6月在四川省龙日种畜场进行,以麦洼牦牛保种群为研究对象,从全黑群(QH)、粉嘴群(FZ)和弗洛群(FL)3个保种群中随机选取406头各年龄段的麦洼牦牛(具体数量见表1)。颈外静脉采血,EDTA抗凝,-80 ℃保存备用。

表1 3个保种群麦洼牦牛数量统计表Table 1 Statistics of the number of Maiwa yaks in the 3 preserved populations

1.2 试验方法

1.2.1 DNA提取 采用天根生物科技(北京)的血液DNA提取试剂盒,对全部血样进行DNA提取,采用Qubit®2.0测定DNA浓度,同时用1%的琼脂糖凝胶电泳检测。

1.2.2 文库构建及测序 参照文库构建方法[19],用MseI限制性核酸内切酶对检测合格的DNA进行酶切处理,在酶切片段两侧加上带有barcode的接头,然后PCR扩增,再混合扩增产物,选取所需片段建立文库。构建的文库使用Qubit® 2.0进行初步定量,稀释至1 ng·μL-1后使用Agilent 2100 Bioanalyzer进行库检[20]。插入片段长度符合预期后,采用q-PCR对文库的有效浓度进行准确定量(文库有效浓度>2 nmol·L-1),保证文库质量[21]。库检合格,依据不同文库的有效浓度及目标下机数据量的需求混池,随后进行Illumina Hi-seq PE150测序[22]。

1.3 数据统计及分析

1.3.1 数据质控及SNP检测 将测序得到的原始图像数据文件(Illumina Hi-seq测序平台)转化为原始测序序列,质控后获得有效序列。从NCBI获得野牦牛Bosmutus(PRJNA221623)参考基因组,使用比对软件Broadband Wireless Access (BWA)(参数为mem-t4-k32-M)[23]将获得的数据与参考基因组进行对比。采用SAMTOOLS[24]软件检测所有SNPs,根据以下条件过滤SNP位点:测序深度>3×,缺失率<20%,最小等位基因频率(minor allele frequency,MAF)>0.05。获得高质量SNP位点,用于后续研究。

1.3.2 群体遗传结构分析 利用GCTA[25]软件和R语言包(ggbiplot)基于个体间SNP的差异,进行PCA分析;采用TreeBest1.9.2[26]计算遗传距离矩阵,构建系统进化树。

1.3.3 近交系数的计算 基于个体间SNP差异,采用GCTA软件计算全部个体的近交系数(FIS)。将个体的近交水平分为以下4种近交程度[27]:①血亲程度FIS≥1/8;②近亲程度1/8≥FIS≥1/32;③中亲程度1/32≥FIS≥1/128;④远亲程度FIS≤1/128。

1.3.4 亲缘系数计算及亲缘关系G矩阵 基于血缘同源(identity-by-descent,IBD)信息,利用GCTA软件计算所有个体间的亲缘系数后用R version 3.6.3[28]进行聚类分析,并用pheatmap 1.0.12软件包绘制出406×406的矩阵。将亲缘系数按其大小分为亲子、全同胞、半同胞等不同的关系。

1.3.5 亲缘关系分析及整理方法 由于麦洼牦牛多为两年一胎或三年两胎[29],当二者年龄差小于3时必为同代,当年龄差大于3则需要进一步亲属关系判定。根据两两之间的亲缘系数和年龄差,初步判定其亲缘关系,具体判定条件见表2。

表2 亲缘关系判定条件Table 2 Kinship determination condition

根据上述判定条件,用Excel 2016整理亲缘系数,并进行性别、年龄、样本编号及亲缘系数等数据的处理。根据表2的条件,在Excel平台中利用Visual Basic编辑器设置一个宏代码,利用该宏代码实现判定和整理归纳亲缘关系信息的自动化。

1.3.6 家系构建 根据《畜禽遗传资源保种场保护区和基因库管理办法》[30]要求,构建家系要求在3个世代内没有血缘关系,则不同家系间亲缘系数应小于0.062 5。通过PLINK v1.90软件计算个体间亲缘关系,并利用最长距离法进行聚类分析。根据群体进化树的结果,利用MEGA v7软件建立体现个体遗传距离的聚类图,以此进行保种群的家系构建。

育种工作中通常以优良种公牛进行较大规模的配种。尤其在推行人工授精的育种场,优良种公牛拥有较多的后代数量往往有利于获得更好的遗传进展。因此,本研究对公牦牛血统进行了系统的研究。从保种群中选取4～8岁的公牦牛,并在保种群进化树图中标出这些公牦牛所在位置。根据MAF和缺失率进行SNP过滤,统计个体间的遗传距离,并根据NJ法进行聚类分析,绘制出公牦牛进化树图。

1.3.7 家系验证分析 为检验构建的家系效果,将个体按照构建的家系分组并进行遗传多样性分析,基于SNP采用GCTA软件计算每个个体的FIS、He、Ho,并将家系内个体的He、Ho、FIS取平均值作为每个家系的He、Ho、FIS。为衡量家系间的遗传分化程度,用R语言包计算家系间的群体遗传分化系数(FST)和群体遗传距离(DR)。FST的取值范围在0～1,FST=0～0.05:群体间遗传分化可忽略不计;FST=0.05～0.15:群体间存在中等程度的遗传分化;FST=0.15～0.25:群体间存在较高程度的遗传分化;FST≥ 0.25:群体间遗传分化程度极高,不同群体间的等位基因已固定。

2 结 果

2.1 基因组DNA提取

提取全部样本血液基因组DNA,DNA浓度范围在40～120 ng·μL-1,经电泳检测(图1)后发现DNA条带整齐,无拖尾。表明提取的DNA符合试验要求。

M为DNA相对分子质量标准;S为标准品;1～15为样品原液M is Trans 2 k plus; S is standard; 1-15 are samples stock solution图1 部分牦牛样本基因组DNA电泳图Fig.1 Genomic DNA electrophoresis of partial yak samples

2.2 测序数据产出与质控

本研究获得有效碱基数279.21～2 456.34 Mb,共产生327.517 Gb原始测序序列,平均为806.69 Mb;过滤后共产生327.512 Gb有效序列,平均为806.68 Mb。碱基检测错误率为0.04%～0.06%,平均为0.04%;Q20质量的碱基占比93.83%～97.58%,平均为95.73%;GC碱基含量为36.56%～40.28%,平均为38.88%;获得有效读长数1.94～17.08 Mb,平均为5.60 Mb;有效读长比对到参考基因上的比例为92.6%～99.68%,平均为98.75%;测序深度为5.65×～31.64×,平均为14.33×,平均覆盖度为3.46%。表明本次测序所获数据可满足后续分析要求。

2.3 SNP检测

利用SAMTOOLS软件检测获得数据的SNP位点,共鉴定出1 628 805个SNPs位点,经质控得到126 122个高质量SNPs位点。获得的高质量SNPs与参考基因组序列进行对比分析,图2为质控SNP位点数目及基因组序列长度,可看出SNP的数目随基因组序列片段长度的增加而增加,SNP分布相对均匀。

竖轴表示基因组序列长度,曲线表示SNPs数目The vertical axis indicates the genome sequence length and the curve indicates the number of SNPs图2 质控SNPs位点数与基因组序列片段长度变化趋势图Fig.2 Variation trend of SNPs number and genome sequence length

2.4 群体遗传结构分析

2.4.1 PCA分析 经计算PC1、PC2的贡献率分别为4.75%、2.72%。由图3可知,406头麦洼牦牛聚成不同类群,FL群聚类紧密,与另外两群部分个体聚在一起;QH群和FZ群存在明显的分层,群内个体分布不均匀,QH群部分个体脱离群体与FZ群聚在一起。

图3 406头麦洼牦牛PCA分析Fig.3 PCA analysis of 406 Maiwa yaks

2.4.2 群体进化树分析 从聚类分析NJ树(图4)可以看出,保种群可分为两大支,其中一大支由QH群聚类形成,有个体与FZ群和FL群聚类在一起;另一大支由FL和FZ群组成,分化过程中FZ群单独形成一支,FL群也分离出一小支,但FL群部分个体与QH群和FZ群聚在一起。

浅灰色.QH群;黑色.FL群;深灰色.FZ群Light grey. QH;Black. FL;Dark grey. FZ图4 406头麦洼牦牛系统进化树Fig.4 Phylogenetic tree of 406 Maiwa yaks

2.5 G矩阵分析

本研究利用GCTA软件计算得到164 836个亲缘关系对,并绘制出保种群的亲缘关系G矩阵,见图5。图中红色对称轴表示每个样本与其本身亲缘系数,白色大多聚在对称轴附近。左上角有两个白色方块构成的矩形区域,表示群落间的亲缘关系较近,发现这两个区域个体均属于QH群和FL群,这与群体进化树分析结果一致;右下角仅有少量白色方块,表示个体间亲缘关系较远。右上角零星散落的少许白色方块可能是近交导致的部分个体亲缘系数过大。

矩阵的右侧和下侧显示样本编号,每个小方块表示两个个体的亲缘关系,红色方块表示该样本与本身的亲缘系数,白色表示个体间亲缘关系较近,蓝色表示个体间无亲缘关系The right and lower sides of the matrix show the samples number, each small square indicates the affinity of two individuals, red squares indicate the affinity coefficient between that sample and itself, white squares indicate close affinity between individuals, blue squares indicate no affinity between individuals图5 保种群G矩阵图Fig.5 Conserved group G matrix visualization results

通过GCTA软件计算所有个体的近交系数(FIS),根据“1.3.3”判定个体的近亲程度,以此验证亲缘关系分析的可靠性。经统计,保种群共有11个血亲程度个体,84个近亲程度个体,121个中亲程度个体,剩下190个远亲程度个体。

2.6 亲缘关系分析及结果整理

按照“1.3.5”中的方法判定保种群的亲缘关系,共鉴定出134个全同胞关系,912个半同胞关系,136个疑似亲子关系或全同胞关系,520个疑似半同胞关系或叔侄关系,205个疑似半同胞关系或叔侄关系或祖孙关系。由于麦洼牦牛多为两年一胎或三年两胎,对年龄在3岁之内(同代)的亲缘关系判定的准确性较高,可直接判定出全同胞和半同胞关系,相关结果可直接用于保种选育。而另外几种亲缘关系(两代及以上)的判定需要进一步的遗传分析。

2.7 家系构建

2.7.1 保种群家系构建 麦洼牦牛保种群的聚类如图6所示。根据《畜禽遗传资源保种场保护区和基因库管理办法》,将龙日种畜场麦洼牦牛划分为12个家系,从右至左分别命名为G1～G12。从图6可以看出,大部分家系中的个体数量及位置分布较为均匀,但G5、G6、G8和G11家系的遗传距离较远,聚在该家系的个体数目较少。

纵坐标表示个体间遗传距离,横坐标表示每个样本编号,序号代表家系编号The vertical coordinate indicates the genetic distance between individuals, the horizontal coordinate indicates each sample number, and the serial number represents the family line number图6 麦洼牦牛家系聚类图Fig.6 Cluster diagram of Maiwa yaks

2.7.2 公牦牛家系构建 公牦牛分布进化树如图7A所示,其黑色标记表示公牦牛。从中可以直观看出麦洼公牦牛在保种群中分布比较均匀,但与群体进化树分析对比,有28头公牦属于QH群,超过半数;有10头公牦牛属于FZ群;FL群仅有4头公牦牛。

A. 42头麦洼公牦牛分布情况;B. 42头麦洼公牦牛聚类图:图右侧为公牦牛个体编号;图左侧两条竖线分别表示亲缘系数为0.062 5、0.125 0,与聚类图交叉即分为一个家系A. Distribution of 42 Maiwa male yaks; B. Pedigree clustering of 42 Maiwa male yaks: On the right side of the picture is the individual number of male yak; The two vertical lines on the left side of the figure indicate the kinship coefficients of 0.062 5 and 0.125 0, respectively, crossing with the cluster diagram is divided into a family图7 42头麦洼公牦牛分析结果Fig.7 Results of analysis of 42 Maiwa male yaks

公牦牛聚类如图7B所示,当亲缘系数为0.062 5时,可将42头麦洼公牦牛划分成3个血统,其中QH118与QH119各成一种血统。当亲缘关系设为0.125 0时,可将42头麦洼公牦牛划分成12个血统。从该结果中可知,FL群的公牦牛有2种血统,FZ群的公牦牛有3种血统,而QH群的公牦牛有7种血统。在排除个体数量较少的血统后,发现QH群的血统组成最完善,可分为6种血统,且个体数量较多。

2.8 家系验证分析

由表3可知,保种群平均FIS为0.021 3,G5、G6、G8、G11家系FIS均高于保种群平均FIS,表明这些家系近交程度较高,其中G8家系的FIS最大,为0.060 7;G10家系的FIS最小,为0.011 1。保种群12个家系的Ho=0.288 9～0.305 5,He=0.307 4～0.307 6,12个家系的Ho均低于He。

表3 12个家系的多样性分析Table 3 Diversity analysis of 12 families

家系间遗传分化(FST)和遗传距离(DR)见表4,可知保种群的FST在0.007～0.092,G5与G8家系的遗传分化程度最小(FST=0.007),G4与G12家系的遗传分化程度最大(FST=0.092)。G5家系的平均FST最小(FST=0.027),与G4成中等遗传分化程度;G4家系的平均FST最大(FST=0.065),与G2、G3家系遗传分化程度较小,与其他各家系均是中等遗传分化程度。

表4 遗传分化系数(FST)(对角线下)和遗传距离(DR)(对角线上)Table 4 F-statistic (below diagonal) and genetic distance (above diagonal)

3 讨 论

研究表明,SNP位点数是亲缘关系估计准确性的主要限制因素[31],而GBS技术可以短时间内获得大量SNPs位点。本研究经GBS测序并严格质控后检测出126 122个SNPs,麦洼牦牛系谱SNP芯片的制备可以此为基础。PCA图中个体空间直线距离与个体遗传关系成正比[32]。PCA分析中FZ和QH群分布不均匀,有分化趋势,这与钟光辉等[33]关于血液蛋白质多态性研究结果不符。其原因可能是龙日种畜场在育种过程中引入优良个体,以缓解麦洼牦牛在封闭选育中出现的近交系数上升。群体系统进化树分析中,大部分QH群个体聚类在一起,FL群小部分个体聚类紧密,表明QH和FL群的保种效果较好。其中FL群与另外两群的部分个体亲缘关系较近,形成较混乱的分支,这可能是由于龙日种畜场会从后备群挑选优秀种公牛补充进核心选育群,出现部分血缘交叉和个体离群现象。该结果与纪会等[34]的研究结果类似,表明GBS测序方法同样适用于麦洼牦牛遗传关系分析,同时表明当前麦洼牦牛不同保种群之间出现了明显的分离,还存在一定的基因交流。

本研究基于获得的SNPs用GCTA软件计算所有个体的近交系数和亲缘系数并构建G矩阵,获得了所有个体的亲属关系。根据近交系数发现,保种群有近1/4个体的近交系数较高,这些个体会严重干扰亲缘关系分析及判定,同时也不利于未来保种选育工作。通过对亲缘关系矩阵的分析,QH群和FL群都存在一个亲缘关系较近的群体,这与群体进化树分析结果一致。张金鑫等[35]和Islam等[36]也构建了G矩阵,说明利用亲缘关系对基因组选择分子育种工作有显著效果。根据获取的164 836个亲缘系数对,运用宏代码依据亲缘系数的大小进行亲属关系的判定,迅速获得了全同胞与半同胞关系和二代及三代内的亲属关系,但年龄段大于3岁以上的亲属关系判定准确率不高,该方法与赵书民等[37]在鉴定人类亲属关系时所使用的方法相似。为提高生产效益,成年牛会与不同异性个体交配形成多个半同胞个体,当半同胞个体数量大于等于4时,同胞关系数量会大于同胞个体数,因此本研究中鉴定出的912个半同胞关系和520个疑似半同胞关系或叔侄关系数量大于保种群数量406。本研究根据亲缘系数和年龄已经能够进行家系分群,并初步判定出群体中的疑似亲子关系或全同胞关系。在后续系谱构建过程中,只需对这些个体进行亲子鉴定即可构建准确的系谱。

本研究根据《畜禽遗传资源保种场保护区和基因库管理办法》的规定将麦洼牦牛划分为12个家系,与郝鑫宇等[38]研究结果相似,表明利用亲缘系数可以划分出不同的家系。进一步解析公牦牛血统分布,当亲缘系数划分为0.062 5时,可划分为3个血统,其中两个血统仅有1头牦牛,这2头牦牛可能是来自其他地区的优良个体。经询问证实龙日种畜场于2018年从群外引进3头优秀牦牛个体,这2头公牦牛可能来自这3头公牛。当亲缘系数划分为0.125 0时,可划分为12个血统,其中QH群可分为7个家系。当亲缘系数增加时血统数量也随之增加,这与王起山等[39]的研究结果相似,表明可根据遗传育种工作的实际需要来划分公牦牛血统。

为验证12个家系的构建效果,本研究对12个家系进行遗传多样性分析。当Ho低于He,则群体存在近交;当Ho高于He,则群体可能曾经导入外血[40]。本研究中12个家系的Ho均比He低,表明12个家系内均存在一定的近交,其中G5、G6、G8的Ho较小,表明近交程度较高。G11家系仅由2头1岁的QH群牦牛组成,G5、G6、G8家系中混杂着不同的表型,家系内存在一定的近交和基因交流,与刘彬等[12]的研究结果相似,表明近交和基因交流是采用群体继代选育法的必然结果。G1和G4家系个体均来自QH群,遗传分化程度较高,可作为QH群保种工作中的重点进行培养,这与马丽娜和马青[41]、赵志达等[42]和王继英等[43]的研究结果类似。

本项目探究了龙日种畜场麦洼牦牛的亲属关系,从基因组层面初步判定了同代、二代和三代的亲属关系,完善了麦洼牦牛系谱。但部分个体间的确切关系未知,本课题组将进一步开发更为简洁的系谱重建方法对牦牛系谱进行重建,以帮助目前开展的牦牛育种工作。

4 结 论

为厘清龙日种畜场选育群的系谱,本研究利用GBS技术获得的SNPs构建了亲缘关系G矩阵,通过Excel软件自动判定了保种群的亲属关系,并划分了12个家系,为后续系谱重建奠定了基础。为准确判定复杂的二代及三代亲属关系,应分析多态信息含量(PIC)、非父排除率(CEP)等数据,进行亲子鉴定,逐步完善个体系谱资料,为遗传资源的开发利用和新品种(品系)的人工选育提供理论基础。