盧 鑫，周靖航，楊朝云，張夢華，葉連萌，李叔臻，黃錫霞，馬 云，王興平，史遠(yuǎn)剛，*

(1.寧夏大學(xué) 農(nóng)學(xué)院，寧夏 銀川 750021； 2.新疆農(nóng)業(yè)大學(xué) 動物科學(xué)學(xué)院，新疆 烏魯木齊 830052)

新疆褐牛是我國自主培育的乳肉兼用型地方品種，具有高乳脂、高乳蛋白的乳用特征和肉質(zhì)優(yōu)良、風(fēng)味鮮美的肉用性能，還具有超強(qiáng)的抗逆性和卓越的適應(yīng)能力，在新疆畜牧業(yè)發(fā)展中有著舉足輕重的地位[1-2]。本研究以新疆褐牛產(chǎn)奶和繁殖性狀相關(guān)候選基因?yàn)檠芯繉ο?，對其進(jìn)行功能注釋和富集分析，旨在為后期的基因功能驗(yàn)證、基因組選擇和分子育種工作奠定理論基礎(chǔ)。目前關(guān)于新疆褐牛分子育種相關(guān)基因的研究已有許多報(bào)道，李娜[3]發(fā)現(xiàn)LEP、FAS、LPL、FTO等4種脂代謝調(diào)控基因的表達(dá)在新疆褐牛中具有組織特異性，劉麗元等[4]通過DNA池重測序等方法得到與產(chǎn)奶性狀和體細(xì)胞評分相關(guān)的單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)位點(diǎn)分別位于HAL、ZNF66、GPIHBP1和PDE9A等8個(gè)基因中，韓麗云[5]對MYF5、PPARγ和AGPAT基因的遺傳多態(tài)性與生長性狀進(jìn)行了關(guān)聯(lián)分析，林嘉鵬等[6]研究POU1F1基因第6外顯子的多態(tài)性與產(chǎn)奶量之間的關(guān)系，發(fā)現(xiàn)在新疆褐牛群體中B等位基因?qū)Ξa(chǎn)奶量具有正效應(yīng)，以上研究均為新疆褐牛產(chǎn)奶性狀和肉質(zhì)性狀的分子標(biāo)記輔助選育提供了理論依據(jù)和參考意見。

實(shí)施標(biāo)記輔助選擇可以加快奶牛分子育種遺傳進(jìn)展，提高育種的準(zhǔn)確性，加快品種的改良和新品系的選育。對控制畜禽數(shù)量性狀基因座和主效基因的精準(zhǔn)定位是分子育種的前期基礎(chǔ)。隨著?；蚪M測序計(jì)劃的完成和商業(yè)化高密度SNP芯片及分型技術(shù)的發(fā)展，全基因組關(guān)聯(lián)分析逐漸取代常規(guī)的數(shù)量性狀座位(quantitative trait locus，QTL)定位方法成為鑒定畜禽數(shù)量性狀功能基因的最有效策略[7]。目前關(guān)于牛的生長性狀[8-9]、繁殖性狀[10-11]、產(chǎn)奶性狀[12-14]全基因組關(guān)聯(lián)分析的報(bào)道已經(jīng)有很多，這些研究為目標(biāo)性狀形成分子機(jī)理的深入研究提供了重要線索。產(chǎn)奶性狀是奶牛最重要的經(jīng)濟(jì)性狀，而奶牛的繁殖效率直接影響生產(chǎn)的經(jīng)濟(jì)效益，二者均是由多個(gè)部分組成的復(fù)雜的生理過程，近些年來國內(nèi)外已逐漸將繁殖性狀納入平衡育種的考量因素[15]。

通過全基因組關(guān)聯(lián)分析(genome-wide association studies， GWAS)來探尋影響復(fù)雜性狀的相關(guān)基因已經(jīng)廣泛應(yīng)用開來，本課題組前期研究中采用150K牛基因芯片設(shè)計(jì)新疆褐牛資源群體，對奶牛的產(chǎn)奶性狀乳脂量(fat yield，F(xiàn)Y)、乳蛋白量(protein yield，PY)、乳脂率(fat percentage，F(xiàn)P)、乳蛋白率(protein percentage，PP)、產(chǎn)奶量(milk yield，MY)、體細(xì)胞評分(somatic cell score，SCS)及青年牛的初配日齡(age at first service，AFS)、初產(chǎn)日齡(age at first calving，AFC)、經(jīng)產(chǎn)牛的妊娠長度(gestation length，GL)、產(chǎn)犢間隔(calving interval，CI)進(jìn)行了全基因組關(guān)聯(lián)分析，共檢測到57個(gè)SNPs與其基因組水平顯著相關(guān)。本研究旨在基于前期新疆褐牛產(chǎn)奶性狀GWAS結(jié)果，進(jìn)一步對產(chǎn)奶性狀和繁殖性狀功能基因進(jìn)行位置候選鑒定及功能注釋，為后期產(chǎn)奶性狀和繁殖性狀的主效基因鑒定及后續(xù)基因功能驗(yàn)證提供研究基礎(chǔ)。

