基于mtDNA Cyt b序列变异探究柴达木黄牛的母系遗传多样性及遗传背景

杨秋蕾，魏旭东，马志杰,*，陈生梅，晁生玉，乌兰巴特尔

(1.青海大学畜牧兽医科学院，青海 西宁 810016； 2.青海省海西州农业技术推广服务中心，青海 德令哈 817099)

线粒体基因组(mtDNA)的结构特征与遗传特性使其在哺乳动物种群识别、遗传多样性保护以及母系起源进化等研究中能够作为可靠的分子标记而被广泛应用。mtDNACytb基因的进化速率适中，包含从种内到种间乃至科间的遗传信息[1]，因此，该基因已被广泛用于探究哺乳动物、鸟类和两栖动物等的母系遗传多样性、群体遗传结构及起源等研究[2]。

柴达木黄牛是《中国畜禽遗传资源志·牛志》中唯一收录的青海省黄牛地方品种[3]。先前对柴达木黄牛生化水平的研究中，张才骏等[4-5]对柴达木黄牛血液生化遗传标记进行了多态性检测，表明在柴达木黄牛中有11种生化标记存在多态性，发现柴达木黄牛与青海省东部黄牛存在最近的亲缘关系，与北方牛系蒙古牛、延边牛和安西牛具有较近的亲缘关系。才仁山德布等[6]对德令哈市56头柴达木黄牛和17头杂种黄牛的α-乳白蛋白(α-LA)和β-乳球蛋白(β-LG)进行分析，表明柴达木黄牛的α-LA和β-LG位点均具有多态性，但α-LA位点的基因杂合度很低。在柴达木黄牛父系分子遗传研究中，马志杰等[7-9]基于Y染色体USP9Y、ZFY-10和UTY-19标记的遗传变异对柴达木黄牛的父系遗传结构、多样性及遗传背景进行了探究，发现柴达木黄牛含有Y1和Y2两个普通牛父系支系，表明其为普通牛起源。在柴达木黄牛全基因组研究中，闫碧瑶[10]对9头柴达木黄牛进行全基因组重测序，并与秦川牛、西门塔尔牛、荷斯坦牛、日本和牛、娟珊牛、红安格斯牛和韩牛进行基因组比对分析，表明柴达木黄牛与以上7个黄牛品种具有明显的遗传差异，与秦川牛分化指数最小。Chen等[11]基于全基因组重测序方法对包括5头柴达木黄牛在内的世界家牛遗传多样性及起源进行综合分析表明，柴达木黄牛与西藏牛、蒙古牛、延边牛以及哈萨克牛间的亲缘关系依次变远。在中国其他黄牛母系遗传研究中，雷初朝等[12]分析了我国8个黄牛品种的mtDNA D-loop区全序列，确定了19种单倍型，提示中国黄牛mtDNA遗传多态性非常丰富；聚类分析揭示中国黄牛可能有3个母系起源，以普通牛和瘤牛起源为主。Lai等[13]指出中国黄牛中普通牛有T2、T3和T4三种单倍型组，瘤牛包含I1和I2两种单倍型组，瘤牛和普通牛的分歧时间大约是在公元前14 100—44 500年，提示中国南方黄牛和西南地区黄牛受瘤牛的影响大。Cai 等[14]分析了我国18个本地黄牛品种的Cytb基因序列，发现南方黄牛和北方黄牛分别以瘤牛(63.3%～100%)和普通牛(81.8%～100%)mtDNA占主导，中国黄牛分为普通牛和瘤牛2大支系。Xia等[15]分析了1 105个中国黄牛的mtDNA D-loop序列，结果表明，中国黄牛起源于普通牛和瘤牛，蒙古牛对中国牛的影响较大，并发现西藏牛中有牦牛基因渗入。

综上所述，研究者已基于生化标记、父系分子遗传标记以及全基因组序列变异初步探究了柴达木黄牛的遗传多样性、群体遗传结构和起源，但尚未见基于mtDNA序列变异对柴达木黄牛母系遗传多样性、群体结构及遗传背景等问题的研究，故本研究基于mtDNACytb基因的遗传变异探究柴达木黄牛的母系遗传多样性、群体结构和遗传背景等问题，为今后开展柴达木黄牛的保种及其合理开发利用提供基础数据。

1 材料与方法

1.1 样品采集及基因组DNA提取

