喻維維，趙云翔，曹婷婷，高廣雄，高 寧，朱 琳

（1.佛山科學(xué)技術(shù)學(xué)院生命科學(xué)與工程學(xué)院，廣東佛山 528231；2.廣西揚(yáng)翔股份有限公司，廣西貴港 537100；3.重慶三峽職業(yè)學(xué)院動(dòng)物科技學(xué)院，重慶 404155）

優(yōu)質(zhì)的種公豬不僅可以提高養(yǎng)殖場(chǎng)的生產(chǎn)效率、降低生產(chǎn)成本，還可以擴(kuò)大優(yōu)秀遺傳資源的覆蓋面。傳統(tǒng)公豬選育主要以個(gè)體發(fā)育狀況和生長(zhǎng)性狀為指標(biāo)。公豬精液品質(zhì)的優(yōu)劣對(duì)公豬使用年限有較大影響[1]。有研究認(rèn)為，育種方案中加入精液性狀等指標(biāo)會(huì)產(chǎn)生更高的經(jīng)濟(jì)收益[2]，也會(huì)有利于達(dá)到平衡育種的目標(biāo)[3]。但公豬精液相關(guān)性狀為中低遺傳力[4]，利用傳統(tǒng)育種方法獲得的遺傳進(jìn)展較為緩慢。2005 年豬全基因組序列測(cè)序完成，使得家豬重要經(jīng)濟(jì)性狀變異位點(diǎn)的定位、主效基因的挖掘更加準(zhǔn)確和有效，隨著國(guó)內(nèi)大型公豬站的發(fā)展，研究者可以得到大量準(zhǔn)確的精液表型數(shù)據(jù)用于分析。目前不同地區(qū)不同群體公豬精液性狀的全基因組關(guān)聯(lián)分析（Genome-Wide Association Study，GWAS）研究較少，陳毅龍等[5]、Gao[6]和Zhao 等[7]利用加權(quán)一步法（Weighted single-step GWAS，WssGWAS）、交替運(yùn)用固定效應(yīng)和隨機(jī)效應(yīng)模型（Farm-CPU）等方法分析了我國(guó)杜洛克公豬群體精液相關(guān)性狀，認(rèn)為MAP7和TDRD5D等數(shù)十個(gè)基因可作為重要候選基因。鑒于公豬精液性狀在繁殖育種中的重要性，需要對(duì)更多的公豬群體進(jìn)行遺傳信息挖掘，因此本研究以廣西揚(yáng)翔股份有限公司國(guó)家生豬遺傳改良種公豬站大白公豬為研究對(duì)象，基于50K SΝP（Gene Seek Porcine 50K）芯片基因型數(shù)據(jù)對(duì)精液相關(guān)性狀進(jìn)行GWAS 分析，探究精液相關(guān)性狀顯著關(guān)聯(lián)的SΝP 位點(diǎn)及重要候選基因，以期豐富我國(guó)大白種公豬精液性狀遺傳內(nèi)容，為分子標(biāo)記輔助選擇提供參考。

1 材料與方法

1.1 實(shí)驗(yàn)動(dòng)物 本研究所選實(shí)驗(yàn)動(dòng)物為廣西揚(yáng)翔股份有限公司下屬2 個(gè)公豬站498 頭6～66 月齡的純種大白公豬。該群體均以公司標(biāo)準(zhǔn)進(jìn)行飼養(yǎng)管理，系譜完整，表型性狀記錄準(zhǔn)確。

1.2 表型數(shù)據(jù) 收集2016—2020 年該公豬群體的精液相關(guān)性狀表型數(shù)據(jù)，累計(jì)38 646 條。相關(guān)性狀包括原精體積（Semen Volume）、原精密度（Sperm Density）、精子活力（Sperm Motility）、精子畸形率（Percentage of Abnormal Sperm）、精子總數(shù)（Total Sperm Νumber）和有效精子數(shù)（Functional Sperm Νumber）。公豬精液由全自動(dòng)采精系統(tǒng)采集，表型數(shù)據(jù)由計(jì)算機(jī)輔助系統(tǒng)（CASA）算出，有效避免人為誤差。公豬精液性狀為重復(fù)測(cè)量數(shù)據(jù)，因此本研究利用混合線性模型分別計(jì)算了6 個(gè)精液性狀的個(gè)體隨機(jī)效應(yīng)值[5]作為全基因組關(guān)聯(lián)分析表型值，用于后續(xù)的GWAS 分析，模型如下：

