基于全基因組重測(cè)序的山羊產(chǎn)羔數(shù)性狀關(guān)鍵調(diào)控基因的篩選

李恒，字向東，王會(huì)，熊燕，呂明杰，劉宇，蔣旭東

基于全基因組重測(cè)序的山羊產(chǎn)羔數(shù)性狀關(guān)鍵調(diào)控基因的篩選

李恒1，字向東1，王會(huì)2，熊燕1，呂明杰1，劉宇1，蔣旭東1

1西南民族大學(xué)動(dòng)物科學(xué)國(guó)家民委重點(diǎn)實(shí)驗(yàn)室，成都 610041；2西南民族大學(xué)青藏高原動(dòng)物遺傳資源保護(hù)與利用教育部重點(diǎn)實(shí)驗(yàn)室，成都 610041

【目的】對(duì)產(chǎn)羔數(shù)不同的山羊進(jìn)行全基因組重測(cè)序分析，挖掘參與調(diào)控川中黑山羊產(chǎn)羔數(shù)性狀關(guān)鍵調(diào)控基因，為解析山羊產(chǎn)羔數(shù)性狀遺傳機(jī)制及分子遺傳改良提供理論依據(jù)?！痉椒ā窟x擇6只產(chǎn)4—6羔的川中黑山羊?yàn)楦叻苯M（high fecundity, HF）和6只產(chǎn)單羔的川中黑山羊?yàn)榈头苯M（low fecundity, LF），采集頸靜脈血液樣本提取基因組DNA，構(gòu)建350 bp雙末端測(cè)序文庫(kù)，利用Illumina HiSeq PE150平臺(tái)對(duì)12個(gè)文庫(kù)進(jìn)行全基因組重測(cè)序。測(cè)序產(chǎn)出的凈數(shù)據(jù)經(jīng)BWA軟件比對(duì)至山羊參考基因組ARS1，所獲得的高質(zhì)量SNPs通過(guò)兩種全基因組掃描分析方法（、）的綜合分析確定候選區(qū)域，候選區(qū)域的注釋基因分別利用g:Profiler和KOBAS在線數(shù)據(jù)庫(kù)進(jìn)行GO分析與KEGG通路分析，以篩選調(diào)節(jié)川中黑山羊產(chǎn)羔數(shù)性狀候選基因。為了進(jìn)一步鑒定調(diào)節(jié)山羊產(chǎn)羔數(shù)目的關(guān)鍵遺傳標(biāo)記，根據(jù)基因組重測(cè)序變異分析，對(duì)繁殖候選基因的同義與非同義SNPs進(jìn)行定位篩選，后續(xù)將12個(gè)山羊樣本的擴(kuò)增產(chǎn)物進(jìn)行Sanger測(cè)序以驗(yàn)證重測(cè)序結(jié)果?！窘Y(jié)果】12只山羊共獲得431.50 Gb 凈數(shù)據(jù)，經(jīng)變異檢測(cè)與注釋發(fā)現(xiàn)，HF組山羊共發(fā)現(xiàn)7 771 417個(gè)單核苷酸多態(tài)性(single nucleotide polymorphism, SNPs)，LF組山羊檢測(cè)到8 935 907個(gè)SNPs，且LF組各類(lèi)SNPs 均多于HF組。設(shè)置同時(shí)達(dá)到top 5%最大值和top 5%最小值的窗口為候選區(qū)域，在低雜合性、高遺傳分化的區(qū)域共注釋130個(gè)強(qiáng)選擇信號(hào)，其中HF組、LF組以及共享窗口的注釋基因分別為84、59和13個(gè)，經(jīng)GO富集與KEGG通路分析發(fā)現(xiàn)，19個(gè)候選基因參與川中黑山羊的繁殖、繁殖過(guò)程和胚胎發(fā)育等調(diào)控，包括11個(gè)HF組特異性候選基因（和），5個(gè)LF組特異性候選基因（、、、和）和3個(gè)HF組與LF組共享窗口基因（、和）。同時(shí)，大多數(shù)GO分析，如G蛋白偶聯(lián)受體活性、激素反應(yīng)和神經(jīng)肽信號(hào)通路等，都包含這19個(gè)候選基因。此外，14個(gè)HF候選基因有9個(gè)顯著富集在代謝途徑、神經(jīng)活性配體-受體相互作用、糖胺聚糖-硫酸乙酰肝素/肝素的生物合成、鈣離子信號(hào)通路、cAMP信號(hào)通路和葉酸生物合成等KEGG通路中（＜0.05）。19個(gè)繁殖候選基因中共有2個(gè)同義突變（，10 A4662G）與2個(gè)非同義突變（G529A，），且僅定位于HF候選基因中。Sanger測(cè)序發(fā)現(xiàn)，、和突變位點(diǎn)均可檢測(cè)到多態(tài)性，與基因組重測(cè)序結(jié)果一致，其中G529A多態(tài)性導(dǎo)致第177位丙氨酸突變?yōu)樘K氨酸，多態(tài)性導(dǎo)致翻譯提前終止?！窘Y(jié)論】本研究共發(fā)現(xiàn)11個(gè)HF組特異性候選基因，推測(cè)是川中黑山羊多羔性狀的關(guān)鍵調(diào)控基因，外顯子G529A與外顯子A281T突變可能是調(diào)控山羊多羔性狀的關(guān)鍵遺傳標(biāo)記，在改良山羊繁殖性能方面具有較大的應(yīng)用價(jià)值。

川中黑山羊；基因組重測(cè)序；多羔性狀；候選基因

0 引言

【研究意義】產(chǎn)羔性狀是山羊的重要經(jīng)濟(jì)性狀，提高窩產(chǎn)羔數(shù)不僅可以提高出欄率，提高經(jīng)濟(jì)效益，而且還可以提高選擇強(qiáng)度，加快山羊育種進(jìn)程。產(chǎn)羔數(shù)是一個(gè)低遺傳力（0.08—0.18）的限性性狀，難以用常規(guī)育種方法進(jìn)行快速改良，但適合用標(biāo)記輔助選擇（marker-assisted selection, MAS）等分子育種技術(shù)來(lái)改良[1]。實(shí)現(xiàn)MAS的先決條件是找到與數(shù)量性狀基因座相連鎖的分子遺傳標(biāo)記，而要找到輔助選擇的分子標(biāo)記，則需要解析山羊高繁殖力的分子遺傳基礎(chǔ)和形成原因?！厩叭搜芯窟M(jìn)展】在過(guò)去20多年里，在一些綿羊品種的突變體中已成功分離鑒定出控制卵泡發(fā)育和排卵數(shù)的多胎基因，主要包括骨形態(tài)發(fā)生蛋白受體-1B（BMP receptor-1B，）基因的1個(gè)突變（FecB）[2]、骨形態(tài)發(fā)生蛋白-15（bone morphogenetic protein-15，）基因的6個(gè)突變（FecX、FecX、FecX、FecX、FecX和FecX）[3-4]和生長(zhǎng)分化因子-9（growth and differentiation factor-9，）基因的2個(gè)突變（FecG和）[5-6]，這些突變與綿羊的多羔性狀密切相關(guān)。但是，迄今為止，在已開(kāi)展過(guò)檢測(cè)的印度[7]、伊朗[8]、泰國(guó)[9]和我國(guó)[10-11]的30多個(gè)山羊品種（類(lèi)群）中均未檢測(cè)到這些突變，排除了這些突變作為山羊品種多羔性狀主基因的可能性，山羊多羔性狀的遺傳機(jī)制有待研究。山羊全基因組測(cè)序技術(shù)的發(fā)展與完善，為發(fā)掘參與調(diào)控川中黑山羊產(chǎn)羔數(shù)性狀關(guān)鍵調(diào)控基因提供了全新的視角[12]。對(duì)嶗山奶山羊高、低繁殖力的兩個(gè)極端種群進(jìn)行基因組重測(cè)序，發(fā)現(xiàn)細(xì)胞周期蛋白B2（cyclin B2,2）雄激素受體（androgen receptor,）腺苷酸環(huán)化酶1（adenylate cyclase 1,）SMAD家庭成員2（SMAD family member 2,）等基因在高繁殖力群體中被特異性選擇[13]。絲氨酸/蘇氨酸激酶3（serine/threonine kinase 3,）蛋白磷酸酶3催化亞基α （protein phosphatase 3 catalytic subunit alpha,）等96個(gè)候選基因與大足黑山羊產(chǎn)羔數(shù)性狀顯著相關(guān)[14]，其中也與阿爾巴斯絨山羊產(chǎn)羔數(shù)性狀相關(guān)[15]。進(jìn)一步分析67個(gè)單核苷酸多態(tài)性（SNP）與大足黑山羊產(chǎn)羔數(shù)的關(guān)聯(lián)，發(fā)現(xiàn)半胱氨酰轉(zhuǎn)運(yùn)RNA合成酶2（cysteinyl-tRNA synthetase 2，）Ⅰ型血小板結(jié)合蛋白基序的解聚蛋白樣金屬蛋白酶（ADAM metallopeptidase with thrombospondin type 1 motif 14，）和甲基轉(zhuǎn)移酶 25（methyltransferase like 25，）基因編碼區(qū)SNP與頭胎產(chǎn)羔數(shù)目顯著相關(guān)[16]。對(duì)頭胎產(chǎn)單羔、雙羔、三羔的三組濟(jì)寧青山羊進(jìn)行全基因組掃描，結(jié)果在雙羔和三羔群體中發(fā)現(xiàn)酪氨酸激酶受體（KIT proto-oncogene, receptor tyrosine kinase，）、P21-活化激酶1P21-activated kinases 1，腺苷活化蛋白激酶α1亞基（protein kinase AMP-activated catalytic subunit alpha 1，）等多個(gè)產(chǎn)羔數(shù)調(diào)節(jié)基因[17]。WANG等[18]在相同飼養(yǎng)管理?xiàng)l件下，對(duì)在5個(gè)連續(xù)分娩記錄中表現(xiàn)出穩(wěn)定的產(chǎn)羔數(shù)差異的12只山羊進(jìn)行基因組重測(cè)序，結(jié)果發(fā)現(xiàn)細(xì)胞分裂周期蛋白25同源蛋白C（cell division cycle 25C，）、核酸內(nèi)切酶G（endonuclease G，）和Nanos同源基因3（nanos C2HC-type zinc finger 3，）的變異與產(chǎn)羔數(shù)性狀顯著相關(guān)，螺旋卷曲過(guò)程對(duì)生殖能力有潛在的調(diào)控作用。這些候選基因的發(fā)現(xiàn)拓展了人們對(duì)繁殖力遺傳基礎(chǔ)的認(rèn)知，為山羊產(chǎn)羔數(shù)性狀的遺傳機(jī)制的解析提供重要線索?！颈狙芯壳腥朦c(diǎn)】川中黑山羊是我國(guó)優(yōu)良地方山羊品種，具有生長(zhǎng)速度快，產(chǎn)肉性能高，產(chǎn)羔率高等優(yōu)良特性[19-20]。我們?cè)诖ㄖ泻谏窖蚝偷那捌谘芯恐形礄z測(cè)到與綿羊多羔性狀相關(guān)突變的存在，但高繁川中黑山羊與單胎藏山羊之間檢測(cè)到5個(gè)新的堿基突變，其中4個(gè)導(dǎo)致氨基酸改變，有2個(gè)堿基突變，其中1個(gè)導(dǎo)致氨基酸突變，兩個(gè)品種的15核苷酸序列則完全一致[21]。說(shuō)明山羊和中氨基酸改變的突變與綿羊中不同，這些突變是否與山羊產(chǎn)羔數(shù)相關(guān)尚有待進(jìn)一步研究。從卵巢轉(zhuǎn)錄組和miRNA水平也不能完全解析川中黑山羊多羔性狀的分子遺傳機(jī)制[22-23]?！緮M解決的關(guān)鍵問(wèn)題】本研究選擇高、低繁殖力的川中黑山羊?yàn)檠芯繉?duì)象，利用基因組重測(cè)序技術(shù)掃描其低雜合性、高遺傳分化區(qū)域，挖掘高、低繁殖力山羊關(guān)鍵基因遺傳變異，并對(duì)候選靶基因進(jìn)行GO富集與KEGG通路分析，以期篩選調(diào)節(jié)川中黑山羊產(chǎn)羔數(shù)的關(guān)鍵候選基因，為山羊產(chǎn)羔數(shù)性狀遺傳機(jī)制的闡釋及后續(xù)山羊分子育種工作提供理論基礎(chǔ)。

