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

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

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

1 材料與方法

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

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

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

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

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

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

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

其中，Yi為第i 個個體精液性狀的個體隨機效應(yīng)值；μi為第i 個個體的總遺傳效應(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 位點1 Mb 區(qū)域的上下游基因，并在已發(fā)表文獻中搜索與精液性狀存在已知關(guān)系的基因，將這些基因確定為影響大白公豬精液性狀的重要候選基因。

2 結(jié)果

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

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

2.2 質(zhì)量控制 質(zhì)控后，剩余493 頭大白公豬和36 935個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)象，這是由于研究對象為丹系、美系和華系的混合群體，遺傳異質(zhì)性較高，易造成實驗誤差或者檢測出假陽性位點，在后續(xù)分析中，模型里加入了前五大主成分效應(yīng)以此來消除群體分層的影響。

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

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

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

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

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

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

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

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

3 討論

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

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

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

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

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

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

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

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

4 結(jié)論

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