基于體細胞評分的中國荷斯坦牛乳房炎抗性全基因組關聯分析

王曉，張勤，俞英

（中國農業大學動物科技學院，北京 100193）

基于體細胞評分的中國荷斯坦牛乳房炎抗性全基因組關聯分析

王曉，張勤，俞英

（中國農業大學動物科技學院，北京 100193）

中國荷斯坦牛；乳房炎抗性；SCC對數轉化；Case-control；全基因組關聯分析

0 引言

【研究意義】目前，全球約有1/3的奶牛患有不同類型的乳房炎，乳房炎給奶牛養殖業造成了巨大的經濟損失[1-2]。由于乳房炎的遺傳力較低[3]，直接選擇低遺傳力的臨床乳房炎并不是遺傳改良的有效方法，而是通過對其他高遺傳力且與乳房炎有高遺傳相關的性狀進行間接選育。標記輔助選擇（marker-assisted selection, MAS）可將乳房炎抗性基因和其它分子標記結合入現有的選育方案進行選育。因此通過全基因組關聯分析（genome wide association study, GWAS）尋找到顯著的單核苷酸突變（single nucleotide polymorphism, SNP），然后通過MAS進行奶牛乳房炎抗性的選育。【前人研究進展】AXFORD等指出，奶牛生產性能測定（dairy herd improvement, DHI）中的體細胞數（somatic cell count, SCC）是用于提高乳房炎抗性的首選性狀[4]，RUPP的研究發現，對低SCC的奶牛進行選育后，其乳房炎的發病率降低[5]。WIJGA通過全基因組關聯分析報道，與泌乳期平均體細胞評分（lactation-average SCS，LASCS）和測定日記錄體細胞評分標準差（test-day SCS standard deviation，SCS-SD）顯著相關的SNPs位于4，6和18號染色體上[6]。【本研究切入點】首先對北京地區中國荷斯坦牛的SCC進行對數轉化（Log- transformed），獲得SCS；再依據LASCS和SCS-SD高低將牛只劃分為乳房炎易感牛（case）及抗性牛（control）；進一步利用Case-control方法，將乳房炎易感性及抗性性狀與Illumina54k芯片的54 001個SNPs進行全基因組關聯分析。【擬解決的關鍵問題】通過全基因組關聯分析找到與荷斯坦牛乳房炎易感性及抗性存在顯著關聯的SNPs，揭示與乳房炎抗性密切相關的基因變異，為發現乳房炎易感性及抗性相關的分子標記及乳房炎分子抗病育種奠定理論依據。

1 材料與方法

試驗于2014年3—5月在中國農業大學動物科技學院分子數量遺傳學實驗室進行。

1.1 表型和基因型

2 093頭中國荷斯坦母牛DHI數據及其54K SNP基因型數據（BovineSNP50, Illumina, USA）。

1.2 數據處理

本研究利用SAS（Statistical Analysis System, Cary, North Carolina, USA）9.1.3軟件進行數據篩選。

原始數據為北京地區36個奶牛場2093頭中國荷斯坦母牛的頭3個胎次的DHI記錄，共計35 522條測定日記錄。其中SCC（單位：×103個/mL）的最小值為0，最大值為9 952，平均值為237.8，標準差為550.0。產犢日期和DHI測定時間段為2000年至2008年，測定天數為5—305d，測定間隔為0—215d，所有母牛來自于14頭公牛家系。

依據LASCS和SCS-SD將泌乳牛劃分為Casecontrol兩個子數據集[7]。劃分標準為：將LASCS和 SCS-SD進行半個標準差（Half of standard deviation, 0.5 SD）和一個標準差（One standard deviation, 1 SD）的劃分，高LASCS和高SCS-SD的牛歸到易感組（Case），低LASCS和低SCS-SD的牛則歸入抗性組（Control）（圖1），用于Case-control關聯分析方法的LASCS和SCS-SD的統計性描述如表1所示。

1.3 質量控制

將54 001個SNPs進行質控，剔除不符合條件的SNPs。剔除條件：（1）SNPs的call rate＜90%；（2）嚴重偏離哈迪-溫伯格平衡（HWE）（P＜ 10E-6）；（3）最小等位基因頻率（MAF）＜0.03。經過質控之后，共有43781/43671（43817/43704）個SNPs分別可用于LASCS（SCS-SD）半個標準差/一個標準差的關聯分析，其所在牛染色體上的分布情況如表2所示。

1.4 關聯分析