1 材料與方法

1.1 樣品的采集

本研究根據(jù)第一、二胎產(chǎn)羔記錄，于2019年9月選擇飼養(yǎng)于四川省樂(lè)至縣天龍育種場(chǎng)，胎產(chǎn)羔數(shù)差異顯著的高繁殖力（high fecundity，HF）和低繁殖力（low fecundity，LF）三歲齡川中黑山羊各6只，其中HF組山羊連續(xù)兩胎胎產(chǎn)羔數(shù)4—6只，LF組連續(xù)兩胎胎產(chǎn)羔數(shù)1只（表1）。每只山羊采集其頸靜脈血液2 mL于EDTA-K2抗凝管中，試劑盒提取基因組DNA后用于基因組重測(cè)序研究。

1.2 山羊血液基因組DNA的提取及質(zhì)量檢測(cè)

使用1%瓊脂糖凝膠電泳與紫外分光光度計(jì)分別鑒定山羊基因組DNA的完整性、濃度和純度，12管質(zhì)檢合格的HF組和LF組DNA樣品各抽取1.5 μg分別建庫(kù)。利用Agilent 2100檢測(cè)文庫(kù)中插入片段的長(zhǎng)度，片段符合預(yù)期后，送至Illumina HiSeq PE150平臺(tái)進(jìn)行雙端測(cè)序，整個(gè)流程委托北京諾禾致源生物信息科技有限公司完成。

表1 12只川中黑山羊的胎產(chǎn)羔數(shù)信息

編號(hào)HF1—6為高繁殖力山羊，編號(hào)LF1—6為低繁殖力山羊

The symbols from HF1 to HF6 were high fecundity goats, and the symbols from LF1 to LF6 were low fecundity goats

1.3 測(cè)序數(shù)據(jù)的變異檢測(cè)與注釋

原始數(shù)據(jù)經(jīng)初步質(zhì)控處理，將接頭序列、未測(cè)出堿基超過(guò)10%和低質(zhì)量堿基數(shù)超過(guò)50%的paired reads去除，獲得凈數(shù)據(jù)。利用BWA[24]軟件將clean reads比對(duì)至山羊參考基因組ARS1[25]，比對(duì)結(jié)果經(jīng)SAMTOOLS軟件去除重復(fù)[26]。利用GATKv3.2-2進(jìn)行SNP鑒定與分型[27]，過(guò)濾后獲得的高質(zhì)量SNPs用ANNOVAR[28]軟件工具進(jìn)行功能注釋。

1.4 群體選擇分析

利用滑動(dòng)窗口法（sliding windows）對(duì)HF、LF組川中黑山羊進(jìn)行全基因組掃描，按照窗口內(nèi)SNPs個(gè)數(shù)不超過(guò)20，確定選擇消除窗口大?。?00 kb）。計(jì)算每個(gè)窗口的固定系數(shù)（Fst）值并轉(zhuǎn)換為值，選擇top 5% 最大t值窗口為選擇區(qū)域。通過(guò)計(jì)算窗口內(nèi)SNPs位點(diǎn)的雜合性（Hp），進(jìn)而對(duì)Selective sweep進(jìn)行評(píng)估。分別計(jì)算窗口的Hp值并轉(zhuǎn)化為值，選擇top 5%最小值窗口作為選擇區(qū)域。在計(jì)算得到每個(gè)窗口Hp和Fst值的基礎(chǔ)上，選擇最小值和最大值均達(dá)到Top 5%的區(qū)域?yàn)楹蜻x區(qū)域，并對(duì)該區(qū)域進(jìn)行基因注釋。

1.5 候選基因的功能注釋

基于Bio-Mart在線數(shù)據(jù)庫(kù)獲得山羊候選區(qū)域內(nèi)的基因注釋信息（基因類(lèi)型和ENSEMBL號(hào)）后，利用g:Profiler網(wǎng)站將山羊基因符號(hào)ID轉(zhuǎn)換為小鼠同源基因。利用g:Profiler[29]網(wǎng)站和KOBAS[30]在線數(shù)據(jù)庫(kù)分別進(jìn)行Gene ontology（GO）富集分析和Kyoto Encyclopedia of Genes and Genomes（KEGG）通路分析。

1.6 SNP位點(diǎn)的驗(yàn)證

為進(jìn)一步鑒定調(diào)節(jié)山羊產(chǎn)羔數(shù)目的關(guān)鍵遺傳標(biāo)記，據(jù)重測(cè)序變異分析報(bào)告對(duì)19個(gè)繁殖候選基因的同義與非同義SNPs進(jìn)行定位篩選。候選位點(diǎn)經(jīng)Primer Premier 5.0 軟件設(shè)計(jì)引物（表2），12個(gè)山羊樣本的擴(kuò)增產(chǎn)物送至成都擎科生物有限公司進(jìn)行Sanger測(cè)序。

2 結(jié)果

2.1 高、低繁殖力組山羊基因組測(cè)序數(shù)據(jù)質(zhì)控

分別完成12只山羊血液樣本基因組重測(cè)序工作，共獲得431.50 Gb 凈數(shù)據(jù)，錯(cuò)誤分布率均為0.03%，測(cè)序質(zhì)量較高（Q20≥96.9%，Q30≥91.69%），所有樣本比對(duì)率在99.61%—99.73%之間，平均測(cè)序深度10.70×，1×覆蓋度在94.91%—95.32%之間，4×覆蓋度在88.70%—93.13%之間。所有原始序列數(shù)據(jù)已上傳至國(guó)家基因組數(shù)據(jù)中心基因組序列檔案庫(kù)（https://bigd.big.ac.cn/gsa.），注冊(cè)號(hào)為CRA003846。

