基于GBS簡化基因組測序評估3個麥洼牦牛保種群的遺傳結構研究

馬士龍，李小偉，李響，謝書瓊，劉益麗，唐嬌，江明鋒*

（1.西南民族大學畜牧獸醫學院，四川 成都 610041；2.青藏高原動物遺傳育種資源保護與利用國家教育部重點實驗室，四川 成都 610041；3.西南民族大學青藏高原研究院，四川 成都 610041；4.四川省龍日種畜場，四川 紅原 624401）

牦牛是非常寶貴的家畜遺傳資源基因庫。中國擁有極其豐富的牦牛資源，總量有1600多萬頭，占世界的90%以上，全國形成了18個地方品種和2個培育品種［1］，保護好牦牛品種資源對于牦牛產業可持續發展具有重要意義。麥洼牦牛，屬于乳肉兼用型地方優良品種，因中心產區是四川省阿壩藏族羌族自治州的紅原縣原麥洼部落而得名；據調查，麥洼牦牛是由阿壩牦牛引入康北牦牛、果洛牦牛的血緣，經長期選育形成［2］；1987年被《四川省家畜家禽品種志》［3］收錄，1988年載入《中國牛品種志》［3］。由于早期阿壩州草場超載過牧、飼養管理不足、牛群近親繁殖等導致麥洼牦牛群體有退化趨勢，國家于1953年建立了四川龍日種畜場，成立麥洼牦牛保種場［4］，并經長期的生產實踐，在2014年組建了純黑（全身被毛為黑色）、粉嘴（全身被毛為黑色，嘴唇毛色為粉色）和弗洛（全身被毛出生時黑色，逐漸變灰色）保種群［5］；其中，弗洛群的產肉性能高于粉嘴群和全黑群，純黑群的產奶性能高于粉嘴群和弗洛群。近年來，基于性能特征劃分的麥洼牦牛群體的生產性能均得到了提高［6］，但3個保種群的遺傳多樣性現狀，是否已經出現遺傳分化的問題一直未展開相關研究；因此，對麥洼牦牛進行遺傳多樣性和遺傳進化分析，評價不同保種群的保種效果顯得尤為重要。

自20世紀以來，國內學者已經通過微衛星標記［7-8］、血液蛋白基因座標記［9］、線粒體DNA標記［10-12］等方法對麥洼牦牛、九龍牦牛、斯布牦牛等品種進行群體遺傳進化分析，然而這些DNA標記僅能覆蓋部分且少量的牦牛基因組，所取得的基因組遺傳變異數據也是非常有限的。而新一代測序技術實現了單核苷酸多態性（single nucleotide polymorphism，SNP）分子標記的高通量檢測，可從整個基因組水平上準確高效地研究牦牛的遺傳進化關系［13］。基因分型測序（genotyping by sequencing，GBS）技術是在二代測序基礎上發展的簡化基因組測序技術，通過酶切加tag的方法進行測序，獲取酶切位點附近基因序列信息，進而檢測大量且準確性高的變異SNP位點信息［14］，被廣泛應用于遺傳圖譜構建、全基因組關聯分析和群體遺傳分析等領域［15-18］。

本研究利用GBS技術研究了麥洼牦牛3個保種群的遺傳多樣性、群體結構和親緣關系，評價保種效果，為今后麥洼牦牛遺傳資源的保護與利用提供了可借鑒方法和理論依據，為下一步進行麥洼牦牛的良種選育和生產性能提高奠定基礎。

1 材料與方法

1.1 試驗動物和方法