本研究的關聯分析方法是基于前期使用的ROADTRIPS軟件（版本 1.2）[7-8]，該軟件可以通過大量基因型信息并結合系譜信息估計出一個可以反應個體間遺傳聯系的經驗協方差矩陣。

ROADTRIPS（Robust Association-Detection Test for Related Individuals with Population Substructure）一共包括3種檢驗：RM（ROADTRIPS-MQLStest）檢驗、RCHI（ROADTRIPS-χ2test）檢驗和RW（ROADTRIPSWQLStest）檢驗。該軟件通過構建服從自由度為1的卡方分布的統計量，計算出每個檢驗的SNP的P值。

1.5 關聯分析結果的顯著性檢驗

為了降低多重檢驗帶來的假陽性升高問題，采用Bonferroni方法對關聯分析結果進行校正。如果設定0.05為顯著性水平，則犯I類錯誤的累積概率就要控制在0.05以內。所以當進行m次顯著性水準為α的假設檢驗時，即有Bonferroni不等式0.05≤mα成立。令各次比較的顯著性水準α=0.05/m，并規定P≤0.05/m時拒絕原假設[9]。在本研究中，針對牛的每條染色體分別制定各條染色體的顯著水平，以0.05分別除每條染色體上的SNPs數目，作為每條染色體水平的顯著性水平[10]。

表1 Case-control關聯分析中LASCS和SCS-SD劃分的描述性統計Table 1 Descriptive statistics of LASCS and SCS-SD for case-control association testing

圖1 用于Case-control分析的LASCS和SCS-SD的統計性描述Fig.1 Descriptive statistics of LASCS and SCS-SD for case-control association testing

表2 LASCS（SCS-SD）半個標準差/一個標準差質控后SNPs的分布和相鄰SNPs的平均距離Table 2 Distribution of SNPs on chromosomes after quality control and the average distances between adjacent SNPs of 0.5/ 1 SD of LASCS（SCS-SD）

2 結果

2.1 全基因組關聯分析

經過對基于LASCS和SCS-SD的乳房炎抗性進行全基因組關聯分析，得到了3種檢驗（RM檢驗、RW檢驗和RCHI檢驗）的全基因組SNPs位點的-log10 (Pvalue)結果，如圖2和圖3所示。通過對基于半個標準差SCS-SD的乳房炎抗性進行全基因組關聯分析發現一個全基因組水平顯著的SNP（Hapmap48573-BTA-104531,P=1.11E-06）位于X染色體上（圖3-B1）。

2.2 顯著SNPs及相關基因

通過Bonferroni校正得到每條染色體0.05顯著性水平的閾值（表3），共發現5個達到染色體水平顯著的SNPs。其中3個SNPs定位到X染色體上，其它2個SNPs分別定位到7和28號染色體上。結果發現，X染色體的顯著性SNPs（Hapmap48573-BTA-104531和Hapmap54175- rs29021817）位于IL1RAPL2基因內，7號染色體的顯著性SNP周圍存在與炎癥反應相關的基因（ILF3）。

圖2 基于LASCS的乳房炎抗性檢驗后的全基因組SNP位點的-log10 (P value)Fig. 2 Manhattan plots (-log10 (P value)) of genome-wide SNPs for mastitis resistance based on LASCS

圖3 基于SCS-SD的乳房炎抗性檢驗后的全基因組SNP位點的-log10 (P value)Fig.3 Manhattan plots (-log10 (P value)) of genome-wide SNPs for mastitis resistance based on SCS-SD

3 討論

本研究通過全基因組關聯分析，獲得5個在染色體水平與中國荷斯坦牛乳房炎抗性顯著關聯的SNPs，其中X染色體上2個顯著性SNPs（Hapmap48573-BTA-104531和Hapmap54175-rs29021817）都位于白介素1受體附屬蛋白2基因（Interleukin 1 receptor accessory protein-like-2，IL1RAPL2）內。7號染色體上的顯著SNP（BTA-78357- no-rs）在100 kp范圍內發現了與炎癥相關的白介素增強結合因子3基因（Interleukin enhancer binding factor 3，ILF3）。