表2 SNP位點(diǎn)驗(yàn)證引物

2.2 變異檢測(cè)與注釋

基于嚴(yán)格的閾值，兩組山羊共檢測(cè)到16 707 324個(gè)高質(zhì)量SNPs位點(diǎn)，其中LF組8 935 907個(gè)SNPs位點(diǎn)，HF組7 771 417個(gè)SNPs位點(diǎn)，且LF組的各類(lèi)SNPs 均多于HF組（表3）。在這些SNPs中，HF組和LF組山羊外顯子SNPs平均各有53 425（0.68%）和60 588（0.69%）個(gè)，基因上游1 kb區(qū)域分別有44 949個(gè)和50 798個(gè)，基因下游1 kb區(qū)域分別分布44 201個(gè)和50 349個(gè)，兩組基因間區(qū)SNPs占比最多，分別占總SNPs的70.01%和70.09%；其次是分布在內(nèi)含子區(qū)域SNPs，HF和LF組分別含有2 145 836個(gè)和2 473 891個(gè)；位于剪切位點(diǎn)的SNPs數(shù)目最少。在編碼區(qū)SNPs 中，同義突變和錯(cuò)義突變分別占外顯子總SNPs的58.66%和40.77%，且兩組相差不大。此外，HF組和LF組在終止密碼子處分別檢測(cè)到300和323個(gè)SNPs，在剪切位點(diǎn)處分別檢測(cè)到214個(gè)和237個(gè)SNPs，這些位點(diǎn)可能會(huì)影響轉(zhuǎn)錄剪接，從而影響蛋白質(zhì)產(chǎn)物及功能。

表3 川中黑山羊SNPs檢測(cè)及注釋結(jié)果統(tǒng)計(jì)

2.3 群體選擇分析

經(jīng)不同大小窗口調(diào)試發(fā)現(xiàn)，窗口大小為100 kb時(shí)，SNPs數(shù)目小于20個(gè)并逐漸趨于穩(wěn)定，因此，選擇100 kb作為選擇信號(hào)檢測(cè)最佳窗口長(zhǎng)度，使用50%的重疊區(qū)為滑動(dòng)尺，計(jì)算川中黑山羊HF組與LF組之間的Fst值，并標(biāo)準(zhǔn)化作圖，選取LF和HF之間的閾值為1.9，共篩選到579個(gè)窗口，注釋623個(gè)基因（圖1）。

分別計(jì)算川中黑山羊HF組與LF組染色體窗口內(nèi)SNP位點(diǎn)的值，Z轉(zhuǎn)換后并標(biāo)準(zhǔn)化作圖，選擇top 5%最小值為閾值，HF群體的閾值為-2.05，篩選到579個(gè)窗口（圖2），共注釋617個(gè)候選基因。 LF群體的閾值為-2.08，篩選到579個(gè)窗口（圖3），注釋到610個(gè)基因。

圖1 高繁組和低繁組川中黑山羊1-29號(hào)常染色體平均固定系數(shù)ZFst的分布

圖2 高繁組川中黑山羊1-29號(hào)常染色體平均雜合度ZHp的分布

為進(jìn)一步篩選川中黑山羊產(chǎn)羔性狀關(guān)鍵調(diào)控基因，本研究初步掃描了具有低雜合性、高遺傳分化的區(qū)域，即候選窗口值均達(dá)到top 5%最大值和top 5%最小值標(biāo)準(zhǔn)。其中HF組中鑒定65個(gè)候選窗口，共注釋84個(gè)候選基因（圖4-A），LF組山羊共鑒定了42個(gè)候選窗口，注釋了59個(gè)候選基因（圖4-B）。因兩組候選基因中存在、等13個(gè)交集基因，所以HF組和LF組共注釋了130個(gè)差異候選基因。

2.4 高、低繁殖力山羊候選基因基因GO富集分析

進(jìn)一步揭示與產(chǎn)羔數(shù)性狀的候選靶基因功能，本研究對(duì)兩群體差異表達(dá)基因（130個(gè)）進(jìn)行生物功能分析，-value經(jīng)Bonferroni校正后，設(shè)置corrected-value ≤0.05為閾值，選擇候選靶基因中顯著富集的GO terms。

圖3 低繁組川中黑山羊1-29號(hào)常染色體平均雜合度ZHp的分布

圖4 川中黑山羊高繁組（A）與低繁組（B）受選擇區(qū)域

GO結(jié)果顯示，川中黑山羊群體中篩選的130個(gè)強(qiáng)選擇性候選基因，共富集到193條顯著性GO功能條目。其中歸類(lèi)為生物學(xué)過(guò)程的GO條目最多，共113條，主要行使生長(zhǎng)發(fā)育和代謝功能，其次是細(xì)胞組成，共50條GO條目，分子功能本體條目最少，其富集的30條GO條目主要參與跨膜信號(hào)受體活性、神經(jīng)肽受體活性和蛋白多糖生物合成的過(guò)程等。

值得注意的是，19個(gè)基因與繁殖緊密相關(guān)，如繁殖、繁殖過(guò)程、胚胎發(fā)育和生殖過(guò)程調(diào)控（表4）。其中HF組特異性強(qiáng)選擇性基因?yàn)?1個(gè)，包含腺苷酸環(huán)化酶10（adenylate cyclase 10，）、多巴胺受體D1（dopamine receptor D1，）、硫酸肝素6-鄰磺酸轉(zhuǎn)移酶1（heparan sulfate 6-O-sulfotransferase 1，）、胰島素樣生長(zhǎng)因子結(jié)合蛋白7（insulin like growth factor binding protein 7，）、促黑激素同源框2（msh homeobox 2，）、頭蛋白（noggin，）、腎單位腎癆4（nephronophthisis 4 (juvenile) homolog，）、妊娠相關(guān)血漿蛋白A（pregnancy- associated plasma protein-A，、睪丸發(fā)育相關(guān)蛋白（testis development-Related Protein，）、木糖轉(zhuǎn)移酶1（xylosyltransferase 1，）和催乳素釋放激素受體（prolactin releasing hormone receptor，），HF和LF組共享窗口基因3個(gè)，如醛酮還原酶家族成員B3（aldo-keto reductase family 1，member B3，）、組蛋白脫乙酰酶4（histone deacetylase 4，）和mu阿片受體（opioid receptor mu 1，），以及LF組特異性候選基因5個(gè)，如膜聯(lián)蛋白5（annexin A5，）、內(nèi)皮素A型受體(endothelin receptor type A，)、FA核心復(fù)合體連接酶（FA complementation group L，）、胰島素樣生長(zhǎng)因子1（insulin-like growth factor1，）和速激肽前體1（tachykinin precursor 1，）。同時(shí)，大多數(shù)GO術(shù)語(yǔ)，如G蛋白偶聯(lián)受體活性、激素反應(yīng)和神經(jīng)肽信號(hào)通路等，都包含這19個(gè)候選基因。

表4 川中黑山羊繁殖相關(guān)GO條目

2.5 高、低繁殖力山羊候選基因KEGG分析

KEGG通路分析結(jié)果顯示，川中黑山羊差異候選基因集共富集112條代通路，選擇前20條顯著的代謝通路作為展示（圖5）。KEGG代謝通路結(jié)果表明，14個(gè)HF山羊候選基因中，9個(gè)基因顯著富集在代謝途徑（、）、糖胺聚糖-硫酸乙酰肝素/肝素的生物合成（、）、神經(jīng)活性配體-受體相互作用（）、鈣離子信號(hào)通路（）、cAMP信號(hào)通路（）和葉酸生物合成（）通路（corrected-value ≤0.05）中，這些通路可能在山羊產(chǎn)羔數(shù)性狀的調(diào)控中極具關(guān)鍵作用。

2.6 SNP位點(diǎn)的驗(yàn)證

為篩選多羔性狀關(guān)鍵遺傳標(biāo)記，本研究掃描LF和HF組外顯子SNPs分布。19個(gè)繁殖候選基因中，共有2個(gè)同義突變（，）與2個(gè)非同義突變（G529A，），且僅定位于HF候選基因中。其中G529A多態(tài)性導(dǎo)致第177位丙氨酸突變?yōu)樘K氨酸，多態(tài)性導(dǎo)致翻譯提前終止。重要的是，4個(gè)SNPs僅有3個(gè)與高通量測(cè)序結(jié)果呈現(xiàn)相同趨勢(shì)（，G529A，），的外顯子SNP測(cè)序結(jié)果與重測(cè)序變異報(bào)告不一致，未發(fā)現(xiàn)該突變位點(diǎn)（圖6）。

3 討論