試驗于2020年6月在四川省阿壩藏族羌族自治州紅原縣國家級肉牛核心育種場—龍日種畜場進行，該地區處于查真梁子山北部區域，海拔3500～4000 m，地勢平坦；只有冷、暖季之分，屬大陸性高原寒溫帶季風氣候，牧草豐茂，是全國著名的地方優良品種-麥洼牦牛、藏綿羊、河曲馬主產區［4］。從3個麥洼牦牛保種群全黑群（QH）、粉嘴群（FZ）、弗洛群（FL）中隨機選取406頭麥洼牦牛（公∶母=1∶3.5）（表1），群體內個體間血緣關系未知，采集血樣樣本。具體步驟為用一次性注射器從牦牛頸外靜脈采集5 mL血液于乙二胺四乙酸（ethylene diamine tetraacetic acid，EDTA）抗凝管，放入冰盒中帶回，放置-80℃冰箱中保存備用。

表1 3個保種群的采集信息Table 1 The sample information of 3 breeding populations

1.1.1 主要試劑及儀器 血液基因組DNA提取試劑盒（D3392-01）購自Omega Bio-Tek公司；異丙醇、無水乙醇、瓊脂糖、50×TAE、6×Loading Buffer、DNA Marker（BM5000）等購自上海生工生物工程技術服務有限公司和愛必信（上海）生物科技有限公司。

主要儀器有離心機（Eppendorf，5804R，德國）、電子天平（Denver，TP-214，美國）、電泳儀（Bio-rad，美國）、凝膠成像系統（Invitrogen，美國）、紫外分光光度計（Thermo Scientific，NanoDrop 2000，美國）；Qubit?2.0熒光計（Invitrogen，美國）。

1.1.2 DNA提取 按照血液基因組DNA提取試劑盒的說明書提取406個血液樣本的基因組DNA，利用1%瓊脂糖凝膠電泳（100 V，40 min）和紫外分光光度計檢測DNA的純度（A260nm/A280nm值）和完整性，使用Qubit?2.0進行DNA濃度測定，樣品純度要求A260nm/A280nm值為1.8～2.0。

1.2 GBS試驗流程

1.2.1 GBS文庫構建 文庫構建參照董世武等［18］的方法，檢測合格的基因組DNA用Mse I限制性核酸內切酶進行酶切，酶切后基因組片段兩端加上帶有barcode的接頭（120 bp）后，利用PCR對每個樣品進行擴增，然后對樣品進行混合，選擇需要的片段（265～315 bp）進行文庫構建。

1.2.2 文庫庫檢及上機測序 構建好的文庫用Qubit?2.0進行初步定量，具體操作為稀釋文庫至1 ng·μL-1，再用Agilent 2100檢測文庫的插入片段大小是否符合265～315 bp，然后為保證文庫質量，使用Q-PCR方法對文庫的有效濃度進行準確定量（文庫有效濃度＞2 nmol·L-1）。庫檢合格后，把不同文庫按照有效濃度及每個樣本可產生10萬個標簽的目標下機數據量的需求進行樣本混合，即混池，然后進行Illumina Hiseq PE150雙末端測序。

1.3 數據統計與分析

1.3.1 測序數據及SNP的質控 Illumina Hiseq測序平臺得到的原始圖像文件經堿基識別分析轉化為測序序列數據（raw reads），這些raw reads需進行質量控制，剔除含有接頭序列的reads和低質量reads之后得到的數據即為有效數據（clean data）。

以野牦牛（Bos mutus）（生物項目數據庫：PRJNA221623）為參考基因組，通過BWA軟件（參數：mem-t 4-k 32-M）［19］將雙端有效數據比對到參考基因組，統計平均比對率、平均測序深度和平均覆蓋度。采用SAMTOOLS［20］軟件進行群體SNP的檢測，得到位點變異文件，以測序深度（sequencing depth）為3，缺失（miss）為0.2，次要等位基因頻率（minor allele frequency，MAF）大于0.05條件過濾SNP位點，最后獲得高質量的SNP位點用于后續分析研究。