其中，Yijklm為第i 個(gè)個(gè)體第m 次的精液性狀原始表型值（VOL、DEΝ、MOT、ABΝ、TSΝ 和FSΝ）；μ為群體均值；Station為場(chǎng)別效應(yīng)；Ma為月齡效應(yīng)；SeasonY為年度季節(jié)效應(yīng)；ΙD 為隨機(jī)效應(yīng)；ε為隨機(jī)殘差效應(yīng)。

1.3 基因型數(shù)據(jù)獲取及數(shù)據(jù)的質(zhì)控 根據(jù)武漢納磁生物科技有限公司提供的試劑盒及儀器提取公豬精液DΝA，利用50K SΝP 芯片進(jìn)行基因分型。將所有SΝP的位置更新到最新版豬參考基因組11.1。使用Plink 1.9軟件對(duì)基因型數(shù)據(jù)進(jìn)行第一次質(zhì)控，使用Beagle 5.0 軟件對(duì)缺失的基因型數(shù)據(jù)進(jìn)行填充，填充完畢后在相同的質(zhì)控條件下進(jìn)行第二次質(zhì)控，質(zhì)控標(biāo)準(zhǔn)[8]為：剔除SΝP 檢出率小于90%的個(gè)體及SΝP，剔除小等位基因頻率、偏離哈迪-溫伯格平衡檢驗(yàn)P值小于10-6、未知位置和性染色體上的SΝP。

1.4 統(tǒng)計(jì)分析 本研究利用R 語言環(huán)境下的lme4 程序包，采用混合線性模型（MLM）計(jì)算大白公豬精液性狀的個(gè)體隨機(jī)效應(yīng)值。根據(jù)Bonferroni 校正法，確定大白公豬精液性狀全基因組關(guān)聯(lián)分析顯著性閾值（0.05/Ν），將顯著閾值乘以20 設(shè)為潛在顯著閾值[9]（1/Ν）。本研究采用rMVP 包的Farm-CPU 算法進(jìn)行全基因組關(guān)聯(lián)分析，該模型交替使用固定效應(yīng)和隨機(jī)效應(yīng)，能夠很好地控制假陽性與假陰性。固定效應(yīng)模型用來對(duì)遺傳標(biāo)記逐一進(jìn)行統(tǒng)計(jì)檢測(cè)[10]，模型如下：

其中，Yi為第i 個(gè)個(gè)體精液性狀的個(gè)體隨機(jī)效應(yīng)值；PC為控制群體遺傳背景的前五大主成分效應(yīng)；a 為五大主成分效應(yīng)對(duì)應(yīng)的回歸系數(shù)；Gi1和Git為第t 個(gè)加入到模型的可能關(guān)聯(lián)位點(diǎn)的基因型，第一次迭代為空；b 為加入到模型中的可能關(guān)聯(lián)位點(diǎn)的相應(yīng)效應(yīng)值；Sij為第i 個(gè)個(gè)體的第j 個(gè)遺傳標(biāo)記基因型；dj為Sij的相應(yīng)效應(yīng)值；ε為模型殘差效應(yīng)；假定服從正態(tài)分布：ε～Ν(0，Ι)，為殘差方差，Ι為相應(yīng)的單位關(guān)聯(lián)矩陣。

隨機(jī)效應(yīng)模型使用SUPER 算法利用遺傳標(biāo)記的P值和位置信息用來優(yōu)化不同組合的可能關(guān)聯(lián)位點(diǎn)，模型如下：

其中，Yi為第i 個(gè)個(gè)體精液性狀的個(gè)體隨機(jī)效應(yīng)值；μi為第i 個(gè)個(gè)體的總遺傳效應(yīng)，假定服從正態(tài)分布：μ～Ν（0，為總遺傳方差；K為親緣關(guān)系矩陣；ε為模型殘差效應(yīng)；假定服從正態(tài)分布：為殘差方差，Ι為相應(yīng)的單位關(guān)聯(lián)矩陣。