IL1RAPL2基因已經被證實在大腦發育過程中有特殊的表達，并在神經功能的發育中起著關鍵作用。目前的研究尚未發現IL1RAPL2基因與特定奶牛疾病相關，但被認為是神經失調的候選基因[11]。ILF3基因是白介素家族中的一個跟炎癥反應相關的因子，其功能與抑制翻譯蛋白有關[12]。其他研究還發現，小鼠中的ILF3基因在劇毒的弗朗西斯菌中曝光4h后，其表達量會增加，此過程與NFAT（nuclear factor of activated T cells）信號通路有關。而當感染弗朗西斯菌急性發病時，很多出現表達增加的免疫相關基因參與了淋巴細胞激活過程，其中就包括NFAT信號通路中的基因[13]。目前關于這兩個基因的研究較少，尚未發現其與奶牛乳房炎抗性有關，但這兩個基因都與白介素有關，而白介素4、5、6、12、13、17、22、23等都參與了不同類型的炎癥反應，并發揮了重要的作用[14-16]，因此應該進一步對這兩個基因進行功能驗證，以作為奶牛乳房炎抗性的候選基因。

表3 牛基因組每條染色體所對應Bonferroni校正的顯著性標準Table 3 Significant P-value of Bonferroni correction on each bovine chromosome

表4 乳房炎抗性全基因組關聯分析顯著SNPs和相關基因Table 4 Significant SNPs and the nearest genes of genome-wide association with mastitis resistance

本研究基于WIJGA的研究[6]，構建了中國荷斯坦牛乳房炎抗性的統計量LASCS和SCS-SD。盡管LASCS和SCS-SD為連續性狀，但奶牛乳房炎抗性作為一個閾性狀，更適合使用Case-control進行全基因組關聯分析[17-18]。最近幾年里，關聯分析研究開始向Case和Control轉變，目的是為了增加檢驗效力[19]。為檢測到與奶牛乳房炎癥反應功能相關的基因，本研究采用Case-control方法研究奶牛乳房炎易感性及抗性，并對LASCS和SCS-SD進行了乳房炎易感性牛（Case）和抗性牛（Control）的劃分。通過全基因組關聯分析發現，顯著SNPs主要集中于X染色體上，其中包括全基因組水平顯著的SNP（Hapmap48573- BTA-104531）。而WIJGA將SCC主要定位到牛的4、6和18號染色體上，其關聯分析顯著SNP的閾值是7.94E-06[6]。通過對比其它研究發現，僅有Abdel-Shaft在德國荷斯坦牛的研究將與SCS有關顯著的SNP（rs41629005, Position: 30639394）定位到X染色體上[20]。雖然其它研究通過常規的單標記回歸分析將影響乳房炎抗性的SNPs定位到常染色體上，但結果多種多樣[21-23]，這也說明奶牛乳房炎抗性受微效多基因控制，常規方法不易定位到效應較大的分子標記。本研究首次通過LASCS和SCS-SD對中國荷斯坦牛乳房炎易感性牛和抗性牛進行劃分，然后進行Case-control關聯分析。本研究的檢驗效力要高于前期的研究，找到的相關基因也與炎癥反應關系緊密[7]。

本研究還通過線性混合模型（linear mixed model，MMA）對LASCS和SCS-SD進行全基因組關聯分析以驗證Case-control關聯分析的結果，共檢測到4個染色體水平顯著的SNPs，并都定位到X染色體上。其中BTA-28466- no-rs被MMA和Case-control方法同時定位到（表5）。對比之前的研究發現[7,24]，用兩種方法對同一性狀進行關聯分析，可以定位到相同的位點，但不同性狀的結果就各不相同。我們前期對體細胞評分育種值（somatic cell score estimated breeding value, SCSEBV）進行關聯分析，兩種方法共同定位到14號染色體上的SNP（ARS- BFGL-NGS-100480），其中線性模型檢驗效力較高（P= 1.24E-10）[24]。然而用Case-control方法對測定日記錄[7]、SCSEBV[24]、LASCS和SCS-SD的Case和Control進行關聯分析，結果不盡相同，可能與SCC受環境影響較大、通過SCC劃分的Case和Control關聯分析檢驗效力較低有關。因此對于奶牛乳房炎抗性性狀，還應直接對乳房炎牛和健康牛進行二分類的Case-control關聯分析。

表5 基于線性混合模型的全基因組關聯分析顯著SNPsTable 5 Significant SNPs of genome-wide association based on MMA method

由于乳房炎易感性及抗性性狀屬于功能性狀，遺傳力偏低，受環境效應影響較大，而3種原理不同的檢驗方法可能導致共同定位的SNP的檢驗效力大小不一。盡管本研究未檢出對兩種性狀（LASCS和SCS-SD）關聯分析均達到顯著的SNP，但X染色體上的SNP（Hapmap48573-BTA-104531）用3種檢測方法均發現其P值較低（表6）。進一步的研究擬增加SNPs標記密度及牛群數量，擴群驗證已獲得的乳房炎易感性及抗性相關分子標記及基因。此外還需對已發現的顯著SNPs及臨近基因進行生物學通路分析及功能驗證，為揭示奶牛乳房炎易感性及抗性的分子機制提供更多數據。