1 材料與方法

1.1 GWAS篩選候選基因

本課題組前期通過家系和表型記錄，從新疆褐牛4個(gè)核心群中挑選出403頭母牛，利用牛150K基因芯片進(jìn)行基因分型，每個(gè)樣本平均檢出139 376個(gè)標(biāo)記。對表型和基因型原始數(shù)據(jù)進(jìn)行質(zhì)量控制后，采用動物模型DMU軟件估計(jì)育種值，用主成分分析方法檢測種群結(jié)構(gòu)，最后用混合線性模型Farm CPU軟件進(jìn)行全基因組關(guān)聯(lián)分析，共計(jì)得到57個(gè)SNP位點(diǎn)分別與產(chǎn)奶性狀和繁殖性狀顯著相關(guān)。

1.2 功能基因組信息學(xué)分析

將顯著SNPs位點(diǎn)與NCBI網(wǎng)站數(shù)據(jù)庫中牛基因組序列數(shù)據(jù)庫(Bos taurus UMD 3.1.1)進(jìn)行比對，根據(jù)SNPs的物理位置來推斷其所在或臨近基因，再根據(jù)NCBI、Ensemble、Panther(http://www.pantherdb.org/)、GeneCards(https://www.genecards.org/)等數(shù)據(jù)庫中牛、人和小鼠對應(yīng)基因編碼產(chǎn)物的生理生化功能，對候選基因進(jìn)行功能和通路富集分析，結(jié)合相應(yīng)文獻(xiàn)報(bào)道，初步確定基因的基本生物學(xué)功能。

2 結(jié)果與分析

2.1 候選基因的鑒定

前期GWAS研究檢測分析得到共計(jì)57個(gè)顯著SNPs位點(diǎn)，分別分散在牛的30個(gè)染色體上。通過SNP位點(diǎn)的物理位置臨近的相關(guān)基因初步得到55個(gè)候選基因。根據(jù)這57個(gè)顯著SNPs位點(diǎn)的物理位置得出有29個(gè)位點(diǎn)位于基因側(cè)翼區(qū)域，其中與相鄰基因最近的距離是882 bp，有28個(gè)位點(diǎn)位于基因的內(nèi)含子區(qū)域。產(chǎn)奶性狀中與乳脂量相關(guān)的2個(gè)SNP位點(diǎn)均位于CNIH3基因；繁殖性狀中EPRS基因中有2個(gè)SNP位點(diǎn)，分別于初產(chǎn)日齡和初配日齡相關(guān)，結(jié)果見表1和表2。

2.2 候選基因功能注釋

通過功能基因組生物信息學(xué)分析，對上述55個(gè)基因進(jìn)行功能注釋。研究結(jié)果發(fā)現(xiàn)，14個(gè)未命名基因中有7個(gè)是根據(jù)序列推測出的未命名的假定基因，有7個(gè)基因編碼長鏈非編碼RNA，但其中有3個(gè)基因由于目前沒有后續(xù)的功能注釋已被NCI召回，LOC101903067、LOC104990303、LOC104969301、LOC104970790這4個(gè)基因可以編碼長鏈非編碼RNA。將被召回的3個(gè)未命名基因、在牛中未發(fā)現(xiàn)的ZDHHC17和RSPH10B基因和編碼tRNA(TRNAF-GAA)基因剔除，共得到49個(gè)基因的染色體定位及功能注釋信息。

表1 產(chǎn)奶性狀相關(guān)SNP位點(diǎn)信息

Table1SNP locus information and candidate genes related to milk production traits

性狀TraitSNP位點(diǎn)名稱SNP Name染色體號Chromosome No.位置Position與候選基因的距離Distance/bp候選基因Candidate geneFPBovineHD0100012734144652580內(nèi)部 WithinEFHBBTB-01890990426186690114212SNX13BovineHD0500001661561941766257824ZDHHC17ARS-BFGL-NGS-82026771141471882NIPAL4Hapmap34648-BES8_Contig530_7431032715286內(nèi)部 WithinMEIS2ARS-BFGL-BAC-61881838850678內(nèi)部 WithinLOC101903067BovineHD210001208521424101884827GPR33BovineHD24000079162429095464內(nèi)部 WithinCDH2FYBovineHD16000079771628501351內(nèi)部 WithinCNIH3BovineHD16000079841628528105內(nèi)部 WithinCNIH3PPHapmap54547-rs290121984110776624120776CNTNAP2BovineHD0500021288574915500內(nèi)部 WithinLOC100847379BTA-103663-no-rs16212438638464LOC538060BovineHD23000145342349827689662TRNAF-GAAPYBovineHD01000367531129629193內(nèi)部 WithinCLSTN2BovineHD0100040938114256936086982LOC104970303BovineHD03000000633343817內(nèi)部 WithinSFT2D2BTA-88698-no-rs74692971041250FSTL4BTB-01731924775830763內(nèi)部 WithinGABRG2ARS-BFGL-NGS-80161281302573內(nèi)部 WithinTMTC4BovineHD13000225391377860915內(nèi)部 WithinARFGEF2BovineHD2000015131205526823381531LOC784462BTA-112312-no-rs2069638890161115IRX1BovineHD25000107402538532556內(nèi)部 WithinRSPH10BMYBovineHD0600014209651560184內(nèi)部 WithinPCDH7BovineHD07000001587778150內(nèi)部 WithinGFPT2Hapmap45084-BTA-1137462719594127121781LOC104976064SCSBovineHD050001329654629133328240DYRK2BovineHD05000332915115103698內(nèi)部 WithinSULT4A1BovineHD0800003484810705865內(nèi)部 WithinSCARA5BovineHD0800007286824250348120194LOC104969301ARS-USDA-AGIL-chr9-30332222-000790930332222內(nèi)部 WithinTBC1D32BovineHD10000093731028618805內(nèi)部 WithinAVENBovineHD110003015611103632442內(nèi)部 WithinNACC2ARS-BFGL-BAC-13745125770736259530LOC782305ARS-BFGL-NGS-1075281525181966內(nèi)部 WithinLOC782610ARS-BFGL-NGS-433022237402411內(nèi)部 WithinPRICKLE2BovineHD22000122612242292699158098FHITBovineHD24000135202448435771163815ZBTB7CARS-BFGL-NGS-11953129460628955539GSTP1

表2 繁殖性狀相關(guān)SNP位點(diǎn)信息及候選基因

Table2SNP locus information and candidate genes related to reproductive traits

性狀TraitSNP位點(diǎn)名稱SNP Name染色體號Chromosome No.位置Position與候選基因的距離Distance/bp候選基因Candidate geneAFCHapmap43251-BTA-376121573304557351578LOC107133188BovineHD16000066911624235446內(nèi)部 WithinEPRSAFSARS-BFGL-NGS-998212676488726427LOC104971350BovineHD030003523731204966613195KIF1ABovineHD1400016327145878179972553RSPO2BovineHD16000066911624235446內(nèi)部 WithinEPRSGLBovineHD1400021729147746414076045LOC786994ARS-BFGL-NGS-24511110588483015267LOC785220ARS-USMARC-5281734752485內(nèi)部 WithinSPRY1BovineHD180001791418621244873301CACNG6ARS-BFGL-NGS-1077982210269204內(nèi)部 WithinSTACCIBovineHD25000034622512378774145129SHISA9BovineHD30000296503010794131644746USP9XBovineHD01000027411874039629865LOC104970790ARS-BFGL-NGS-14087166721297內(nèi)部 WithinHCLS1BovineHD190000200719755725034950ANKFN1BTB-01898603203957663761382RAI14

2.2.1 乳脂率候選基因

乳脂率候選基因包括EFHB、SNX13、NIPAL4、MEIS2、LOC101903067、GPR33、CDH2等7個(gè)基因。EFHB基因有7個(gè)可變剪接數(shù)，編碼EF手性結(jié)構(gòu)家族蛋白，可與鈣離子結(jié)合；SNX13基因有13個(gè)可變剪接數(shù)，能編碼分選蛋白SNX家族的PHOX結(jié)構(gòu)域和G蛋白信號轉(zhuǎn)導(dǎo)調(diào)節(jié)因子RGS家族的RGS結(jié)構(gòu)域，前者是磷脂酰肌結(jié)合結(jié)構(gòu)域，作為異源三聚體G蛋白Gα亞基的GTP酶激活蛋白的調(diào)節(jié)分子，后者參與細(xì)胞內(nèi)轉(zhuǎn)運(yùn)；NIPAL4基因有5個(gè)可變剪接數(shù)，可編碼膜受體蛋白，與鎂離子結(jié)合；MEIS2基因有27個(gè)可變剪接數(shù)，可編碼TALE家族同源異形蛋白，是高度保守的轉(zhuǎn)錄調(diào)節(jié)因子；LOC101903067基因有1個(gè)可變剪接數(shù)，可編碼長鏈非編碼RNA；GPR33基因有1個(gè)可變剪接數(shù)，被鑒定為孤獨(dú)基因趨化因子GPCR(G蛋白偶聯(lián)受體)的假基因，參與信號轉(zhuǎn)導(dǎo)；CDH2基因有5個(gè)可變剪接數(shù)，可編碼鈣黏蛋白，與鈣離子結(jié)合，參與細(xì)胞骨架蛋白合成，其中一個(gè)轉(zhuǎn)錄本可以編碼細(xì)胞黏附分子和糖蛋白，促進(jìn)神經(jīng)發(fā)育和骨骼形成。

2.2.2 乳脂量候選基因

乳脂量候選基因?yàn)镃NIH3基因，有14個(gè)可變剪接數(shù)，可編碼Cornichon家族AMPA受體蛋白，參與蛋白的代謝，協(xié)助蛋白轉(zhuǎn)運(yùn)到高爾基體并參與隨后的修飾。