山羊產(chǎn)羔數(shù)是復(fù)雜的表型特征，其不僅由排卵率、胚胎或胎兒成活率、子宮容受性等多個(gè)次級(jí)性狀調(diào)控，也受遺傳、季節(jié)、降雨、氣溫和營(yíng)養(yǎng)等因素影響[31-33]。GWAS雖是檢測(cè)復(fù)雜性狀的遺傳學(xué)基礎(chǔ)的最有力的工具[34]，但深入挖掘GWAS的統(tǒng)計(jì)能力需要巨大的樣本量和適當(dāng)?shù)膶?duì)照設(shè)計(jì)[35]。因此，大樣本、家系信息詳盡等要求限制了GWAS在山羊群體中的應(yīng)用?；蚪M重測(cè)序仍是目前分析山羊復(fù)雜性狀的強(qiáng)有力手段。該技術(shù)不僅可發(fā)掘出特定性狀的多個(gè)候選基因，其公布的變異數(shù)據(jù)也為其他性狀研究提供思考。例如最新發(fā)現(xiàn)的山羊產(chǎn)羔數(shù)候選基因、等[36-37]，遺傳標(biāo)記信息均是來(lái)自先前的山羊基因組重測(cè)序報(bào)道。

圖5 川中黑山羊選擇信號(hào)KEGG通路分析

本研究對(duì)高產(chǎn)山羊品種川中黑山羊不同產(chǎn)羔數(shù)群體進(jìn)行了全基因組測(cè)序，獲得了大量遺傳變異，為山羊高繁殖力遺傳機(jī)制解析提供了重要的遺傳資料。為闡明川中黑山羊產(chǎn)羔性狀的選擇特征，設(shè)置top 5%為閾值[13-14]，在川中黑山羊群體的低雜合性、高遺傳分化的區(qū)域共注釋130個(gè)候選基因。

通過(guò)GO富集進(jìn)一步發(fā)現(xiàn)了19個(gè)候選基因可能參與川中黑山羊產(chǎn)羔性狀的調(diào)節(jié)，參與繁殖、繁殖過(guò)程、胚胎發(fā)育等生殖過(guò)程調(diào)控。在這19個(gè)候選基因中，11個(gè)候選基因在HF組中特異性選擇，可能與川中黑山羊的多羔性狀有關(guān)。這些基因的功能主要與生殖細(xì)胞發(fā)育及胚胎發(fā)育等生殖過(guò)程相關(guān)。其中編碼的可溶性腺苷酸環(huán)化酶（SAC），對(duì)男性生殖至關(guān)重要[38]。SAC調(diào)控cAMP/PKA通路與AMPK的活性，其中cAMP/PKA通路在顆粒細(xì)胞FSH的功能調(diào)節(jié)和E2產(chǎn)生過(guò)程中發(fā)揮重要作用[39]，AMPK的活性影響卵母細(xì)胞成熟[40]。推測(cè)SAC通過(guò)調(diào)節(jié)生殖激素分泌水平進(jìn)而調(diào)控動(dòng)物繁殖過(guò)程，其表達(dá)的豐度制約了雌性個(gè)體排卵率。與與家禽產(chǎn)卵相關(guān)[41-42]，可能在山羊群體的排卵機(jī)制中發(fā)揮特殊作用。可能是卵巢早衰病因?qū)W的候選基因[43]，該基因突變也被認(rèn)為與特發(fā)性低促性腺激素減退癥有關(guān)[44]。牛繁殖候選基因下調(diào)促進(jìn)了細(xì)胞凋亡，阻礙了細(xì)胞增殖，增加了孕酮和雌二醇的產(chǎn)生[45]，其還與[46]功能類(lèi)似，對(duì)胚胎著床、妊娠的建立與維持至關(guān)重要[47]。抑制IGFBP7表達(dá)可顯著降低植入胚胎數(shù)與妊娠率[47]，這意味著除排卵數(shù)外，早期胚胎的丟失也是山羊窩產(chǎn)羔數(shù)目減少的誘因。目前并未發(fā)現(xiàn)功能突變與山羊早期胚胎丟失有關(guān)，其調(diào)控機(jī)制還需深入研究。在卵巢中表達(dá)并與BMP蛋白相互作用，參與正常胚胎的發(fā)育調(diào)控[48]。PAPAA蛋白的缺乏會(huì)下調(diào)小鼠卵巢類(lèi)固醇生成和損害女性生育能力[49]，其mRNA水平與妊娠期間血清水平分別是識(shí)別卵母細(xì)胞能正常發(fā)育至胚胎和評(píng)估子宮內(nèi)膜容受性的潛在重要標(biāo)志[50-51]。PAPAA蛋白具有成為預(yù)測(cè)高繁殖力個(gè)體診斷工具的潛力。和是參與精子發(fā)育的關(guān)鍵基因[52-53]。是一種木糖基轉(zhuǎn)移酶，參與對(duì)著床期囊胚的附著和生長(zhǎng)有直接作用的蛋白多糖肝素/硫酸肝素合成過(guò)程[54]。這些證據(jù)均表明了這11個(gè)繁殖調(diào)控基因可能是調(diào)控山羊多羔性狀的關(guān)鍵基因，在改良山羊繁殖性能方面具有較大的應(yīng)用價(jià)值，但其在產(chǎn)羔調(diào)節(jié)網(wǎng)絡(luò)中的具體功能還需進(jìn)一步研究。

A：PRLHR，B：MSX2，C：DRD1，D：ADCY10

LF組差異表達(dá)基因（和）在川中黑山羊排卵、植入后胚胎發(fā)育等生殖過(guò)程可能呈現(xiàn)負(fù)調(diào)控作用。研究表明，可抑制應(yīng)激條件引發(fā)的牛卵丘細(xì)胞的凋亡[55]，對(duì)牛卵母細(xì)胞發(fā)育可能具有增益作用。其5′側(cè)翼區(qū)雜合SNPs雖不影響其在卵巢中的表達(dá)，但對(duì)山羊的多產(chǎn)性具有顯著影響[56]，且編碼區(qū)同義突變還可通過(guò)調(diào)控基因表達(dá)影響與受體結(jié)合的親和力[57]，進(jìn)一步影響動(dòng)物繁殖。因此，后續(xù)擬定分析部分SNPs是否與山羊排卵率相關(guān)，對(duì)于解析川中黑山羊產(chǎn)羔性狀遺傳基礎(chǔ)可能具有理論指導(dǎo)意義。啟動(dòng)子區(qū)域內(nèi)的SNPs可下調(diào)其基因的表達(dá)[58]，進(jìn)而導(dǎo)致胎盤(pán)介導(dǎo)的妊娠并發(fā)癥及女性的反復(fù)流產(chǎn)[59-60]。EDNRA作為一種內(nèi)皮素受體，不僅參與山羊毛色調(diào)控[61]，還在肥胖引起的小鼠排卵缺失過(guò)程中發(fā)揮重要作用[62]。突變可破壞DNA損傷修復(fù)過(guò)程，可能是引發(fā)人類(lèi)卵巢早衰的潛在病因[63-64]，具體的生殖功能調(diào)控機(jī)制還有待探索?；蚓幋a速激肽家族成員P（SP）和神經(jīng)激肽A（Neurokinin A，NKA）[65]，參與調(diào)控哺乳動(dòng)物生殖激素分泌[66]。速激肽信號(hào)對(duì)排卵前促黃體素分泌激增具有促進(jìn)作用，也被認(rèn)為是調(diào)節(jié)排卵的一個(gè)負(fù)面因素，其功能障礙導(dǎo)致雌鼠排卵障礙[67]。此外，川中黑山羊LF組次優(yōu)選擇基因與美姑山羊的繁殖性狀相關(guān)[68]。這些發(fā)現(xiàn)均為產(chǎn)羔數(shù)調(diào)控網(wǎng)絡(luò)解析提供了新視角。

HF與LF組間存在多個(gè)共享窗口基因，正如之前嶗山奶山羊與大足黑山羊的相關(guān)報(bào)道[13-14]，川中黑山羊19個(gè)繁殖候選基因中也存在3個(gè)HF和LF組的交集基因（）。其中調(diào)控前列腺素水平[69]，參與睪酮的生成[70]，存在于卵母細(xì)胞和顆粒細(xì)胞中，體外成熟液中其激動(dòng)劑的添加可提高卵母細(xì)胞的囊胚率[71]，然用其激動(dòng)劑孵育精子，囊胚率的結(jié)果相反[72]。這些基因的功能與動(dòng)物繁殖緊密相關(guān)，其功能的深入挖掘?qū)ι窖蚍敝硨W(xué)研究是極為重要的，但該相關(guān)報(bào)道并未對(duì)共享窗口基因功能定位做深入剖析，這3個(gè)共享基因在川中黑山羊產(chǎn)羔性狀調(diào)控機(jī)制中發(fā)揮主控或微效調(diào)控作用還有待探查。