1.3.2 群體結構與遺傳多樣性分析 基于個體間SNP差異，通過GCTA［21］軟件計算特征向量和特征值，用R語言繪制第一主成分與第二主成分的主成分分析（principal component analysis，PCA）二維圖。使用TreeBest［22］軟件中的P-distances模型計算遺傳距離矩陣并構建系統進化樹，引導值（bootstrap values）設置為1000，檢驗可靠性。采用ADMIXTURE軟件構建群體遺傳結構和群體世系信息，假定的祖先群體個數為2～20，根據交叉驗證錯誤值（cross validation error，CV）確定最優K值（一般CV最小的為最優K值）。使用PopLDdecay（v 3.40）［23］進行連鎖不平衡分析（linkage disequilibrium，LD），并用R軟件繪制LD衰減圖。

利用PLINK v 1.90［24］計算觀測雜合度（observed heterozygosity，Ho）和期望雜合度（expected heterozygosity，He）等遺傳多樣性指標。使用R包genepop統計分析每個SNP位點在所有個體中的近交系數（inbreeding coefficients，Fis）和群體遺傳分化系數（genetic differentiation index，Fst）。基于篩選后高質量SNP信息，進行選擇消除分析，采用PLINK軟件進行選擇信號檢測，以100 Kb為窗口、10 Kb為step進行滑動的核苷酸多樣度（π）比率和Fst計算，分別以π比率和Fst值top 10%作為閾值，對群體選擇信號進行分析，其top 10%的交集為候選位點，并對受選擇位點所在基因進行GO富集分析。

2 結果與分析

2.1 血液基因組DNA提取

提取406頭牦牛樣本血液基因組DNA，經1%瓊脂糖凝膠電泳檢測發現，條帶清晰明亮（圖1），各項參數均符合測序所需標準。

圖1 部分樣本血液基因組DNA電泳圖Fig.1 Genome DNA electrophoresis of some blood samples

2.2 測序數據產出與質控

對406份牦牛血液基因組進行GBS測序和質量控制（詳見1.3.1），獲得有效堿基數279.211～2456.350 Mb，共產生327.517 Gb原始數據，平均為806.690 Mb；過濾后共產生327.512 Gb有效數據，有效率達99.990%，平均為806.680 Mb。堿基檢測錯誤率為0.040%～0.060%，平均僅有0.044%（表2）；達到Q20質量的堿基占比為93.83%～97.58%，平均為95.73%；GC堿基含量為36.56%～40.28%，平均為38.88%；獲得有效讀長數為1.94～17.08 Mb，平均為5.60 Mb；有效讀長比對到參考基因上的占比為92.60%～99.68%，平均比對率為98.75%；測序堿基總量與基因組的比值，即測序深度最小為5.65×，測序深度最大為31.64×，平均測序深度為14.33×，平均覆蓋度為3.46%（表3）。表明測序質量高，GC分布正常，樣本均沒有被污染，此次建庫測序成功。

表2 麥洼牦牛測序數據統計Table 2 Summary of sequencing data of Maiwa yak

表3 測序數據質量評估Table 3 Quality assessment of sequencing

2.3 SNP檢測

應用SAMTOOLS等軟件對測序有效數據進行SNP檢測，共鑒定出1628805個SNP位點，進行質控后獲得高質量SNP位點126122個。質控后的SNP經過與牦牛參考基因組（生物項目數據庫：PRJNA 221623）比對分析，基因組的SNPs數目及基因組序列長度如圖2所示，總體上來看SNPs的數目與基因組序列長度總體上呈正相關；序列號為NW_005393843.1的片段最大（8785736 bp），所含SNP數目最多，為463個SNP位點；而由于基因組序列長度差異，不同基因組序列的SNP數目變化較大，NW_005393917.1的 片 段 較 大（542197 bp），卻 僅 有1個SNP位點。

圖2 質控SNPs位點數與基因組序列片段長度變化趨勢Fig.2 Variation trend of SNPs number and genome sequence length

2.4 遺傳多樣性和群體結構分析