2.2.3 乳蛋白率候選基因

乳蛋白率候選基因包括CNTNAP2、LOC100847379、LOC538060等3個(gè)。CNTNAP2基因有25個(gè)可變剪接數(shù)，可編碼接觸蛋白相關(guān)蛋白，含有表皮生長因子重復(fù)序列和層粘連蛋白G結(jié)構(gòu)域，可能有PDZ結(jié)合位點(diǎn)，在神經(jīng)系統(tǒng)發(fā)育中介導(dǎo)神經(jīng)元與膠質(zhì)相互作用，參與分化軸突內(nèi)鉀通道的定位；LOC100847379基因有1個(gè)可變剪接數(shù)，是編碼載脂蛋白L2類的假定基因；LOC538060基因有1個(gè)可變剪接數(shù)，可編碼二氫二醇脫氫酶1和20-α-(3-α)-羥基類固醇脫氫酶。

2.2.4 乳蛋白量候選基因

乳蛋白量候選基因包括CLSTN2、LOC104970303、SFT2D2、FSTL4、GABRG2、TMTC4、ARFGEF2、LOC784462、IRX1這9個(gè)。CLSTN2基因有2個(gè)可變剪接數(shù)，編碼鈣同線蛋白，可與鈣離子結(jié)合，調(diào)節(jié)鈣介導(dǎo)的突觸后信號；LOC104970303基因有1個(gè)可變剪接數(shù)，可編碼長鏈非編碼RNA；SFT2D2基因有4個(gè)可變剪接數(shù)，可編碼STF結(jié)構(gòu)域蛋白2，參與胞內(nèi)小泡介導(dǎo)途徑，將其逆向轉(zhuǎn)運(yùn)至高爾基復(fù)合體；FSTL4基因有7個(gè)可變剪接數(shù)，編碼卵泡抑素蛋白，可與鈣離子結(jié)合，對腦源性神經(jīng)營養(yǎng)因子受體信號通路有負(fù)調(diào)節(jié)作用；GABRG2基因有29個(gè)可變剪接數(shù)，編碼γ-氨基丁酸(GABA)受體蛋白，可與氯離子結(jié)合，是配體門控離子通道Cys環(huán)家族的成員，負(fù)責(zé)協(xié)調(diào)GABA的抑制效應(yīng)，參與葡萄糖、膽鹽和有機(jī)酸、金屬離子和胺類化合物的運(yùn)輸；TMTC4基因有11個(gè)可變剪接數(shù)，編碼含有N-端跨膜結(jié)構(gòu)域和C-端四肽重復(fù)序列結(jié)構(gòu)域的轉(zhuǎn)運(yùn)蛋白家族，介導(dǎo)了突觸囊泡融合、蛋白質(zhì)折疊和蛋白質(zhì)移位中的蛋白質(zhì)相互作用；ARFGEF2基因有2個(gè)可變剪接數(shù)，編碼核糖基化因子鳥嘌呤核苷酸交換因子2，參與胞內(nèi)囊泡轉(zhuǎn)運(yùn)，可以加速GTP取代結(jié)合GDP而參與ARFs的激活，并參與高爾基體小泡的運(yùn)輸；LOC784462基因有1個(gè)可變剪接數(shù)，是編碼Aurora激酶B的假定基因；IRX1基因有1個(gè)可變剪接數(shù)，編碼易洛魁家族同源家族蛋白，可與DNA結(jié)合。