進(jìn)一步的功能分析表明，神經(jīng)活性配體-受體相互作用、糖胺聚糖-硫酸乙酰肝素/肝素的生物合成、鈣離子信號(hào)通路、cAMP信號(hào)通路和葉酸生物合成通路，在動(dòng)物繁殖活動(dòng)中也發(fā)揮重要功能。例如硫酸肝素蛋白聚糖是在GDF9信號(hào)通路中起著重要作用，并參與排卵周期卵泡中卵母細(xì)胞信號(hào)和卵丘細(xì)胞功能的形成[73]。鈣離子信號(hào)通路參與調(diào)控哺乳動(dòng)物卵母細(xì)胞發(fā)育[74-75]。cAMP通路調(diào)節(jié)卵巢類(lèi)固醇P4、E2的產(chǎn)生[76]，參與卵泡顆粒細(xì)胞分化和成熟[77]。神經(jīng)活性配體-受體相互作用途徑不僅調(diào)節(jié)家禽卵泡發(fā)育和卵子生產(chǎn)[78-79]，還參與軟體動(dòng)物性腺發(fā)育和產(chǎn)卵[80-81]。推測(cè)這些通路可能參與川中黑山羊卵母細(xì)胞發(fā)育及排卵過(guò)程，并且對(duì)排卵率至關(guān)重要。此外，葉酸(維生素B9)是機(jī)體健康和發(fā)育的最佳營(yíng)養(yǎng)必需物質(zhì)。葉酸缺乏會(huì)引發(fā)貧血、生殖健康和胎兒發(fā)育受損[82-83]。多個(gè)HF特異性強(qiáng)選擇候選基因在這些通路中顯著富集，暗示這些通路可能是川中黑山羊產(chǎn)羔數(shù)性狀關(guān)鍵調(diào)控途徑，對(duì)其多羔性狀調(diào)節(jié)起促進(jìn)作用。

本研究還特別分析19個(gè)候選基因外顯子區(qū)域的SNPs，以期進(jìn)一步篩選山羊多羔機(jī)制關(guān)鍵分子標(biāo)記。發(fā)現(xiàn)導(dǎo)致氨基酸序列變化的外顯子突變只在HF組中，并且僅定位于（G529A）與（A281T）候選基因中。測(cè)序結(jié)果表明，外顯子G529A與A281T突變?cè)谌斯みx擇條件下被強(qiáng)烈選擇，表現(xiàn)出較高的遺傳分化，改變了、的翻譯，可能在山羊繁殖力中發(fā)揮關(guān)鍵作用，有望作為高產(chǎn)標(biāo)記用于山羊分子育種。后續(xù)在川中黑山羊群體中擴(kuò)群驗(yàn)證G529A與A281T突變位點(diǎn)，將成為本課題組下一步工作的方向與重點(diǎn)。

總而言之，本研究在川中黑山羊群體中發(fā)現(xiàn)了許多新的變異體。這些可能是調(diào)控生殖的關(guān)鍵的基因。其中、、、、、和可能是調(diào)控山羊多羔性狀的關(guān)鍵基因，同時(shí)這些候選基因大部分顯著富集在神經(jīng)活性配體-受體相互作用、糖胺聚糖-硫酸乙酰肝素/肝素的生物合成、鈣離子信號(hào)通路、cAMP信號(hào)通路和葉酸生物合成等潛在的生殖相關(guān)通路中。外顯子G529A與突變可能是川中黑山羊多羔性狀的關(guān)鍵遺傳標(biāo)記，在改良山羊繁殖性能方面可能具有較大的應(yīng)用價(jià)值。

4 結(jié)論

本研究篩選到19個(gè)與產(chǎn)羔性狀相關(guān)的新候選基因，包括11個(gè)HF組特異性候選因（和），5個(gè)LF組特異性候選因（和）和3個(gè)LF組和HF組共享窗口基因（和）。

[1] DE LIMA L G, DE SOUZA N O B, RIOS R R, DE MELO B A, DOS SANTOS L T A, DE MORAES SILVA K, MURPHY T W, FRAGA A B. Advances in molecular genetic techniques applied to selection for litter size in goats (): a review. Journal of Applied Animal Research, 2020, 48(1): 38-44. doi:10.1080/09712119.2020.1717497.

[2] MULSANT P, LECERF F, FABRE S, SCHIBLER L, MONGET P, LANNELUC I, PISSELET C, RIQUET J, MONNIAUX D, CALLEBAUT I, CRIBIU E, THIMONIER J, TEYSSIER J, BODIN L, COGNIé Y, CHITOUR N, ELSEN J M. Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola Mérino ewes. Proceedings of the National Academy of Sciences of the United States of America, 2001, 98(9): 5104-5109. doi:10.1073/pnas.091577598.

[3] GALLOWAY S M, MCNATTY K P, CAMBRIDGE L M, LAITINEN M P E, JUENGEL J L, JOKIRANTA T S, MCLAREN R J, LUIRO K, DODDS K G, MONTGOMERY G W, BEATTIE A E, DAVIS G H, RITVOS O. Mutations in an oocyte-derived growth factor gene () cause increased ovulation rate and infertility in a dosage- sensitive manner. Nature Genetics, 2000, 25(3): 279-283. doi:10. 1038/77033.

[4] MARTINEZ-ROYO A, JURADO J J, SMULDERS J P, MARTí J I, ALABART J L, ROCHE A, FANTOVA E, BODIN L, MULSANT P, SERRANO M, FOLCH J, CALVO J H. A deletion in the bone morphogenetic protein 15 gene causes sterility and increased prolificacy in Rasa Aragonesa sheep. Animal Genetics, 2008, 39(3): 294-297. doi:10.1111/j.1365-2052.2008.01707.x.

[5] HANRAHAN J P, GREGAN S M, MULSANT P, MULLEN M, DAVIS G H, POWELL R, GALLOWAY S M. Mutations in the genes for oocyte-derived growth factors GDF9 and BMP15 are associated with both increased ovulation rate and sterility in Cambridge and belclare sheep (). Biology of Reproduction, 2004, 70(4): 900-909. doi:10.1095/biolreprod.103.023093.

[6] NICOL L, BISHOP S C, PONG-WONG R, BENDIXEN C, HOLM L E, RHIND S M, MCNEILLY A S. Homozygosity for a single base-pair mutation in the oocyte-specific GDF9 gene results in sterility in Thoka sheep. Reproduction (Cambridge, England), 2009, 138(6): 921-933. doi:10.1530/rep-09-0193.

[7] AHLAWAT S, SHARMA R, MAITRA A. Screening of indigenous goats for prolificacy associated DNA markers of sheep. Gene, 2013, 517(1): 128-131. doi:10.1016/j.gene.2012.12.015.

[8] TEJANGOOKEH H M, SHAHNEH A Z, ZAMIRI M J, DALIRI M, KOHRAM H, JAVAREMI A N. Study of BMP15 gene polymorphism in Iranian goats. African Journal of Biotechnology, 2009, 8(13), 2929-2932.

[9] SUPAKORN C, PRALOMKARN W. Sheep FecB gene polymorphism role in Thai meat goat proliferation rate//Proceedings of 9th World Congress Genetics Applied to Livestock Production, Leipizig Germany, 2010.

[10] HUA G H, CHEN S L, AI J T, YANG L G. None of polymorphism of ovine fecundity major genes FecB and FecX was tested in goat. Animal Reproduction Science, 2008, 108(3/4): 279-286. doi:10.1016/ j.anireprosci.2007.08.013.

[11] HE Y Q, MA X K, LIU X Y, ZHANG C X, LI J. Candidate genes polymorphism and its association to prolificacy in Chinese goats. Journal of Agricultural Science, 2010, 2(1): 88–92. doi:10.5539/jas. v2n1p88.

[12] 李恒, 字向東. 全基因組測(cè)序在山羊上的研究進(jìn)展. 中國(guó)畜牧雜志, 2021, 57(10): 29-34. doi:10.19556/j.0258-7033.20200910-01.

LI H, ZI X D. Research progress on whole-genome sequencing on goat. Chinese Journal of Animal Science, 2021, 57(10): 29-34. doi:10.19556/j.0258-7033.20200910-01. (in Chinese)

[13] LAI F N, ZHAI H L, CHENG M, MA J Y, CHENG S F, GE W, ZHANG G L, WANG J J, ZHANG R Q, WANG X, MIN L J, SONG J Z, SHEN W. Whole-genome scanning for the litter size trait associated genes and SNPs under selection in dairy goat (). Scientific Reports, 2016, 6: 38096. doi:10.1038/srep38096.

[14] GUANG-XIN E, ZHAO Y J, HUANG Y F. Selection signatures of litter size in Dazu black goats based on a whole genome sequencing mixed pools strategy. Molecular Biology Reports, 2019, 46(5): 5517-5523. doi:10.1007/s11033-019-04904-6.

[15] ISLAM R, LIU X X, GEBRESELASSIE G, ABIED A, MA Q, MA Y H. Genome-wide association analysis reveals the genetic locus for high reproduction trait in Chinese Arbas Cashmere goat. Genes & Genomics, 2020, 42(8): 893-899. doi:10.1007/s13258-020-00937-5.

[16] E G X, ZHOU D K, YANG B G, DUAN X H, NA R S, HAN Y G, ZENG Y. Association analysis of sixty-seven single nucleotide polymorphisms with litter size in Dazu Black goats. Animal Genetics, 2020, 51(1): 151-152. doi:10.1111/age.12879.

[17] WANG J J, ZHANG T, CHEN Q M, ZHANG R Q, LI L, CHENG S F, SHEN W, LEI C Z. Genomic signatures of selection associated with litter size trait in Jining gray goat. Frontiers in Genetics, 2020, 11: 286. doi:10.3389/fgene.2020.00286.