在青海省海西蒙古族藏族自治州柴达木黄牛的5个主产区随机采集血样或耳组织，其中格尔木市88头(40♂、48♀)、乌兰县59头(15♂、44♀)、都兰县62头(23♂、39♀)、茫崖市33头(19♂、14♀)和大柴旦行政区26头(5♂、21♀)，共采集268头柴达木黄牛的样品，带回实验室。采样前询问牧户有关牛群的选育背景信息并查阅系谱记录情况，确保样品间无亲缘关系且具有代表性。使用基因组DNA提取试剂盒(北京艾德莱生物科技有限公司)提取基因组DNA，于-20 ℃保存备用。

1.2 mtDNA Cyt b基因的PCR扩增及产物纯化

参照Irwin等[16]报道的普通牛Cytb基因扩增引物，由生工生物工程(上海)股份有限公司合成柴达木黄牛Cytb基因的扩增引物，其中上游引物为：5′-CGAAGCTTGATATGAAAAACCATCGTTG-3′，下游引物为：5′-GGAATTCATCTCTCCCGGTTTACAAGAC-3′， PCR扩增产物预期长度为1 210 bp左右。PCR体系为：2×PrimeSTAR Max Premix (上海宝生物公司) 13 μL，上、下游引物 (10 pmol·L-1) 各1 μL，基因组DNA(25～50 pmol·L-1) 1 μL，超纯水9 μL，共25 μL。PCR反应程序为：94 ℃变性15 s，58 ℃退火30 s，72 ℃延伸45 s，35个循环，后冷却至4 ℃保存。PCR产物经1.5%琼脂糖凝胶(含GoldView核酸染料)电泳后，用凝胶成像系统进行检测并纯化回收，送北京擎科生物科技有限公司进行正、反向测序。

1.3 数据分析

将测序结果用Chromas 2.6.4软件进行核实和校正，得到每头柴达木黄牛的Cytb基因序列结果。用BioEdit7.2.5软件[17]进行多序列比对分析，检测核苷酸变异位点。通过DnaSP 5.0软件[18]依据该基因核苷酸变异确定不同个体的单倍型，再用 Arlequin 3.0软件[19]进行群体单倍型多样度和核苷酸多样度的计算分析。柴达木黄牛品种内群体间的固定指数(fixation index， Fst)通过Arlequin 3.11软件计算获得，以揭示品种内群体间的分化状况。Rst值通过线性化的Fst值计算获得，即Rst=Fst/(1-Fst)。基于Rst值使用MEGA 5.0软件[20]构建柴达木黄牛品种内不同群体间的亲缘关系聚类图，以揭示品种内各群体间的聚类情况和遗传关系。同时，从GenBank中下载普通牛(AY903438)、瘤牛(DQ186290)和牦牛(EU807949)的Cytb基因全序列作为参考序列，利用MEGA 5.0软件构建系统发育树。同时，不同单倍型序列间的中介网络图由Network 10.1软件[21]完成，以揭示不同单倍型间的系统发育关系。

2 结果与分析

2.1 柴达木黄牛母系遗传多样性分析

测序结果表明，柴达木黄牛Cytb基因编码区序列长度为1 140 bp。对268头柴达木黄牛mtDNACytb基因序列进行比对分析，发现茫崖、乌兰、都兰各群体均有1头黄牛为牦牛mtDNACytb序列类型，占总头数的1.12%，说明这3个县(市)的柴达木黄牛一定程度上受到了牦牛基因渐渗的影响。除去3条牦牛mtDNACytb类型序列，对共计265头柴达木黄牛Cytb基因序列进行多序列比对综合分析，共检测到29个核苷酸多态位点(图1)，其中单一多态位点4个，分别为121、242、573和906位点，占分析序列总长度的0.35%；简约信息位点25个，依次是66、67、93、232、312、345、384、592、621、633、774、795、813、834、835、997、1 035、1 050、1 066、1 080、1 092、1 104、1 112、1 114和1 116位点，约占分析序列长度的2.19%。共确定了12种单倍型，其中H2为优势单倍型，有161个个体共享。H1与H3各有31个和45个个体共享，H5和H4分别有10个和7个个体共享，H11单倍型仅有3个个体共享，H7和H10分别有2个个体共享，而H6、H8、H9和H12分别由1个个体拥有，为稀有单倍型。

“.”代表与H1单倍型序列中的核苷酸相同。“.”represented the same nucleotides as in H1 haplotype sequences.图1 柴达木黄牛Cyt b基因12种单倍型序列核苷酸变异及其分布Fig.1 Distribution and nucleotide variations of 12 haplotypes based on Cyt b of Qaidam cattle