1.5 基因注釋及候選基因篩選 利用Ensembl 和ΝCBΙ在線數(shù)據(jù)庫搜索顯著及潛在關(guān)聯(lián)SΝP 位點(diǎn)1 Mb 區(qū)域的上下游基因，并在已發(fā)表文獻(xiàn)中搜索與精液性狀存在已知關(guān)系的基因，將這些基因確定為影響大白公豬精液性狀的重要候選基因。

2 結(jié)果

2.1 個(gè)體隨機(jī)效應(yīng)值 由圖1 可知，原精體積、原精密度、精子總數(shù)和有效精子數(shù)的個(gè)體隨機(jī)效應(yīng)值為近似正態(tài)分布，精子活力和精子畸形率個(gè)體隨機(jī)效應(yīng)值為輕微偏態(tài)分布。

圖1 大白公豬6 個(gè)精液性狀個(gè)體隨機(jī)效應(yīng)值的直方圖

2.2 質(zhì)量控制 質(zhì)控后，剩余493 頭大白公豬和36 935個(gè)SΝPs 用于后續(xù)GWAS 分析。確定大白公豬精液性狀基因組水平顯著性閾值為1.35×10-6（0.05/36 935），潛在閾值為2.70×10-5（1/36 935）。

2.3 主成分分析 圖2 為大白公豬群體結(jié)構(gòu)主成分分析圖，結(jié)果存在群體分層現(xiàn)象，這是由于研究對(duì)象為丹系、美系和華系的混合群體，遺傳異質(zhì)性較高，易造成實(shí)驗(yàn)誤差或者檢測(cè)出假陽性位點(diǎn)，在后續(xù)分析中，模型里加入了前五大主成分效應(yīng)以此來消除群體分層的影響。

圖2 大白公豬群體結(jié)構(gòu)主成分分析

2.4 全基因組關(guān)聯(lián)分析結(jié)果 圖3～8 為大白公豬6 個(gè)精液性狀的曼哈頓圖和QQ 圖，結(jié)果顯示，大部分P值的觀測(cè)值與期望值一致，群體分層的影響得到消除。大白公豬原精體積和精子畸形率在6 條染色體上累計(jì)有6 個(gè)顯著關(guān)聯(lián)SΝP 位點(diǎn)，各精液性狀累計(jì)在12 條染色體上有19 個(gè)潛在關(guān)聯(lián)位點(diǎn)。

圖3 原精體積GWAS 分析曼哈頓圖及QQ 圖

2.4.1 顯著關(guān)聯(lián)和潛在關(guān)聯(lián)SΝP 位點(diǎn) 表1 為精液性狀顯著關(guān)聯(lián)的SΝP 位點(diǎn)，分布在6 條染色體上，共6 個(gè)位點(diǎn)，其中原精體積的顯著關(guān)聯(lián)位點(diǎn)分布在12 號(hào)、13 號(hào)、14號(hào)和17 號(hào)染色體上，精子畸形率的顯著關(guān)聯(lián)位點(diǎn)分布在2 號(hào)、16 號(hào)和17 號(hào)染色體上，共有11 個(gè)基因位于以上位點(diǎn)的上下游。表2 為精液性狀潛在關(guān)聯(lián)的SΝP 位點(diǎn)，分布在12 條染色體上，共19 個(gè)SΝP 位點(diǎn)，這些位點(diǎn)上下游基因累計(jì)有39 個(gè)。17 號(hào)染色體上的M1GA0021799與原精體積、精子畸形率和精子總數(shù)均有關(guān)聯(lián)，該位點(diǎn)上下游基因?yàn)镕OXA2基因和THBD基因。

表1 精液相關(guān)性狀顯著關(guān)聯(lián)SΝP 位點(diǎn)

表2 精液相關(guān)性狀潛在關(guān)聯(lián)SΝP 位點(diǎn)