2.2.5 產(chǎn)奶量候選基因

產(chǎn)奶量候選基因包括PCDH7基因和GFPT2基因。PCDH7基因有8個(gè)可變剪接數(shù)，編碼非成簇原鈣黏蛋白家族蛋白，與鈣離子結(jié)合，在細(xì)胞識別和粘附中起作用；GFPT2基因有8個(gè)可變剪接數(shù)，是果糖-6-磷酸酶酰胺轉(zhuǎn)移酶和氨基己糖生物合成途徑的限速酶，參與能量代謝。

2.2.6 體細(xì)胞數(shù)候選基因

體細(xì)胞數(shù)候選基因包括DYRK2、SULT4A1、SCARA5、LOC104969301、TBC1D32、AVEN、NACC2、LOC782305、LOC782610、PRICKLE2、FHIT、ZBTB7C、GSTP1這13個(gè)候選基因。

DYRK2基因有6個(gè)可變剪接數(shù)，是雙特異性酪氨酸磷酸化調(diào)節(jié)激酶2，參與細(xì)胞生長和發(fā)育過程，可調(diào)節(jié)細(xì)胞凋亡以應(yīng)對DNA損傷，其磷酸化可通過泛素-蛋白酶體系誘發(fā)蛋白降解，可與鎂離子結(jié)合；SULT4A1基因有4個(gè)可變剪接數(shù)，編碼腦特異性硫轉(zhuǎn)移酶，參與神經(jīng)遞質(zhì)代謝；SCARA5基因有4個(gè)可變剪接數(shù)，是A類清道夫受體5型，可借助特異性配體和受體分子識別和清楚外來物質(zhì)和體內(nèi)廢物，參與免疫防御；LOC104969301基因有1個(gè)可變剪接數(shù)，編碼長鏈非編碼RNA；TBC1D32基因有9個(gè)可變剪接數(shù)，編碼TBC結(jié)構(gòu)域包含蛋白，與CDK20一起，通過協(xié)調(diào)纖毛膜和軸突的組裝來控制初級纖毛的結(jié)構(gòu)，參與SHH信號通路，與細(xì)胞周期相關(guān)激酶互相作用；AVEN基因有3個(gè)可變剪接數(shù)，是凋亡和Caspase激活抑制因子，能通過結(jié)合凋亡蛋白激活因子-1干擾器自交連能力抑制細(xì)胞凋亡蛋白酶的水解作用，進(jìn)而抗細(xì)胞凋亡；NACC2基因有3個(gè)可變剪接數(shù)，能夠抑制細(xì)胞增殖，對DNA損傷致使細(xì)胞凋亡敏感，是P53通路上重要的轉(zhuǎn)錄調(diào)節(jié)因子；LOC782305基因有1個(gè)可變剪接數(shù)，編碼小泛素相關(guān)修飾物類蛋白；LOC782610基因有1個(gè)可變剪接數(shù)，編碼煙酰胺N-甲基轉(zhuǎn)移酶類蛋白；PRICKLE2基因有9個(gè)可變剪接數(shù)，編碼PET家族蛋白，有1個(gè)PET和3個(gè)半胱氨酸-組氨酸、鋅配位的LIM結(jié)構(gòu)域、N-糖基化位點(diǎn)、環(huán)磷酸腺苷依賴的蛋白激酶磷酸化位點(diǎn)、3個(gè)核定位信號和1個(gè)C端異戊二烯化化標(biāo)記，可以促進(jìn)非典型Wnt信號傳遞，又可促進(jìn)無序降解而抑制典型Wnt通路，可與金屬離子結(jié)合；FHIT基因有8個(gè)可變剪接數(shù)，是脆性組氨酸三聯(lián)體，編碼二腺苷P1-P3雙三磷酸水解酶，含有FRA3B脆性位點(diǎn)，脆性區(qū)域有較多Alu序列，參與嘌呤代謝，參與細(xì)胞增殖和凋亡；ZBTB7C基因有28個(gè)可變剪接數(shù)，是BTB-POZ家族轉(zhuǎn)錄因子，參與糖異生的調(diào)節(jié)因子，可以被PIAS1糖基化，通過泛素介導(dǎo)的蛋白酶體途徑促進(jìn)Kr-Pok的降解；GSTP1基因有9個(gè)可變剪接數(shù)，編碼谷胱甘肽S-轉(zhuǎn)移酶P1，通過催化疏水和親電化合物與還原型谷胱甘肽結(jié)合，發(fā)揮解毒作用。