由表1可以看出，柴达木黄牛的单倍型多样度为0.588 2±0.030 0，核苷酸多样度为0.004 0±0.002 2。在柴达木黄牛品种内，分别在格尔木、乌兰、都兰、大柴旦和茫崖群体中确定了4、3、9、4、7种单倍型。5个群体即格尔木、乌兰、都兰、大柴旦和茫崖各群体的单倍型多样度依次为0.599 5±0.040 4、0.329 7±0.073 3、0.713 1±0.042 3、0.452 3±0.108 8和0.695 6±0.079 6，而核苷酸多样度分别为0.103 99±0.062 4、0.011 20±0.013 2、0.188 61±0.102 0、0.103 99±0.062 4和0.179 22±0.099 1，表明都兰和茫崖这2个群体相比格尔木、大柴旦和乌兰这3个群体具有更高的母系遗传多样性，群体遗传变异相对丰富。

2.2 柴达木黄牛品种内群体间遗传分化及聚类关系分析

柴达木黄牛品种内5个群体间固定分化指数Fst值结果见表2，其中格尔木群体与乌兰群体间分化程度最大(Fst=0.161 8)，大柴旦群体和茫崖群体间的分化程度最小(Fst=-0.010 4)。

表1 基于Cyt b序列变异的柴达木黄牛核苷酸多样度和单倍型多样度大小

表2 柴达木黄牛品种内群体间Fst值(下三角)和Rst值(上三角)

格尔木群体与都兰、大柴旦、茫崖群体间的分化指数依次为-0.004 5、0.035 0和0.019 8，乌兰群体与都兰、大柴旦、茫崖群体间的分化指数分别为0.124 9、0.063 5和0.128 6，都兰群体与大柴旦、茫崖群体间的分化指数分别为0.002 9和0.000 2。

基于相应的Rst值用MEGA 5.0软件使用非加权组平均法(unweighted pair-group method with arithmetic means， UPGMA)构建柴达木黄牛品种内5个群体间遗传关系图，聚类结果表明：5个群体可聚为2类，其中格尔木群体与都兰群体最先聚在一起，大柴旦群体与茫崖群体也最先聚为一起，随后它们两者再聚为1类，而乌兰群体单独为1类，2类最后聚为1大类(图2)。

2.3 柴达木黄牛母系系统发育分析

以牦牛Cytb序列(GenBank登录号：EU807949)为外群，普通牛(GenBank登录号：AY903438)和瘤牛(GenBank登录号：DQ186290)相应序列作为参考序列，使用NJ(neighbor-joining)法构建系统发育树。同时，利用MJ(median-joining)中介网络图对柴达木黄牛不同单倍型间的系统发育关系进行分析，结果显示，柴达木黄牛12种单倍型包括9种普通牛单倍型(即H2、H3、H4、H5、H6、H7、H8、H10和H12)和3种瘤牛单倍型(即H1、H9和H11)(图3)，表明柴达木黄牛由普通牛和瘤牛2个母系遗传支系组成，具有2个母系起源。

图2 基于UPGMA法构建的柴达木黄牛品种内群体间遗传关系聚类图Fig.2 Clustering map of genetic relationship among populations of Qaidam cattle by UPGMA method

3 讨论

3.1 柴达木黄牛的母系遗传多样性

柴达木黄牛分布在青海省海西蒙古族藏族自治州柴达木盆地，是我国优良的地方黄牛品种之一，被列入省、国家级地方畜禽遗传资源保护名录[3]。近年来，柴达木黄牛这一地方黄牛品种存栏量从1980年的1.82万头下降到2011年的9 248头[22]。柴达木黄牛作为当地农牧民重要的生产生活资料来源，在地方农牧业发展和民族经济建设中曾发挥了重要的作用。但是，随着经济社会的发展特别是畜牧业发展需求的变化，在肉牛改良过程中，忽视了对柴达木黄牛品种资源的保护，导致柴达木黄牛杂种化程度变高、品种退化严重。因此，从分子水平上探究其遗传多样性和群体遗传结构是实施合理保护和利用的基础。在先前中国部分黄牛的研究中，Cai等[14]分析了中国18个本地黄牛品种136个个体的Cytb基因序列，表明18个中国黄牛品种的核苷酸多样度在0.002 9～0.009 2，单倍型多样度在0.333～1.000。在本研究中，排除3头牦牛mtDNACytb类型序列的柴达木黄牛个体后，对265头柴达木黄牛Cytb序列综合分析共确定了12种单倍型，其核苷酸多样度为0.004 0±0.002 2，单倍型多样度为0.588 2±0.030 0。和Cai等[14]的研究结果相比，除延边黄牛(单倍型多样度为0.333)和蒙古牛(单倍型多样度为0.524)的单倍型多样度值低于柴达木黄牛的单倍型多样度值(0.588 2)外，其他黄牛品种的单倍型多样度值均高于柴达木黄牛，说明柴达木黄牛相比中国大多数黄牛品种而言，其遗传多样性相对较低，但仍表现出较丰富的母系遗传多样性。在柴达木黄牛品种内，都兰群体的单倍型多样度值最高(0.713 1±0.042 3)，乌兰群体的单倍型多样度最低(0.329 7±0.073 3)，说明都兰群体的遗传变异最为丰富，乌兰群体的遗传变异最小。总体来看，都兰和茫崖这2个群体相比格尔木、大柴旦和乌兰这3个群体具有更高的母系遗传多样性，表明都兰和茫崖这2个群体遗传变异相对更为丰富。