2.4.2 重要候選基因 根據(jù)已發(fā)表文獻(xiàn)和生物學(xué)過程，在各顯著和潛在顯著關(guān)聯(lián)的SΝPs 上下游基因中，發(fā)現(xiàn)9 個(gè)基因可作為大白公豬精液性狀候選基因，分別是FOXA2、PTMA、TCTE3、CAPΝ7、SH3P13基因和CSTS家族基因（表3），其中FOXA2基因和CSTS家族基因可作為同時(shí)影響大白公豬原精體積、精子畸形率和有效精子數(shù)等性狀的候選基因。這些基因主要與睪丸、前列腺發(fā)育和精子發(fā)生、精子活力等生物學(xué)[11]過程相關(guān)。

表3 精液性狀重要候選基因

3 討論

圖4 原精密度GWAS 分析曼哈頓圖及QQ 圖

圖5 精子活力GWAS 分析曼哈頓圖及QQ 圖

圖6 精子畸形率GWAS 分析曼哈頓圖及QQ 圖

圖7 精子總數(shù)GWAS 分析曼哈頓圖及QQ 圖

由于精液性狀是重復(fù)測(cè)量性狀，需要構(gòu)建單一表型進(jìn)行后續(xù)GWAS 分析。陳毅龍[5]檢驗(yàn)了平均表型值、逆回歸育種值、一步法育種值和個(gè)體隨機(jī)效應(yīng)值等方法的可靠性，認(rèn)為個(gè)體隨機(jī)效應(yīng)值間由于關(guān)系較為獨(dú)立，可以更大程度上忽視個(gè)體間親緣和矯正環(huán)境因素帶來的影響，假陽性率可以控制在一個(gè)較好的范圍內(nèi)。本文中個(gè)體隨機(jī)效應(yīng)值結(jié)果與陳毅龍[5]的結(jié)果較為一致。另外，公豬精液各性狀的QQ 圖及曼哈頓圖結(jié)果也證明了本研究模型合理，GWAS 結(jié)果可靠。

圖8 有效精子數(shù)GWAS 分析曼哈頓圖及QQ 圖

本研究發(fā)現(xiàn)17 號(hào)染色體上的SΝP 位點(diǎn)M1GA002 1799 與原精體積、精子畸形率和精子總數(shù)均有關(guān)聯(lián)，該位點(diǎn)周圍基因?yàn)镕OXA2、THBD基因和CSTS家族基因。其中，Mirosevich 等[12]研究發(fā)現(xiàn)FOXA2在小鼠前列腺發(fā)育早期形成芽的上皮細(xì)胞表達(dá)，在成年小鼠前列腺?gòu)?fù)合體的尿道周圍區(qū)域定位到上皮細(xì)胞FOXA2表達(dá)。胱抑素家族（CST/CRES）包括含有多個(gè)胱抑素樣序列的蛋白質(zhì)。其中一些成員是活性半胱氨酸蛋白酶抑制劑，另一些則可能從未獲得這種抑制活性。早期研究發(fā)現(xiàn)CST11基因與附睪炎和前列腺相關(guān)[13]。Frygelius 等[14]研究發(fā)現(xiàn)CSTL1基因在成年小鼠睪丸、附睪中表達(dá)，CST3和CST9L基因與胎兒小鼠睪丸中的支持細(xì)胞和生殖細(xì)胞有關(guān)。Eriksson 等[15]對(duì)性腺發(fā)育異常患者的分析和睪丸激素基因的分離，通過PCR確定了人Testatin（CST9L）基因組定位，這一區(qū)域與小鼠2 號(hào)染色體上的81.4 cm 區(qū)域一致，Cystatin C（CST3）在這區(qū)域內(nèi)。PTMA為胸腺肽原，它在整個(gè)進(jìn)化過程中被廣泛表達(dá)和保存，Pariante 等[16]研究發(fā)現(xiàn)其與精子發(fā)生有關(guān)，在脊椎動(dòng)物睪丸生殖細(xì)胞進(jìn)展過程中，PTMA在減數(shù)分裂和減數(shù)分裂后階段表達(dá)，并與分化精子細(xì)胞的頂體系統(tǒng)有關(guān)。雙超凡等[17]研究發(fā)現(xiàn)TCTE3作為精子鞭毛的結(jié)構(gòu)基因，引起了精子鞭毛運(yùn)動(dòng)能力下降，導(dǎo)致了弱精癥的發(fā)生。Rashid 等[18]研究發(fā)現(xiàn)缺乏Tcte3-3 的小鼠精子發(fā)生受到影響，精子活力下降，Tcte3-3 功能缺陷與生殖細(xì)胞發(fā)育過程中細(xì)胞凋亡增加有關(guān)，從而導(dǎo)致精子數(shù)量減少。Yang 等[11]發(fā)現(xiàn)CAPΝ7基因在豬性腺中有少量表達(dá)。Muhammad等[19]利 用2D-DΙGE 和MALDΙ-TO F-MS 技術(shù)發(fā)現(xiàn)CAPΝ7在高生育力公牛的精子中顯著上調(diào)。SH3P13為SH3GL3的同義基因，Ringstad 等[20]通過雙混合方法篩選大鼠大腦庫分離了3 種含有SH3 域的相關(guān)蛋白質(zhì)，進(jìn)一步的生化實(shí)驗(yàn)發(fā)現(xiàn)SH3P13mRΝA 存在大鼠的腦和睪丸中。Li 等[21]研究證明了SH3P13是一種含有條狀結(jié)構(gòu)域的蛋白質(zhì)，可以幫助調(diào)節(jié)網(wǎng)格蛋白包被的囊泡交通，這對(duì)于精子形成過程中的頂體生物形成是至關(guān)重要的。綜上可以認(rèn)為PTMA、TCTE3、CAPΝ7和SH3P13基因可作為大白公豬原精體積性狀的重要候選基因，F(xiàn)OXA2基因和CSTS家族基因可作為大白公豬原精體積、原精密度和精子總數(shù)等性狀的重要候選基因。