2.4.1 遺傳多樣性分析 根據變異位點分析QH、FZ、FL群的雜合度（Ho、He）、核苷酸多樣度（π）、近交系數（Fis），結果見表4。3個群的觀測雜合度分別為0.3029、0.3042、0.3044，觀測雜合度均與期望雜合度接近，只有FZ的觀測雜合度稍高于期望雜合度，總體平均觀測雜合度與期望雜合度接近。3個群的核苷酸多樣度分別為0.3030、0.3020、0.3009，FL群最低；而與QH、FZ群的近交系數相比，FL群是最高的，為0.0209。

表4 麥洼牦牛3個育種群遺傳多樣性參數Table 4 Genetic diversity parameters of three preserved populations of Maiwa yak

連鎖不平衡（LD）模式見圖3，3個群體LD均隨距離增長逐漸平滑下降，各等位基因不存在強烈的連鎖不平衡水平（r2≤0.5），SNP距離 在0～100 kb的LD衰減 速度較快，LD系數r2約從0.5降至0.2以 下，SNP距離在100～300 kb的LD水平較低。當r2=0.3時，各群體對應LD衰減距離基本相等，說明各群體LD衰減速度差異不大。

圖3 LD隨距離增長的衰減Fig.3 Attenuation of LD with distance

由表5可以看出，3個群體遺傳分化系數Fst為0.02504～0.03513（介于0～0.05），其中QH和FZ群之間的遺傳分化程度最大（0.03513），兩群的遺傳距離同樣是最遠（0.0358）。

表5 群體間的遺傳分化系數（Fst）（對角線下）和遺傳距離（DR）（對角線上）Table 5 Genetic differentiation index（below diagonal）and genetic distance（above diagonal）among three populations of Maiwa yak

2.4.2 群體遺傳結構和親緣關系分析 基于個體基因組SNP差異程度，通過PCA二維圖的形式來展示個體的聚類情況，可以大致了解群體的遺傳結構。由圖4可知，麥洼牦牛群體的PC1、PC2的貢獻率分別為4.75%、2.72%，FL群緊密聚在一起，與QH、FZ群部分個體關系近；FZ群個體較分散，QH群中有個體偏離群體與FZ群聚在一起，兩群存在明顯的分層。

圖4 PCA分析二維圖Fig.4 PCA analysis of two-dimensional graph

群體遺傳結構指遺傳變異在物種或群體中的一種非隨機分布。按照地理分布或親緣關系可將一個群體分為若干亞群，處于同一亞群內的不同個體親緣關系較近，而亞群與亞群之間則親緣關系稍遠。從圖5可以看出，3個保種群的血統構成，當假設祖先群為2（K=2）時，QH、FZ、FL群都攜帶兩個祖代血緣；當最 佳K值 為12時 ，FL群 可 分 為 兩 個群 ，FZ群 部 分 個體（棕色）可單獨劃分為一個群，QH群部分個體（粉色）可單獨劃分一個群，其他需根據系譜記錄進行劃分。

圖5 麥洼牦牛structure群體遺傳結構Fig.5 Genetic structure of Maiwa yak structure populations

406頭麥洼牦牛聚類的系統進化樹（圖6）直觀地以圈圖的形式展示它們之間的親緣關系遠近，結果可以看出，QH群聚成一個大支（綠色），群體存在離群個體；FL和FZ群組成另一支，FZ群逐漸分離出來，最后形成一小支（藍色），FL群有一支最后被分離出來（紅色），但整體有一半的個體與QH、FZ群聚在一起。

圖6 系統進化樹Fig.6 Phylogenetic tree