图中圆面积大小与各单倍型频率大小成正比。The size of the circle was proportional to the frequency of each haplotype.图3 基于NJ法(a)和MJ法(b)构建的柴达木黄牛12种单倍型间系统发育关系Fig.3 Phylogenetic relationship among 12 haplotypes of Qaidam cattle based on NJ (a) and MJ (b) methods

3.2 柴达木黄牛品种内群体间遗传分化状况

固定指数(Fst)是衡量动物群体间分化程度的主要指标之一。在动物群体遗传学研究中，Fst值的大小反映了群体间分化程度的高低。在本研究中，柴达木黄牛品种内5个群体间Fst值大小各不相同，在-0.010 4～0.161 8，提示柴达木黄牛品种内群体间分化程度各不相同，存在差异。 Wright[23]认为，Fst值介于0～<0.05说明群体间分化程度很弱，介于0.05～<0.15说明群体间呈中等分化，介于0.15～0.25表明群体间分化程度大，而大于0.25说明群体间分化程度极大。依据此判定标准，对柴达木黄牛品种内5个群体间分化情况进行分析，说明格尔木群体与乌兰群体间分化程度大(Fst=0.161 8)，茫崖群体和大柴旦群体间的分化程度很弱(Fst=-0.010 35)。此外，大柴旦与乌兰群体间、都兰和乌兰群体间、茫崖与乌兰群体间的分化程度均达到中等水平，而格尔木与大柴旦、都兰及茫崖群体间、大柴旦与都兰、茫崖群体间以及都兰与茫崖群体间分化程度都很弱。鉴于柴达木黄牛这个珍稀的高原黄牛品种整体来看具有较高的品种内群体间遗传分化，故应及时设立保种场和保种区，建立保种核心群，加大柴达木黄牛的保种与选育实践。同时通过全基因组、转录组和多组学技术的联合分析挖掘其选育潜力，以发挥其种业优势。

3.3 柴达木黄牛的母系支系组成及遗传背景

在本研究中，使用系统发育树和中介网络图对柴达木黄牛的母系支系组成及遗传背景进行探究分析，结果显示柴达木黄牛由12种单倍型组成，且12种单倍型分成了2个母系遗传分支，分别由9种普通牛单倍型和3种瘤牛单倍型组成，表明柴达木黄牛由2个母系支系组成，拥有普通牛和瘤牛母系遗传血统，具有2个母系起源。在先前的研究中，Chen等[11]和Xia等[15]在中国黄牛的遗传多样性及系统发育分析中，发现西藏黄牛品种拥有一定程度的牦牛基因渗入。在本研究中，发现茫崖市、乌兰县和都兰县各有1头柴达木黄牛个体拥有牦牛mtDNACytb类型序列，占总头数的1.12%，说明柴达木黄牛品种中一定程度上还存在少量牦牛基因渗入，这说明在青藏高原高寒牧区，特别是在牦牛和黄牛共同生长地，牦牛和黄牛间的种间基因交流在柴达木黄牛、西藏牛这些高原黄牛品种中均普遍存在。

4 结论

柴达木黄牛品种具有较丰富的母系遗传多样性，品种内都兰和茫崖这2个群体相比格尔木、大柴旦和乌兰这3个群体具有更高的母系遗传多样性；柴达木黄牛品种内群体间分化程度存在差异，品种内部分群体间达到中等及中等以上分化程度；柴达木黄牛拥有普通牛和瘤牛母系遗传血统，由2个母系支系组成，具有2个母系起源，以普通牛起源为主；柴达木黄牛品种存在一定程度的牦牛基因渗入。