表6 對SNP（Hapmap48573-BTA-104531）用3種方法檢驗后的P值Table 6 The P values of SNP (Hapmap48573-BTA-104531)by three tests

4 結論

本研究利用Case-control方法對荷斯坦牛乳房炎易感性及抗性進行了全基因組關聯分析，共檢測5個顯著性單核苷酸突變及2個與炎癥反應密切相關的基因（IL1RAPL2和ILF3基因），并發現X染色體上的SNP（Hapmap48573-BTA-104531）在全基因組水平顯著（P＜1.14E-06）。本研究首次通過對泌乳期平均體細胞評分和測定日記錄體細胞評分標準差進行Case-control方法關聯分析，并找到了與炎癥反應相關的基因，為奶牛乳房炎抗性功能相關基因的研究提供數據支持。

[1] HALASA T, HUIJPS K, OSTERAS O, HOGEVEEN H. Economic effects of bovine mastitis and mastitis management: A review.Veterinary Quarterly, 2007, 29(1):18-31.

[2] 張沅，張勤，孫東曉. 奶牛分子育種技術研究. 北京：中國農業大學出版社, 2012.

ZHANG Y, ZHANG Q, SUN D X.The Technical Research on Molecular Breeding in Dairy Cattle.Beijing: China AgriculturalUniversity Press, 2012. (in Chinese)

[3] HERINGSTADB, KLEMETSDAL G, RUANE J. Selection for mastitis resistance in dairy cattle: a review with focus on the situation in the Nordic countries.Livestock Production Science, 2000, 64(2-3): 95-106.

[4] AXFORD R F E, BISHOP S C, NICHOLAS F W, OWEN J B.Breeding for Disease Resistance in Farm Animals. CABI Publishing, Oxon, 2000, 42:407-424.

[5] RUPP R, BOICHARD D. Relationship of early first lactation somatic cell count with risk of subsequent first clinical mastitis.Livestock Production Science, 2000, 62:169-180.

[6] WIJGA S, BASTIAANSEN J W M., WALL E, STRANDBERG E, DE HASS Y, GIBLIN L, BOVENHUIS H. Genomic association with somatic cell score in first-lactation Holstein cows.Journal of Dairy Science, 2012, 95(2):899-908.

[7] 王曉，解小莉，王勝，錢夢櫻，馬裴裴，張沅，孫東曉，劉劍鋒，丁向東，姜力，王雅春，張毅，張勝利，張勤，俞英. 中國荷斯坦牛乳房炎易感性及抗性的全基因組關聯分析. 畜牧獸醫學報, 2013, 44 (12):1907-1912.

WANG X, XIE X L, WANG S, QIAN M Y, MA P P, ZHANG Y, SUN D X, LIU J F, DING X D, JIANG L, WANG Y C, ZHANG Y, ZHANG S L, ZHANG Q, YU Y. Genome-wide Association Study for Mastitis Susceptibility and Resistance in Chinese Holsteins.Acta Veterinaria et Zootechnica Sinica, 2013, 44(12):1907-1912. (in Chinese)

[8] THORNTON T, MCPEEK M S. ROADTRIPS: Case-control association testing with partially or completely unknown population and pedigree structure.American Journal of Human Genetics, 2010, 86(2):172-184.

[9] BLAND J M, ALTMAN D G. Multiple significance tests: the Bonferroni method.British Medical Journal, 1995, 310(6973):170.

[10] 盧昕. 對豬部分免疫性狀的QTL定位及全基因組關聯分析[D]. 北京：中國農業大學，2010.

LU X. Quantitative trait loci (QTL) mapping and genome-wide association study for some immune traits in swine[D]. Beijing: China Agricultural University, 2010. (in Chinese)

[11] FERRANTE M I, GHIANI M, BULFONE A, FRANCO B.IL1RAPL2maps to Xq22 and is specifically expressed in the central nervous system.Gene, 2001, 275:217-221.

[12] LI H, SUN L, CHEN X, XIONG W, HU D, JIE S. Microvesicle microRNA profiles and functional roles between chronic hepatitis B and hepatocellular carcinoma.Gene, 2014, 16:315-321.