2.4.3 群體選擇信號分析 Fst和π分析是兩種有效地檢測選擇消除區域的方法，特別是在挖掘與生存環境密切相關的基因功能區時可以得到較強的選擇信號，便于目標基因的篩選。根據麥洼牦牛群體間Fst和π進行選擇信號分析，篩選出3513個差異較大的SNP位點，共有104個差異基因受到選擇，進行GO功能注釋分析，結果表明（圖7），這些受選擇基因主要富集在有絲分裂、氯離子運輸、節律過程、刺激反應等生物學過程；主要涉及的細胞組成：氯離子通道復合體、細胞間連接、膜的錨定組分等；主要涉及的分子功能包括：氯離子通道活性、酶結合、蛋白質二聚化活性等。KEGG結果富集到了28條通路，主要集中在生殖激素信號、內/外分泌、鈣離子信號、cGMP-PKG信號、血管平滑肌收縮等調控通路上（表6）。

表6 受選擇基因的KEGG通路Table 6 KEGG pathway of differential expression genes

圖7 受選擇基因的GO富集條目Fig.7 GO enrichment items of selected gene

3 討論

牦牛是我國青藏高原珍稀而寶貴的遺傳資源，在高原畜牧業發展和生態系統平衡方面起著至關重要的作用。牦牛遺傳多樣性水平和遺傳結構的研究不僅是保護牦牛遺傳資源的主要內容，也是實施科學有效的品種保護措施的前提條件。隨著分子生物技術的快速發展，基于SNP的新一代測序技術成為主流方法，使得從基因組層面上研究牦牛遺傳結構和遺傳多樣性成為可能，更全面精確地監測牦牛保種效果。自麥洼牦牛品種被挖掘以來，其遺傳多樣性研究是保護這個地方牦牛品種資源的核心問題之一；因此，本研究對406頭麥洼牦牛進行GBS測序探究其群體的遺傳多樣性及遺傳結構，測序結果得到了806.68 Mb的有效數據量，共鑒定出126122個SNP位點，這些數據將為未來制備麥洼牦牛育種SNP芯片提供基礎。

遺傳雜合度是度量群體遺傳多樣性的最適參數，其觀測雜合度和期望雜合度的接近程度可反映群體遺傳多樣性的高低。本研究中麥洼牦牛3個保種群的觀測雜合度為0.3029～0.3044，屬于中度雜合群體，整個群的平均觀測雜合度與平均期望雜合度接近，這表明麥洼牦牛品種具有較高純度和豐富的遺傳多樣性，這與趙芳芳［25］的微衛星法和師方等［11］的隨機擴增多態性脫氧核糖核酸（random amplified polymorphic DNA，RAPD）法結果相似；龍日種畜場為高原淺丘地貌，可利用草地面積1.667萬hm2，在這樣的相對隔離地理環境中麥洼牦牛得到了有效保護，較低的近交水平也說明保種群并沒有受外來品種及近交的影響。另外，粉嘴群的觀測雜合度高于全黑群，結果與侯孟典等［26］和柴志欣等［27］的報道一致；這表明粉嘴群所經過的長期定向人工選擇強度高于全黑群，有研究［5］指出粉嘴個體相較于全黑個體的乳脂率較高，而擠奶量較低，推測粉嘴群已經在乳質等生產性狀上受到了強烈的選擇。群體間的遺傳分化與選育的方式和世代長短有關，四川省龍日種畜場的牦牛原種場采取的開放式群體繼代選育法，在世代選育過程中，按需地從麥洼牦牛主產區的選育群引進優良非近親種公牦牛進行血緣更新［3］，使得場內的麥洼牦牛群體內個體間遺傳變異度變高；因此粉嘴和全黑群的遺傳分化指數和遺傳距離較遠，分化呈現出不穩定的混亂狀態，有部分血緣交叉和個體離群，這提示在保種過程中，應該加強對種公牛的嚴格選留工作，避免近親交配的發生。另外，早期的弗洛個體其實并沒有從總群中單獨劃分出來，而是后來依據其毛色特征和產肉性能才組建的［5］，這印證了弗洛群是粉嘴群和全黑群的祖代個體雜交所形成的群落的可能；進一步說明弗洛群含有部分全黑群和粉嘴群血統，有助于解釋弗洛群的出生毛色變化現象。龍日種畜場內麥洼牦牛群體的遺傳結構隨著時間的推移也發生了較大變化，推測是由較小的保種群體規模限制導致，群體結構圖和系統進化樹都表明全黑群大部分個體聚為一類，弗洛群小部分個體聚為一類，這說明全黑、弗洛群的保種效果相對較好，全黑群和弗洛群出現血統較純正的個體群，這些個體和粉嘴群一樣，它們的某些經濟性狀可能已經得到了固定。

