基于体细胞评分的中国荷斯坦牛乳房炎抗性全基因组关联分析

王晓，张勤，俞英

（中国农业大学动物科技学院，北京 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