2.2.7 初產(chǎn)日齡候選基因

初產(chǎn)日齡候選基因?yàn)镋PRS基因，有6個(gè)可變剪接數(shù)，編碼谷氨酰-脯氨酰-tRNA合成酶，催化谷氨酸和脯氨酸分別結(jié)合不同tRNA，通過選擇性沉默阻遏血漿銅藍(lán)蛋白的表達(dá)參與炎癥反應(yīng)調(diào)控過程，通過SLC27A1促進(jìn)脂肪細(xì)胞攝取長鏈脂肪酸，發(fā)揮在mTORC 1信號通路中的作用。

2.2.8 初配日齡候選基因

初配日齡候選基因包括KIF1A、RSPO2、EPRS等3個(gè)基因。KIF1A基因有31個(gè)可變剪接數(shù)，編碼驅(qū)動蛋白家族3，參與將前突觸囊泡從胞體沿著微管軌道向軸突末端的運(yùn)輸，也參與高爾基體到ER的逆行轉(zhuǎn)運(yùn)以及巨核細(xì)胞發(fā)育和血小板生成；RSPO2基因和EPRS基因分別有8個(gè)和6個(gè)可變剪接數(shù)，RSPO2基因編碼富含半胱氨酸的分泌蛋白，由Pspo1、Pspo2、Pspo3、Pspo4組成，具有50%左右的序列同源性和相似的結(jié)構(gòu)域，是Wnt信號通路的胞外激動劑，能與LGR4/5/6受體相互作用，調(diào)控Wnt/β-catenin信號通路和非典型Wnt信號通路。