[18] WANG K, LIU X F, QI T, HUI Y Q, YAN H L, QU L, LAN X Y, PAN C Y. Whole-genome sequencing to identify candidate genes for litter size and to uncover the variant function in goats (). Genomics, 2021, 113(1): 142-150. doi:10.1016/j.ygeno.2020.11.024.

[19] ZI X D, MU X K, LU J Y, MA L, WANG Y. Polymorphisms of growth hormone(GH) and insulin-like growth factor I(IGF-I) genes in prolific Lezhi Black Goat: Possible association with litter size. Journal of Southwest University for Nationalities (Natural Science Edition), 2014, 40(3): 344-349.

[20] Lü M J, LI H, ZI X D. Assessment of estrous synchronization protocols and pregnancy specific protein B concentration for the prediction of kidding rate in Lezhi black goats. Small Ruminant Research, 2021, 195: 106299. doi:10.1016/j.smallrumres.2020.106299.

[21] YANG C X, ZI X D, WANG Y, YANG D Q, MA L, LU J Y, NIU H R, XIAO X. Cloning and mRNA expression levels of GDF9, BMP15, and BMPR1B genes in prolific and non-prolific goat breeds. Molecular Reproduction and Development, 2012, 79(1): 2. doi:10. 1002/mrd.21386.

[22] ZI X D, LU J Y, MA L. Identification and comparative analysis of the ovarian microRNAs of prolific and non-prolific goats during the follicular phase using high-throughput sequencing. Scientific Reports, 2017, 7: 1921. doi:10.1038/s41598-017-02225-x.

[23] ZI X D, LU J Y, ZHOU H, MA L, XIA W, XIONG X R, LAN D L, WU X H. Comparative analysis of ovarian transcriptomes between prolific and non-prolific goat breeds via high-throughput sequencing. Reproduction in Domestic Animals, 2018, 53(2): 344-351. doi:10. 1111/rda.13111.

[24] LI H, DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25(14): 1754- 1760. doi:10.1093/bioinformatics/btp324.

[25] BICKHART D M, ROSEN B D, KOREN S, SAYRE B L, HASTIE A R, CHAN S, LEE J, LAM E T, LIACHKO I, SULLIVAN S T, BURTON J N, HUSON H J, NYSTROM J C, KELLEY C M, HUTCHISON J L, ZHOU Y, SUN J J, CRISà A, PONCE DE LEóN F A, SCHWARTZ J C, HAMMOND J A, WALDBIESER G C, SCHROEDER S G, LIU G E, DUNHAM M J, SHENDURE J, SONSTEGARD T S, PHILLIPPY A M, VAN TASSELL C P, SMITH T P L. Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome. Nature Genetics, 2017, 49(4): 643-650. doi:10.1038/ng.3802.

[26] LI H, HANDSAKER B, WYSOKER A, FENNELL T, RUAN J, HOMER N, MARTH G, ABECASIS G, DURBIN R. 1000 genome project data processing subgroup. The sequence alignment/ map format and SAMtools. Microbiology Spectrum, 2009, 25(16): 2078-2079. doi:10.1093/bioinformatics/btp352.

[27] MCKENNA A, HANNA M, BANKS E, SIVACHENKO A, CIBULSKIS K, KERNYTSKY A, GARIMELLA K, ALTSHULER D, GABRIEL S, DALY M, DEPRISTO M A. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Cell Reports, 2010, 20(9): 1297-1303. doi:10.1101/ gr.107524.110.

[28] YANG H, WANG K. Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR. Nature Protocols, 2015, 10(10): 1556-1566. doi:10.1038/nprot.2015.105.

[29] RAUDVERE U, KOLBERG L, KUZMIN I, ARAK T, ADLER P, PETERSON H, VILO J. G: Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Research, 2019, 47(W1): W191-W198. doi:10.1093/ nar/gkz369.

[30] BU D C, LUO H T, HUO P P, WANG Z H, ZHANG S, HE Z H, WU Y, ZHAO L H, LIU J J, GUO J C, FANG S S, CAO W C, YI L, ZHAO Y, KONG L. KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Research, 2021, 49(W1): W317-W325. doi:10.1093/nar/gkab447.

[31] MOKHTARI M S, ASADI FOZI M, GUTIERREZ J P, NOTTER D R. Genetic and phenotypic aspects of early reproductive performance in Raeini Cashmere goats. Tropical Animal Health and Production, 2019, 51(8): 2175-2180. doi:10.1007/s11250-019-01915-0.

[32] ?URI?I? D, BENI? M, ?AJA I ?, VALPOTI? H, SAMARD?IJA M. Influence of season, rainfall and air temperature on the reproductive efficiency in Romanov sheep in Croatia. International Journal of Biometeorology, 2019, 63(6): 817-824. doi:10.1007/s00484-019- 01696-z.

[33] ASTUTI D A, KHOTIJAH L, MAIDIN M S, NUGROHO P. Reproductive profile of etawah crossbred does fed Flushing diet containing different kinds of plant oil and animal fat. Pakistan Journal of Biological Sciences, 2020, 23(5): 650-657. doi:10.3923/pjbs.2020. 650.657.

[34] MCCARTHY M I, ABECASIS G R, CARDON L R, GOLDSTEIN D B, LITTLE J, IOANNIDIS J P A, HIRSCHHORN J N. Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nature Reviews Genetics, 2008, 9(5): 356-369. doi:10. 1038/nrg2344.

[35] HONG E P, PARK J W. Sample size and statistical power calculation in genetic association studies. Genomics & Informatics, 2012, 10(2): 117-122. doi:10.5808/gi.2012.10.2.117.

[36] LIU N, CUI W B, CHEN M Y, ZHANG X L, SONG X Y, PAN C Y. A 21-bp indel within thegene is significantly associated with litter size in goat. Animal Biotechnology, 2021, 32(2): 213-218. doi:10.1080/10495398.2019.1677682.

[37] JIANG E H, KANG Z H, WANG X Y, LIU Y, LIU X F, WANG Z, LI X C, LAN X Y. Detection of insertions/deletions (InDels) within the goat Runx2 gene and their association with litter size and growth traits. Animal Biotechnology, 2021, 32(2): 169-177. doi:10.1080/10495398. 2019.1671858.

[38] BALBACH M, FUSHIMI M, HUGGINS D J, STEEGBORN C, MEINKE P T, LEVIN L R, BUCK J. Optimization of lead compounds into on-demand, nonhormonal contraceptives: Leveraging a public- private drug discovery institute collaboration. Biology of Reproduction, 2020, 103(2): 176-182. doi:10.1093/biolre/ioaa052.

[39] CHEN H, CHAN H C. Amplification of FSH signalling by CFTR and nuclear soluble adenylyl cyclase in the ovary. Clinical and Experimental Pharmacology & Physiology, 2017, 44(Suppl 1): 78-85. doi:10.1111/1440-1681.12756.

[40] JAYARAJAN V, APPUKUTTAN A, ASLAM M, REUSCH P, REGITZ-ZAGROSEK V, LADILOV Y. Regulation of AMPK activity by type 10 adenylyl cyclase: Contribution to the mitochondrial biology, cellular redox and energy homeostasis. Cellular and Molecular Life Sciences, 2019, 76(24): 4945-4959. doi:10.1007/s00018-019-03152-y.

[41] WANG C, LI S J, LI C, FENG Y P, PENG X L, GONG Y Z. Molecular cloning, expression profile, polymorphism and the genetic effects of the dopamine D1 receptor gene on duck reproductive traits. Molecular Biology Reports, 2012, 39(9): 9239-9246. doi:10.1007/ s11033-012-1797-3.

[42] LIU Z, YANG N, YAN Y, LI G, LIU A, WU G, SUN C. Genome-wide association analysis of egg production performance in chickens across the whole laying period. BMC Genetics, 2019, 20(1): 67. doi:10.1186/ s12863-019-0771-7.

[43] BARONCHELLI S, VILLA N, REDAELLI S, LISSONI S, SACCHERI F, PANZERI E, CONCONI D, BENTIVEGNA A, CROSTI F, SALA E, BERTOLA F, MAROZZI A, PEDICINI A, VENTRUTO M, POLICE M A, DALPRà L. Investigating the role of X chromosome breakpoints in premature ovarian failure. Molecular Cytogenetics, 2012, 5(1): 32. doi:10.1186/1755-8166-5-32.

[44] FESTA A, UMANO G R, MIRAGLIA DEL GIUDICE E, GRANDONE A. Genetic evaluation of patients with delayed puberty and congenital hypogonadotropic hypogonadism: Is it worthy of consideration? Frontiers in Endocrinology, 2020, 11: 253. doi:10.3389/fendo.2020. 00253.

[45] LI J, LIU J, CAMPANILE G, PLASTOW G, ZHANG C, WANG Z, CASSANDRO M, GASPARRINI B, SALZANO A, HUA G, LIANG A, YANG L. Novel insights into the genetic basis of buffalo reproductive performance. BMC Genomics, 2018, 19(1): 814. doi:10. 1186/s12864-018-5208-6.