[13] WALTERS K, OLSUFKA R, KUESTNER R E, CHO H J, LI H, ZORNETZER G A, WANG K, SKERRETT S J, OZINSKY A.Francisella tularensissubsp.tularensisinduces a unique pulmonary inflammatory response: Role of bacterial gene expression in temporal regulation of host defense responses.PLoS ONE, 2013, 8(5):e62412.

[14] BRADDING P, FEATHER I H, WILSON S, BARDIN P G, HEUSSER C H, HOLGATE S T, HOWARTH P H. Immunolocalization of cytokines in the nasal mucosa of normal and perennial rhinitic subjects. The mast cell as a source of IL-4, IL-5, and IL-6 in human allergic mucosal inflammation.The Journal of Immunology, 1993, 151(7):3853-3865.

[15] MURPHY C A, LANGRISH C L, CHEN Y, BLUMENSCHEIN W, MCCLANAHAN T, KASTELEIN R A, SEDGWICK J D, CUA D J. Divergent Pro- and Antiinflammatory Roles for IL-23 and IL-12 in Joint Autoinmmune Inflammation.The Journal of Experimental Medicine, 2003, 198(12):1951-1957.

[16] ZHENG Y, DANILENKO D M, VALDEZ P, KASMAN I, EASTHAM-ANDERSON J, WU J, OUYANG W. Interleukin-22, a TH17 cytokine, mediates IL-23-induced dermal inflammation and acanthosis.NatureLetter, 2007, 445:648-651.

[17] DE HAAS Y, OUWELTJES W, TEN NAPEL J, WINDIG J J, DE JONG G. Alternative somatic cell count traits as mastitis indicators for genetic selection.Journal of Dairy Science, 2008, 91(6):2501-2511.

[18] HERINGSTAD B, REKAYA R, GLANOLA D, KLEMETSDAL G, WEIGEL K A. Genetic change for clinical mastitis in Norwegian cattle: a threshold model analysis.Journal of Dairy Science, 2003, 86(1):369-375.

[19] FERGUSON-SMITH A C, GREALLY J M, MARTIENSSEN R A.Epigenomics.Springer, Dordrecht, 2009.

[20] ABDEL-SHAFY H, BORTFELDT R H, REISSMANN M, BROCKMANN G A. Short communication: Validation of somatic cell score-associated loci identified in a genome-wide association study in German Holstein cattle.Journal of Dairy Science, 2014, 97(4): 2481-2486.

[21] MEREDITH B K, KEARNEY F J, FINLAY E K, BRADLEY D G, FAHEY A G, BERRY D P, LYNN D J. Genome-wide associations for milk production and somatic cell score in Holstein-Friesian cattle in Ireland.BMC Genetics, 2012, 13:21.

[22] SODELAND M, KENT M P, OLSEN H G, OPSAL M A, SVENDSEN M, SEHESTED E, HAYES B J, LIEN S. Quantitative trait loci for clinical mastitis on chromosomes 2, 6, 14 and 20 in Norwegian Red cattle.Animal Genetics, 2011, 42(5):457-465.

[23] SAHANA G, GULDBRANDTSEN B, THOMSEN B, LUND M S. Confirmation and fine-mapping of clinical mastitis and somatic cell score QTL in Nordic Holstein cattle.Animal Genetics, 2013, 44(6): 620-626.

[24] WANG X, MA P P, LIU J F, ZHANG Q, ZHANG Y, DING X D, JIANG L, WANG Y C, ZHANG Y, SUN D X, ZHANG S L, SU G S, YU Y. Genome-wide association study in Chinese Holstein cows reveal two candidate genes for somatic cell score as an indicator for mastitis susceptibility.BMC genetics, 2015, 16:111.

（責任編輯 林鑒非）

Genome-Wide Association Study on Mastitis Resistance Based on Somatic Cell Scores in Chinese Holstein Cows

WANG Xiao, ZHANG Qin, YU Ying
(College of Animal Science and Technology, China Agricultural University, Beijing 100193)

Chinese Holstein cows; mastitis resistance; log-transformed SCC; Case-control; GWAS

10.3864/j.issn.0578-1752.2017.04.015

2016-03-24；接受日期：2016-12-20

國家自然科學基金（31272420）、農業部奶業體系項目（CARS-37-04B）、十二五國家科技支持項目（2011BAD28B02）、“863”重大項目（2008AA101002）、教育部基本科研項目（2011JS006）、長江學者與創新團隊發展計劃（IRT1191）

聯系方式：王曉，E-mail：wangxiao880923@gmail.com。通信作者俞英，E-mail：yuying@cau.edu.cn