2.2.9 妊娠長度候選基因

妊娠長度候選基因包括LOC786994、LOC785220、SPRY1、CACNG6、STAC這5個(gè)基因。LOC786994基因有1個(gè)可變剪接數(shù)，編碼核仁RNA解旋酶2假基因；LOC785220基因有1個(gè)可變剪接數(shù)，編碼小著絲粒相關(guān)蛋白假基因；SPRY1基因有9個(gè)可變剪接數(shù)，編碼軟脂?；椎鞍譙prouty1，負(fù)反饋調(diào)節(jié)受體酪氨酸激酶(RTK)信號通路，調(diào)控細(xì)胞增殖和凋亡，受生長因子的調(diào)節(jié)；CACNG6基因有3個(gè)可變剪接數(shù)，編碼鈣離子電壓門控通道輔助亞基γ6蛋白，參與電壓門控離子通道的調(diào)節(jié)；STAC基因有8個(gè)可變剪接數(shù)，編碼SH3富含半胱氨酸結(jié)構(gòu)域蛋白，可與金屬離子結(jié)合，參與信號轉(zhuǎn)導(dǎo)和電壓門控離子通道的調(diào)節(jié)。

2.2.10 產(chǎn)犢間隔候選基因

產(chǎn)犢間隔候選基因包括SHISA9、USP9X、LOC104970790、HCLS1、ANKFN1、RAI14等6個(gè)。SHISA9基因有5個(gè)可變剪接數(shù)，編碼Shisa家族蛋白，通過阻滯Wnt和Fgf通路的分子成熟或?qū)⑺麄兊氖荏w轉(zhuǎn)運(yùn)到細(xì)胞表面；USP9X基因有8個(gè)可變剪接數(shù)，編碼X染色體連鎖的泛素特異性蛋白酶，可通過囊泡網(wǎng)絡(luò)定向蛋白質(zhì)運(yùn)輸參與信號轉(zhuǎn)導(dǎo)，參與細(xì)胞增殖和細(xì)胞凋亡；LOC104970790基因有1個(gè)可變剪接數(shù)，編碼長鏈非編碼RNA；HCLS1基因有12個(gè)可變剪接數(shù)，是造血細(xì)胞特異性酪氨酸激酶的作用底物，參與胞外刺激因子介導(dǎo)的免疫應(yīng)答以及生長因子信號傳導(dǎo)通路；ANKFN1基因有9個(gè)可變剪接數(shù)，編碼錨蛋白重復(fù)序列和纖維連接蛋白Ⅲ型結(jié)構(gòu)域蛋白1；RAI14基因有28個(gè)可變剪接數(shù)，編碼一種肌動蛋白結(jié)合蛋白，參與絲狀肌動蛋白的重構(gòu)過程，參與建立精子極性和正常精子細(xì)胞的黏附作用，促進(jìn)支持細(xì)胞在血-睪丸屏障的緊密連接的完整性。

2.3 候選基因的GO分類