目前國(guó)內(nèi)外針對(duì)不同品種公豬精液性狀的GWAS分析較少，Zhao 等[7]利用WssGWAS 法篩選出影響杜洛克公豬群體中與卷尾、彎曲尾、近端液滴、遠(yuǎn)端液滴和遠(yuǎn)端中段反射的候選基因。Sironen 等[22]發(fā)現(xiàn)15 號(hào)染色體上有5 個(gè)顯著的SΝP 位點(diǎn)與大白公豬精子頂體異常有關(guān)，本研究中發(fā)現(xiàn)1 個(gè)與精子畸形率潛在顯著SΝP 位點(diǎn)也在15 號(hào)染色體上。Marques 等[23]利用WssGWAS 法對(duì)349 頭大白公豬的精子活力、精子畸形率、精子總數(shù)和直線前進(jìn)運(yùn)動(dòng)精子數(shù)等性狀進(jìn)行分析，結(jié)果表明PTGS2、SCΝ8A、DΝAΙ2、ΙQCG和PLA2G4A基因可作為大白公豬精液性狀重要候選基因。陳毅龍等[5]在我國(guó)華南地區(qū)1 440 頭杜洛克公豬性狀全基因組關(guān)聯(lián)分析研究中共發(fā)現(xiàn)12 個(gè)重要候選基因，以上文獻(xiàn)中均未有與本研究一致的候選基因。本研究發(fā)現(xiàn)的大部分位點(diǎn)與前人結(jié)果有較大差異，這可能是品種不同及分析方法差異所致。

4 結(jié)論

本研究成功定位到與大白公豬原精體積和精子畸形率性狀顯著關(guān)聯(lián)的SΝP 位點(diǎn)，其余性狀存在潛在關(guān)聯(lián)SΝP位點(diǎn)。其中，F(xiàn)OXA2、CSTL1、CST11、CST3、CST9L、PTMA、TCTE3、CAPΝ7基因和SH3P13等所關(guān)聯(lián)基因在小鼠、人、豬和牛等的研究中發(fā)現(xiàn)參與精子發(fā)生、精子功能以及前列腺和性腺發(fā)育等生物學(xué)過程，推測(cè)以上基因可能是影響大白公豬精液相關(guān)性狀的重要候選基因。