[46] NALLASAMY S, KAYA OKUR H S, BHURKE A, DAVILA J, LI Q X, YOUNG S L, TAYLOR R N, BAGCHI M K, BAGCHI I C. Msx homeobox genes act downstream of BMP2 to regulate endometrial decidualization in mice and in humans. Endocrinology, 2019, 160(7): 1631-1644. doi:10.1210/en.2019-00131.

[47] LIU Z K, WANG R C, HAN B C, YANG Y, PENG J P. A novel role of IGFBP7in mouse uterus: regulating uterine receptivity through Th1/Th2 lymphocyte balance and decidualization. PLoS ONE, 2012, 7(9): e45224. doi:10.1371/journal.pone.0045224.

[48] GERHART J, SCHEINFELD V L, MILITO T, PFAUTZ J, NEELY C, FISHER-VANCE D, SUTTER K, CRAWFORD M, KNUDSEN K, GEORGE-WEINSTEIN M. Myo/Nog cell regulation of bone morphogenetic protein signaling in the blastocyst is essential for normal morphogenesis and striated muscle lineage specification. Developmental Biology, 2011, 359(1): 12-25. doi:10.1016/j.ydbio. 2011.08.007.

[49] NYEGAARD M, OVERGAARD M T, SU Y Q, HAMILTON A E, KWINTKIEWICZ J, HSIEH M, NAYAK N R, CONTI M, CONOVER C A, GIUDICE L C. Lack of functional pregnancy- associated plasma protein-A (PAPPA) compromises mouse ovarian steroidogenesis and female Fertility1. Biology of Reproduction, 2010, 82(6): 1129-1138. doi:10.1095/biolreprod.109.079517.

[50] KORDUS R J, HOSSAIN A, CORSO M C, CHAKRABORTY H, WHITMAN-ELIA G F, LAVOIE H A. Cumulus cell pappalysin-1, luteinizing hormone/choriogonadotropin receptor, amphiregulin and hydroxy-delta-5-steroid dehydrogenase, 3 beta- and steroid delta- isomerase 1 mRNA levels associate with oocyte developmental competence and embryo outcomes. Journal of Assisted Reproduction and Genetics, 2019, 36(7): 1457-1469. doi:10.1007/s10815-019- 01489-8.

[51] YU M, WANG J, LIU S, WANG X Q, YAN Q. Novel function of pregnancy-associated plasma protein A: promotes endometrium receptivity by up-regulating N-fucosylation. Scientific Reports, 2017, 7: 5315. doi:10.1038/s41598-017-04735-0.

[52] WON J, DE EVSIKOVA C M, SMITH R S, HICKS W L, EDWARDS M M, LONGO-GUESS C, LI T S, NAGGERT J K, NISHINA P M. NPHP4is necessary for normal photoreceptor ribbon synapse maintenance and outer segment formation, and for sperm development. Human Molecular Genetics, 2010, 20(3): 482-496. doi:10.1093/hmg/ ddq494.

[53] MAO S, WU F, CAO X, HE M, LIU N, WU H, YANG Z, DING Q, WANG X. Tdrp deficiency contributes to low sperm motility and is a potential risk factor for male infertility. American Journal of Translational Research, 2016, 8(1): 177-187.

[54] FARACH M C, TANG J P, DECKER G L, CARSON D D. Heparin/heparan sulfate is involved in attachment and spreading of mouse embryos. Developmental Biology, 1987, 123(2): 401-410. doi:10.1016/0012-1606(87)90398-8.

[55] 潘陽(yáng)陽(yáng), 王萌, 芮弦, 王立斌, 何翃閎, 王靖雷, 馬睿, 徐庚全, 崔燕, 樊江峰, 余四九. IGF-1調(diào)控RBM3表達(dá)抑制低溫應(yīng)激誘導(dǎo)牦牛卵丘細(xì)胞凋亡. 中國(guó)農(nóng)業(yè)科學(xué), 2020, 53(11): 2285-2296.

PAN Y Y, WANG M, RUI X, WANG L B, HE H H, WANG J L, MA R, XU G Q, CUI Y, FAN J F, YU S J. RNA-binding motif protein 3(RBM3) expression is regulated by insulin-like growth factor(IGF-1) for protecting yak() cumulus cells from apoptosis during hypothermia stress. Scientia Agricultura Sinica, 2020, 53(11): 2285-2296. (in Chinese)

[56] THOMAS N, VENKATACHALAPATHY T, ARAVINDAKSHAN T, RAGHAVAN K C. Molecular cloning, SNP detection and association analysis of 5' flanking region of the goat IGF1gene with prolificacy. Animal Reproduction Science, 2016, 167: 8-15. doi:10.1016/j. anireprosci.2016.01.016.

[57] CHENG Y Y, LIU S C, WANG G, WEI W Z, HUANG S, YANG R, GENG H W, LI H Y, SONG J, SUN L D, YU H, HAO L L. Porcine IGF1synonymous mutation alter gene expression and protein binding affinity with IGF1R. International Journal of Biological Macromolecules, 2018, 116: 23-30. doi:10.1016/j.ijbiomac.2018.05.022.

[58] INAGAKI H, OTA S, NISHIZAWA H, MIYAMURA H, NAKAHIRA K, SUZUKI M, NISHIYAMA S, KATO T, YANAGIHARA I, KURAHASHI H. Obstetric complication-associatedpromoter polymorphisms may affect gene expression via DNA secondary structures. Journal of Human Genetics, 2019, 64(5): 459-466. doi:10.1038/s10038-019-0578-4.

[59] ARANDA F, UDRY S, PERéS WINGEYER S, AMSHOFF L C, BOGDANOVA N, WIEACKER P, LATINO J O, MARKOFF A, LARRA?AGA G. Maternal carriers of the ANXA5 M2 haplotype are exposed to a greater risk for placenta-mediated pregnancy complications. Journal of Assisted Reproduction and Genetics, 2018, 35(5): 921-928. doi:10.1007/s10815-018-1142-4.

[60] DRYLLIS G, GIANNOPOULOS A, ZOI C, POULIAKIS A, LOGOTHETIS E, VOULGARELIS M, ZOI K, KOUSKOUNI E, DINOU A, STAVROPOULOS-GIOKAS C, KREATSAS G, KONSTANTOPOULOS K, POLITOU M. Correlation of single nucleotide polymorphisms in the promoter region of the(annexin A5) gene with recurrent miscarriages in women of Greek origin. The Journal of Maternal-Fetal & Neonatal Medicine, 2020, 33(9): 1538-1543. doi:10.1080/14767058.2018.1521799.

[61] DI GERLANDO R, MASTRANGELO S, MOSCARELLI A, TOLONE M, SUTERA A M, PORTOLANO B, SARDINA M T. Genomic structural diversity in local goats: Analysis of copy-number variations. Animals, 2020, 10(6): E1040. doi:10.3390/ani10061040.

[62] HOHOS N M, ELLIOTT E M, GIORNAZI A, SILVA E, RICE J D, SKAZNIK-WIKIEL M E. High-fat diet induces an ovulatory defect associated with dysregulated endothelin-2 in mice. Reproduction (Cambridge, England), 2021, 161(3): 307-317. doi:10.1530/rep-20- 0290.

[63] YANG Y, GUO T, LIU R, KE H, XU W, ZHAO S, QIN Y.gene mutations in premature ovarian insufficiency. Human Mutation, 2020, 41(5): 1033-1041. doi:10.1002/humu.23997.

[64] YANG Y, ZHAO S, QIN Y. Response to “Should FANCL heterozygous pathogenic variants be considered as potentially causative of primary ovarian insufficiency? ”. Human Mutation, 2020, 41(9): 1700-1701. doi:10.1002/humu.24073.

[65] FERGANI C, NAVARRO V M. Expanding the role of tachykinins in the neuroendocrine control of reproduction. Reproduction (Cambridge, England), 2016, 153(1): R1-R14. doi:10.1530/rep-16-0378.

[66] LEóN S, FERGANI C, TALBI R, SIMAVLI S, MAGUIRE C A, GERUTSHANG A, NAVARRO V M. Characterization of the role of NKA in the control of puberty onset and gonadotropin release in the female mouse. Endocrinology, 2019, 160(10): 2453-2463. doi:10. 1210/en.2019-00195.

[67] LEóN S, FERGANI C, TALBI R, MAGUIRE C A, GERUTSHANG A, SEMINARA S B, NAVARRO V M. Tachykinin signaling is required for induction of the preovulatory luteinizing hormone surge and normal luteinizing hormone pulses. Neuroendocrinology, 2021, 111(6): 542-554. doi:10.1159/000509222.

[68] GUO J Z, TAO H X, LI P F, LI L, ZHONG T, WANG L J, MA J Y, CHEN X Y, SONG T Z, ZHANG H P. Whole-genome sequencing reveals selection signatures associated with important traits in six goat breeds. Scientific Reports, 2018, 8: 10405. doi:10.1038/s41598-018- 28719-w.