分別對上述新疆褐牛產(chǎn)奶和繁殖性狀的47個(gè)候選基因進(jìn)行GO分類，主要通過分子功能(molecular function，MF)、生物學(xué)過程(biological process BP)和細(xì)胞組分(cellular component)3個(gè)方面進(jìn)行生物信息學(xué)分析。在奶牛產(chǎn)奶性狀33個(gè)候選基因中，共計(jì)分布于18個(gè)GO terms中(圖1)；繁殖性狀14個(gè)候選基因中，共計(jì)有12個(gè)GO terms(圖2)。在產(chǎn)奶性狀和繁殖性狀中，富集基因最多的GO條目都是細(xì)胞進(jìn)程，產(chǎn)奶性狀的候選基因參與其中的有GFPT2、GPR33、CNIH3、CLSTN2、TBC1D32、CDH2、GABRG2及他們的可變剪接體，繁殖性狀的候選基因參與其中的有STAC、HCLS1、SPRY1和KIF1A，其次是產(chǎn)奶性狀中的細(xì)胞有8個(gè)基因分別是GABRG2、GPR33、CDH2、LOC538060、CLSTN2和TBC1D32等，繁殖性狀中的催化活性有5個(gè)基因分別是EPRS、SPRY1、KIF1A、USP9X和其可變剪接。

圖1 新疆褐牛產(chǎn)奶性狀候選基因GO分類Fig.1 GO classification of candidate genes for milk production traits in Xinjiang Brown cattle

2.4 候選基因的KEGG pathway富集分析

進(jìn)一步對新疆褐牛候選基因進(jìn)行KEGG通路富集分析，結(jié)果顯示，產(chǎn)奶性狀候選基因富集于4 pathways，繁殖性狀富集于3個(gè)pathways(圖3)，其中Wnt信號通路和鈣黏素信號通路基因數(shù)最多，且CDH2和PCDH7基因都有參與其中，表明此次GWAS篩選出的候選基因中大多數(shù)參與其中。

3 討論

本研究基于前期全基因組關(guān)聯(lián)分析得到的與產(chǎn)奶和繁殖性狀顯著關(guān)聯(lián)的SNPs位點(diǎn)，通過比對各多肽位點(diǎn)所在染色體的物理位置可知有29個(gè)位點(diǎn)位于基因間隔區(qū)，有28個(gè)位點(diǎn)位于基因的內(nèi)含子區(qū)域。雖然商用基因芯片的密度和數(shù)量很高，但其采用的SNP位點(diǎn)多位于基因非編碼區(qū)和基因間隔區(qū)，大多不是和性狀存在直接關(guān)聯(lián)的致因變異，而是通過與致變異相連鎖而顯示出關(guān)聯(lián)性[16]。因此，GWAS分析也只是鑒定出和目標(biāo)性狀相關(guān)染色體片段，但不同群體中可能具有不同的等位基因頻率和不同的連鎖不平衡區(qū)域，導(dǎo)致非功能性的標(biāo)簽SNP在不同群體中的應(yīng)用具有不確定性[17]。但是顯著的SNP位點(diǎn)很可能與真正影響產(chǎn)奶和繁殖性狀的SNP緊密連鎖，雖然得到的大部分候選基因從后續(xù)功能注釋的結(jié)果中并不能明確參與到泌乳和繁殖過程，只能判斷他們所調(diào)控的轉(zhuǎn)錄和翻譯等生物過程以及生成的蛋白可能與產(chǎn)奶和繁殖性狀相關(guān)聯(lián)，與畜禽上的其他重要性狀GWAS研究的結(jié)果相類似[18]。近年來全基因組關(guān)聯(lián)分析在人類復(fù)雜疾病、豬和奶牛生產(chǎn)性狀上的大量應(yīng)用證實(shí)了在全基因組范圍內(nèi)發(fā)掘目標(biāo)性狀顯著相關(guān)新位點(diǎn)上的高效性[19-21]，實(shí)現(xiàn)對重要性狀的遺傳標(biāo)記和基因的鑒別，實(shí)現(xiàn)個(gè)體基因組選擇。

圖2 新疆褐牛繁殖性狀候選基因GO分類Fig.2 GO classification of candidate genes for reproductive traits in Xinjiang Brown cattle

圖3 新疆褐牛產(chǎn)奶性狀、繁殖性狀候選基因KEGG通路分析Fig.3 KEGG pathway analysis of candidate genes for milk and reproductive traits in Xinjiang Brown cattle