為鑒別3個保種群之間在基因組特定位點上存在的分化，對麥洼牦牛群體104個受到強烈選擇的基因進行GO和KEGG分析，發現這些基因廣泛參與生殖機能調控、生長發育、刺激反應等生物學相關過程。其中，特異型鈣調磷酸酶PPP3CC基因參與精子發生和睪丸成熟，缺失該基因的小鼠精子尾部中段會發生僵直、超活化能力降低［28］，該基因也已經在牛體內被檢測到［29］；鉀離子通道宿主因子KCNMA1與牛的繁殖性能相關，可增強動物性行為［30］。神經生長因子NGF能有效促進卵泡形成，通過調控原始卵泡中扁平細胞分化成次級卵泡中包裹卵母細胞的立方形顆粒細胞過程，在牛孤雌胚胎的早期發育起關鍵作用［31］；鳥嘌呤核苷酸結合蛋白GNAQ基因是通過影響小腦分泌激素從而影響動物毛色形成的重要調控基因［32］，也參與動物性腺的發育和生殖細胞的成熟［33］。這表明3個保種群在龍日種畜場的人工選擇條件下，公母牛的繁殖力均得到了增強，尤其公牛的性行為方面較為明顯。肌細胞增強因子MEF2C主要參與調控肌細胞分化，特別是在從胎牛到成年牛發育過程中，介導骨骼肌、心肌和平滑肌的細胞分化［34］；小GTP結合蛋白Rho下游效應因子ROCK2對低氧環境誘導的肺動脈高血壓發揮著重要作用，并在成肌細胞分化中持續上調表達，參與肌肉形成過程［35］；另外，KCNJ3基因與心律發生有關［36］，這3個基因可能與麥洼牦牛的高原適應性和肌肉發育調節相關，說明保種群在強烈的人工選擇中有向肉質性狀方面發生了定向選擇的趨勢。谷氨酸突觸通路中的ITPR2基因是鈣離子跨膜轉運活性的關鍵調控因子，內含子區存在動物抗環境壓力的潛在靶點［37］，這表明保種群在龍日種畜場的天然草場放牧和一定的補飼管理下，群體的抗應激能力也得到了強化。酪氨酸蛋白激酶受體KIT基因是調控哺乳動物白色皮膚及被毛的基因，KIT基因突變會導致動物由花斑到全白不同程度的白色被毛，而牦牛的白色斑點就是由KIT基因在6號染色體上的易位形成的［38］，提示該基因可作為弗洛群的灰色表型發生的候選基因；另外，仍有許多受選擇功能基因有待進一步研究注釋。本研究在基因組水平上進一步說明龍日種畜場內麥洼牦牛不同保種群經歷了定向選擇，篩選的受選擇功能基因對了解麥洼牦牛的適應性變化并獲得其優良的經濟性狀具有非常深遠的意義。

4 結論

本研究利用GBS技術揭示了麥洼牦牛3個保種群的遺傳多樣性和遺傳結構，粉嘴群和全黑群有遺傳分化趨勢，分群的效果得到了體現，3個保種群的保種效果較好，可在世代之間進行比較，進一步分析保種效果。測序出3513個顯著SNP位點和104個差異基因；3個保種群在毛色、肉質、繁殖性狀上經歷了人工選擇，這可為麥洼牦牛的保種選育工作提供建議。