[69] PASTEL E, POINTUD J C, LOUBEAU G, DANI C, SLIM K, MARTIN G, VOLAT F, SAHUT-BARNOLA I, VAL P, MARTINEZ A, LEFRAN?OIS-MARTINEZ A M. Aldose reductases influence prostaglandin F2α levels and adipocyte differentiation in male mouse and human species. Endocrinology, 2015, 156(5): 1671-1684. doi:10. 1210/en.2014-1750.

[70] ZHANG Q, PEI L G, LIU M, LV F, CHEN G H, WANG H. Reduced testicular steroidogenesis in rat offspring by prenatal nicotine exposure: Epigenetic programming and heritability via nAChR/ HDAC4. Food and Chemical Toxicology, 2020, 135: 111057. doi:10.1016/j.fct.2019.111057.

[71] OLABARRIETA E, TOTORIKAGUENA L, AGIRREGOITIA N, AGIRREGOITIA E. Implication of mu opioid receptor in thematuration of oocytes and its effects on subsequent fertilization and embryo development in mice. Molecular Reproduction and Development, 2019, 86(9): 1236-1244. doi:10.1002/mrd.23248.

[72] OLABARRIETA E, TOTORIKAGUENA L, ROMERO- AGUIRREGOMEZCORTA J, AGIRREGOITIA N, AGIRREGOITIA E. Mu opioid receptor expression and localisation in murine spermatozoa and its role in IVF. Reproduction Fertility and Development, 2020, 32(4): 349-354. doi:10.1071/rd19176.

[73] WATSON L N, MOTTERSHEAD D G, DUNNING K R, ROBKER R L, GILCHRIST R B, RUSSELL D L. Heparan sulfate proteoglycans regulate responses to oocyte paracrine signals in ovarian follicle morphogenesis. Endocrinology, 2012, 153(9): 4544-4555. doi:10. 1210/en.2012-1181.

[74] TIWARI M, PRASAD S, SHRIVASTAV T G, CHAUBE S K. Calcium signaling during meiotic cell cycle regulation and apoptosis in mammalian oocytes. Journal of Cellular Physiology, 2017, 232(5): 976-981. doi:10.1002/jcp.25670.

[75] STEWART T A, DAVIS F M. An element for development: Calcium signaling in mammalian reproduction and development. Biochimica et Biophysica Acta (BBA) - Molecular Cell Research, 2019, 1866(7): 1230-1238. doi:10.1016/j.bbamcr.2019.02.016.

[76] LU N S, LI M J, LEI H L, JIANG X Y, TU W L, LU Y, XIA D. Butyric acid regulates progesterone and estradiol secretion via cAMP signaling pathway in porcine granulosa cells. The Journal of Steroid Biochemistry and Molecular Biology, 2017, 172: 89-97. doi:10. 1016/j.jsbmb.2017.06.004.

[77] JOZKOWIAK M, HUTCHINGS G, JANKOWSKI M, KULCENTY K, MOZDZIAK P, KEMPISTY B, SPACZYNSKI R Z, PIOTROWSKA- KEMPISTY H. The stemness of human ovarian granulosa cells and the role of resveratrol in the differentiation of MSCs-A review based on cellular and molecular knowledge. Cells, 2020, 9(6): E1418. doi:10.3390/cells9061418.

[78] ZHANG T, CHEN L, HAN K P, ZHANG X Q, ZHANG G X, DAI G J, WANG J Y, XIE K Z. Transcriptome analysis of ovary in relatively greater and lesser egg producing Jinghai Yellow Chicken. Animal Reproduction Science, 2019, 208: 106114. doi:10.1016/j.anireprosci. 2019.106114.

[79] CHEN X, SUN X, CHIMBAKA I M, QIN N, XU X, LISWANISO S, XU R, GONZALEZ J M. Transcriptome analysis of ovarian follicles reveals potential pivotal genes associated with increased and decreased rates of chicken egg production. Frontiers in Genetics, 2021, 12: 622751. doi:10.3389/fgene.2021.622751.

[80] XU R Y, PAN L Q, YANG Y Y, ZHOU Y Y. Characterizing transcriptome in female scallopfarreri provides new insights into the molecular mechanisms of reproductive regulation during ovarian development and spawn. Gene, 2020, 758: 144967. doi:10. 1016/j.gene.2020.144967.

[81] HUANG D X, ZHANG B, HAN T, LIU G B, CHEN X, ZHAO Z H, FENG J Q, YANG J W, WANG T M. Genome-wide prediction and comparative transcriptomic analysis reveals the G protein-coupled receptors involved in gonadal development of. Genomics, 2021, 113(1): 967-978. doi:10.1016/j.ygeno.2020.10.030.

[82] NADERI N, HOUSE J D. Recent developments in folate nutrition// Advances in Food and Nutrition Research. Amsterdam: Elsevier, 2018: 195-213. doi:10.1016/bs.afnr.2017.12.006.

[83] BROWN L L, COHEN B E, EDWARDS E, GUSTIN C E, NOREEN Z. Physiological need for calcium, iron, and folic acid for women of various subpopulations during pregnancy and beyond. Journal of Womens Health (Larchmt), 2021, 30(2): 207-211. doi:10.1089/jwh. 2020.8873.

Screening of Key Regulatory Genes for Litter Size Trait Based on Whole Genome Re-Sequencing in Goats ()

LI Heng1, ZI XiangDong1, WANG Hui2, XIONG Yan1, Lü MingJie1,LIU Yu1, JIANG XuDong1

1Key Laboratory of Animal Science of National Ethnic Affairs Commission, Southwest Minzu University, Chengdu 610041;2Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Ministry of Education, Southwest Minzu University, Chengdu 610041

【Objective】The purpose of this study was to analyze the genome of different fecundity populations of goats (s) and to explore the key regulatory genes involved in the regulation of litter size traits of Chuanzhong black goats (CBGs), and to provide the theoretical reference for analyzing the genetic mechanism of litter size traits and molecular genetic improvement of fecundity in goats. 【Method】The high fecundity (HF) CBG does (= 6) that produced 4-6 kids per doe kidding and low fecundity (LF) does (= 6) that produced only one kid per doe kidding were chosen in this study. The jugular blood samples were collected to extract genomic DNA. The 350 bp double-terminal sequencing library was constructed, and then 12 whole genome libraries were resequenced by IlluminaHiSeqPE150 platform. The clean data from sequencing were mapped to goat reference genome ARS1 by using BWA software, and two whole-genome scanning analysis methods (and) were used to comprehensively analyze the high-quality SNPs obtained to identify candidate regions. GO analysis and KEGG pathway analysis were performed on the G:Profiler and KOBAS online databases, respectively, to screen candidate genes for regulating the number of kids in CBGs. To further identify the key genetic markers that regulate the number of kids, the synonymous and non-synonymous single nucleotide polymorphisms (SNPs) of reproductive candidate genes were mapped and screened according to the variation analysis report of genome resequencing. The amplified products of 12 goat samples were sequenced by Sanger sequencing to verify the resequencing results.【Result】A total of 431.50 Gb clean data were obtained from the genome resequencing study of 12 CBGs. Through mutation detection and annotation, 7 771 417 SNPs were detected in HF group and 8 935 907 SNPs were detected in LF group, and all types of the LF group SNPs were more than those in HF group. The windows that reach the maximumvalue of top 5% and the minimumvalue of top 5% were set as candidate regions. A total of 130 strong selection signals were annotated in the regions with low heterozygosity and high genetic differentiation, of which 84, 59 and 13 genes were annotated in HF group, LF group and shared window, respectively. GO enrichment analysis and KEGG pathway showed that 19 candidate genes were involved in the regulation of reproduction, reproduction and embryonic development of CBG, including 11 HF group-specific candidate genes (,,,,,,,,,, and), and five strong selection signal genes (,,,, and) in LF group, and three window genes (,and) in HF group shared with LF group. The most GO terms, such as G-protein-coupled receptor activity, hormone response and neuropeptide signal pathway, contained these 19 candidate genes. In addition, nine of the 14 HF candidate genes were significantly enriched in metabolic pathway, neuroactive ligand-receptor interaction, glycosaminoglycan-heparan sulfate/heparin biosynthesis, calcium signal pathway, cAMP signal pathway and folate biosynthesis KEGG pathways (＜0.05). Among the 19 reproductive candidate genes, there were two synonymous mutations (G771T,G529,,andgene mutations could be detected, and this result was consistent with the results of genome resequencing, in whichG529A polymorphism led to alanine mutation to threonine, and【Conclusion】A total of 11 HF group-specific candidate genes were found in this study, which were speculated to be the key regulatory genes for fecundity trait. The mutations ofgene exon G529A andexon A281T might be the key genetic markers for regulating prolificacy traits in goats, which had great application value in improving reproductive performance of goats.

Chuanzhong black goat; genome resequencing; fecundity; candidate genes

10.3864/j.issn.0578-1752.2022.23.015

2021-08-17；

2022-10-12

國(guó)家自然科學(xué)基金（31902154）、西南民族大學(xué)中央高校基本科研業(yè)務(wù)專(zhuān)項(xiàng)（2021PTJS26）

李恒，E-mail：lih199501 @sina.com。通信作者字向東，E-mail：zixd@sina.com

（責(zé)任編輯 林鑒非）