通過功能基因組信息學(xué)分析GWAS結(jié)果得到的47個(gè)候選基因共分布于30個(gè)GO terms，富集在7個(gè)pathways上，GO分析得到富集基因數(shù)最多的是細(xì)胞進(jìn)程、細(xì)胞和催化活性，表明大部分候選基因參與了新疆褐牛的生產(chǎn)發(fā)育和代謝。KEGG代謝通路富集發(fā)現(xiàn)候選基因參與最多的是Wnt通路，調(diào)控細(xì)胞增殖、分化和遷移過程[22]，對脂肪細(xì)胞的影響表現(xiàn)為對間充質(zhì)干細(xì)胞(MSC)轉(zhuǎn)化為脂肪前體細(xì)胞和脂肪細(xì)胞終末分化為成熟脂肪細(xì)胞的兩個(gè)階段有調(diào)控功能，抑制脂肪細(xì)胞的分化。細(xì)胞黏附分子是調(diào)節(jié)細(xì)胞之間相互作用的膜表面糖蛋白，鈣黏蛋白是最早發(fā)現(xiàn)得介導(dǎo)細(xì)胞間相互聚集的黏附分子，神經(jīng)鈣粘素是其中的一員，在CD抗原中被命名為CD325，編碼基因?yàn)镃DH2，由于該基因同時(shí)參與了Wnt通路和鈣黏素通路，并且GWAS得到SNP位點(diǎn)位于該基因內(nèi)，但是該基因與牛產(chǎn)奶性狀的相關(guān)研究很少，故將其篩選為乳脂率性狀的候選基因。同時(shí)，GABRG2基因主要參與γ-氨基丁酸門控氯離子通道的活動，并參與GABAA受體活性，其所在的QTL區(qū)域與乳蛋白量有關(guān)。EPRS基因可以通過SLC27A1促進(jìn)脂肪細(xì)胞攝取長鏈脂肪酸，發(fā)揮在mTORC1信號通路中的作用。mTORC1通路處于生長調(diào)節(jié)的中心環(huán)節(jié)，可以調(diào)控蛋白質(zhì)和核糖體的生物合成、營養(yǎng)物質(zhì)的吸收和細(xì)胞自噬等過程，初步認(rèn)為該基因與初配日齡有關(guān)。因此，通過生物學(xué)功能、生理生化分析初步確定產(chǎn)奶性狀中的CDH2、GABRG2基因和繁殖性狀中的EPRS基因?yàn)楹蜻x基因，為后續(xù)試驗(yàn)驗(yàn)證分子標(biāo)記對特定性狀的效應(yīng)提供前期理論基礎(chǔ)。

由于可變剪接的機(jī)制復(fù)雜多樣，所以在進(jìn)行功能注釋的同時(shí)也對基因的可變剪接數(shù)進(jìn)行統(tǒng)計(jì)。首先調(diào)控mRNA剪接位點(diǎn)的有順式和反式作用元件，其活性受細(xì)胞內(nèi)剪接因子之間的互作和競爭調(diào)控來進(jìn)一步影響機(jī)體剪切位點(diǎn)的選擇，其次不同啟動子的選擇和結(jié)構(gòu)，以及募集剪接因子的能力和調(diào)控轉(zhuǎn)錄速度的差異都會影響可變剪接。此外，RNA的二級結(jié)構(gòu)、細(xì)胞因子、激素等也能影響剪接位點(diǎn)的識別。目前研究可變剪接的方法有實(shí)時(shí)熒光定量RCR法、表達(dá)標(biāo)簽序列分析法以及生物芯片法。mRNA的可變剪接增加了蛋白質(zhì)組的多樣性，反映出遺傳信息在基因組水平外的變異和重組，增加蛋白質(zhì)表達(dá)的復(fù)雜性，使機(jī)體利用有限的基因滿足多樣性功能蛋白的需求。根據(jù)GWAS發(fā)掘功能位點(diǎn)所得到的候選基因需要進(jìn)一步通過基因表達(dá)譜分析、功能試驗(yàn)等方法進(jìn)一步確定其具體的基因功能。

4 小結(jié)

本研究基于課題組前期GWAS研究結(jié)果，對57個(gè)顯著SNPs位點(diǎn)進(jìn)行功能基因組生物信息學(xué)分析，經(jīng)過篩選后得到47個(gè)候選基因，隨后進(jìn)行了GO和KEGG分析和功能注釋。根據(jù)GWAS得到顯著的位點(diǎn)與候選基因的位置和候選基因本身的功能信息，對產(chǎn)奶性狀中CDH2、GABRG2基因和繁殖性狀中EPRS基因可進(jìn)一步進(jìn)行系列的基因功能研究。