基于全基因组重测序筛选关中奶山羊产羔数性状相关候选基因研究

任毅杰,贤 明,倪 洁,李 宇,郭松茂,冉本康,胡建宏

(西北农林科技大学 动物科技学院,陕西 杨凌 712100)

繁殖性状作为重要的经济性状直接反映了母羊的生产能力,是影响山羊养殖业发展的重要因素。山羊的繁殖性状是由多基因、多位点控制并受多种因素影响的复杂数量性状,其涉及卵巢卵泡发育、卵母细胞成熟、排卵、受精、胚胎发育等多个生物过程,且该性状在生产中存在遗传力较低和育种进展缓慢等问题[1]。目前,研究人员已发现有些基因能够影响母羊的繁殖性能,如BMPR1B、GDF9和BMP15等基因通过调控下丘脑-垂体-性腺轴的生殖激素分泌、排卵率、胚胎存活率和子宫内膜容受性等生物过程,进而影响母羊繁殖力[1-4]。最近,基于山羊全基因重测序的相关研究发现,KIT、PAK1、KCNH7、SMAD9和PRKKA1等基因很可能是影响济宁青山羊产羔数性状的关键基因,而AURKA、ENDOG、SOX2、RORA、GJA10、RXFP2、CDC25C和NANOS3等是陕北白绒山羊产羔数性状的候选基因[5-6]。

关中奶山羊作为中国著名优良乳用奶山羊品种,其适应性强、耐粗饲、抗病力强,产奶性能稳定、产奶量高、产奶周期长,体格健壮、结构匀称,具有身长、颈长和腿长“三长”特征[7]。关中奶山羊母羊头胎多产单羔,二胎及二胎以上多产双羔或多羔,平均产羔率为188%[8]。研究发现,INHA、PRLR、NGF、KISS1、KITLG和KIT等基因的多态性与关中奶山羊产羔数具有显著性关系[9-14]。An等[15]通过对多胎和单胎关中奶山羊发情期卵巢中的miRNAs进行鉴定和分析,发现miRNAs可能通过降低GnRH信号通路和卵巢类固醇激素合成的靶基因表达来调节卵泡发育和排卵,且ggo-miR-4488-p3_1ss10CG、bta-miR-2892-p5_1ss8CG和hsa-miR-4532_L+1R-3与多胎性状密切相关。张效川[16]通过研究关中奶山羊发情期和间情期的卵巢全基因组甲基化差异,发现NR5A2、STAR、FGF2和BMP5基因可能调控关中奶山羊的发情行为。除发情行为外,母羊妊娠阶段子宫内的变化也可影响其产羔数。妊娠的维持取决于母体识别妊娠和胚胎植入期间基因表达的动态变化,然而在妊娠建立的敏感时期,子宫内膜基因的异常表达能够导致胚胎植入失败和不孕。Song等[17]通过使用全基因组亚硫酸氢盐测序和生物信息学分析研究了关中奶山羊妊娠第5天和第15天子宫内膜的全基因组DNA甲基化模式,发现了一些与胚胎的子宫容受性相关的差异甲基化基因,如IGF2BP2、PTGDS、VEGFB、PGR、IGFBP3和IGF1R等。尽管目前已有部分研究探究了调控关中奶山羊繁殖相关性状的基因和机理,但迄今仍未发现决定关中奶山羊繁殖性状的主效基因和分子标记,且尚未有研究从全基因组层面挖掘关中奶山羊产羔数性状的候选基因。

为了阐明关中奶山羊产羔数性状的遗传基础,本研究根据生产记录将关中奶山羊母羊分为高繁殖力组和低繁殖力组两个不同产羔数群体,分析不同产羔数群体的变异模式并进行选择消除分析,筛选与关中奶山羊产羔数性状相关的候选基因和变异位点,以期为关中奶山羊群体内优质高繁个体的选育以及未来关中奶山羊的生物学和遗传育种等研究提供一定的试验支撑。

1 材料与方法

1.1 试验动物

从陕西澳尼克奶山羊育种有限公司选取16只产奶状况正常、乳品质良好、无流产、难产、死胎情况的母羊,根据母羊的前两胎产羔记录将其分为低繁殖力组(Low fecundity, LF)和高繁殖力组(High fecundity, HF),表型统计及分组信息见表1。

表1 关中奶山羊表型统计及分组信息Table 1 Phenotypic statistics and grouping information of Guanzhong dairy goats

1.2 样品采集及测序

使用一次性真空采血管(EDTA K2)采集16只关中奶山羊的颈静脉血,轻轻颠倒使抗凝剂和血液混匀后,将抗凝血转移至NEST冻存管并通过液氮保存。样品运回实验室后干冰保温寄送至北京诺禾致源科技股份有限公司,对全血样品进行基因组DNA的提取、质量检测和文库构建,通过Illumina HiSeq平台对构建好的文库进行测序。

1.3 变异检测及注释

根据比对结果,采用SAMTOOLS(mpileup -m 2 -F 0.002)进行群体SNP的检测。为了降低SNP检测的错误率,以SNP的reads支持数不低于4和SNP的质量值不低于20为标准进行过滤。利用ANNOVAR软件对检测出的基因组变异进行注释。

1.4 选择信号检测

利用滑动窗口法对关中奶山羊LF组和HF组进行全基因组扫描,采用核苷酸多样性(Pi或θπ)、Tajima's D检测和Fst &θπ3种方法进行选择信号检测,Fst &θπ联合筛选的top 5%区域为受选择区域。其中,HF组受选择区域为log2(Pi ratio(Pi LF/Pi HF))>0.56,且Fst>0.12;LF组受选择区域为log2(Piratio (PiHF/PiLF))>0.42,且Fst>0.13。

1.5 基因功能富集分析

由于山羊富集分析数据库较少,因此将通过g:Profiler在线工具将山羊基因symbol转换为小鼠同源基因后再进行基因功能的富集分析;使用DAVID和KOBAS对HF组和LF组的受选择基因分别进行GO (gene ontology)功能注释和KEGG富集分析。

2 结果与分析

2.1 全基因组重测序数据质控与比对

本研究共获得原始数据(raw bases)494.82 Gb,过滤后的纯化数据(clean bases)491.16 Gb,测序质量高,GC分布正常。所有样本的比对率在99.67%～99.76%之间,对参考基因组的平均覆盖深度在7.90×～9.34×之间,1×覆盖度(至少有1个碱基的覆盖)在94.86%以上,4×覆盖度(至少有4个碱基的覆盖)在85.49%以上。

2.2 不同繁殖力组的基因组变异模式

关中奶山羊不同繁殖力组的基因组变异模式见图1,LF组的平均测序深度略高于HF组(图1A);从突变信息看,在LF组和HF组中分别鉴定出17492669和17742455个SNPs位点,与LF组相比,HF组拥有更多的变异位点(图1B)。随着产羔数的增加,LF组和HF组的SNP杂合与纯合位点的突变比率逐渐下降(图1C)。

2.3 不同繁殖力组受选择区域筛选

本研究在HF组和LF组分别筛选到1 754和865个选择信号窗口(图2A、图2C),其中HF组中最多有227个(12.94%)选择窗口位于9号染色体,29号染色体仅筛选到6个(0.34%)选择信号窗口(图2B);LF组中最多有95个(12.94%)选择窗口位于2号染色体,25号染色体仅筛选到5个(0.58%)选择信号窗口(图2D)。同时,对HF组和LF组选择信号窗口的基因进行注释,分别注释到437和230个基因,HF组中与繁殖相关的基因有SP1、FGFR1、AMHR2、SMAD9、OXTR、ACTG1、PCGF1、STAT3、ACVR1B、JAK2、SGK1、LTBP1、ACVR1B、CREBRF和OPRK1等;LF组中与繁殖相关的基因有TGFBR3、BMPR2、SOX11、BCL2、THBS1、ITGA11、SOX17、ROCK2、FZD7、EIF2AK3、MAP2K5、CASP8、ILK、HESX1和LIFR等。此外,除HF组和LF组特有的受选择基因外,两组的受选择区域有8个共同的受选择基因,分别为ABI3BP、ST5、TAOK3、SERPINB11、PGBD5、C1GALT1、GABRG2和PRR16(图3)。

图1 不同繁殖力组关中奶山羊的基因组变异模式Fig. 1 Genomic variation patterns of Guanzhong dairy goats with different fecundity groups

图 2 关中奶山羊不同繁殖力组基于Fst &θπ筛选的受选择区域A. HF组受选择区域;B. HF组基于Fst &θπ筛选的1 754个选择信号窗口在染色体上的分布;C. LF组受选择区域;D. LF组基于Fst &θπ筛选的865个选择信号窗口在染色体上的分布Fig. 2 The selected regions of Guanzhong dairy goats with different litter sizes screened based on Fst &θπA. Selected region of HF group; B. The chromosome distribution of 1754 selective signal Windows screened by Fst &θπ in HF group;C. Selected region of LF group; D. The chromosome distribution of 865 selective signal Windows screened by Fst &θπ in LF group

图3 HF组和LF组受选择区域基因韦恩图Fig. 3 Venn plots of selected regional genes in HF and LF groups

2.4 候选基因的Tajima's D检验与核苷酸多态性检验

对HF组与LF组中的部分基因(HF组:SP1和FGFR1;LF组:TGFBR3和BMPR2)进行Tajima's D检验与核苷酸多态性检验,结果如图4所示。SP1、FGFR1、TGFBR3和BMPR2等4个基因在

图4 HF组和LF组产羔数性状候选基因的Tajima's D值和核苷酸多样性值分布灰色区域为受选择区域,红线为HF组,蓝线为LF组Fig. 4 The distribution of Tajima's D value and nucleotide diversity of candidate genes for litter size traits inHF and LF groups Gray area is selected region, red line is HF group, blue line is LF group

HF组与LF组中均受到不同程度的选择。其中,在SP1和FGFR1基因的受选择区域,LF组中的Tajima's D值与核苷酸多态性均高于HF组(图4A、图4B);在TGFBR3和BMPR2基因的受选择区域,HF组中的Tajima's D值与核苷酸多态性均高于LF组(图4C、4D)。

2.5 不同繁殖力组受选择基因的功能富集

HF组的GO富集分析结果显示,生物学过程分析在与物质运输和蛋白质磷酸化的相关条目中具有较多基因;细胞成分分析在涵盖细胞胞内区域的相关条目中具有较多基因;分子功能分析在与蛋白或核苷酸结合和酶活性的相关条目中具有较多基因(图5A)。HF组受选择基因显著富集到与繁殖性能相关的通路,如MAPK信号通路、TGF-β信号通路、调节干细胞多能性的信号通路、黏附、细胞凋亡、mTOR信号通路、Jak-STAT信号通路和催产素信号通路等(图5B)。LF组的GO富集分析结果显示,生物学过程分析在与蛋白质磷酸化和转运的相关条目中具有较多基因;细胞成分分析在涵盖细胞胞内区域的相关条目中具有较多基因;分子功能分析在与蛋白或核苷酸结合和转移酶活性的相关条目中具有较多基因(图5C)。LF组受选择基因显著富集到与繁殖性能相关的通路,如PI3K-Akt信号通路、内质网中的蛋白质加工、p53信号通路、调节干细胞多能性的信号通路、ECM-受体相互作用和Wnt信号通路等(图5D)。

图5 HF组和LF组受选择基因富集分析A. HF组的前30条GO条目条形图;B. HF组的前20个KEGG通路气泡图;C. LF组的前30条GO条目条形图;D. LF组的前20个KEGG通路气泡图Fig. 5 Enrichment analysis of selected genes in HF and LF groupsA. Bar plot of the first 30 GO terms in HF group; B. Bubbles of the top 20 KEGG pathways in HF group;C. Bar plot of the first 30 GO terms in LF group; D. Bubbles of the top 20 KEGG pathways in LF group

3 讨 论

在家畜育种过程中,因需满足人类的不同生活需求,持续的人工选择在家畜基因组中留下相应的变异信息和选择信号。选择消除能够识别基因组中被消除多态性的重要区域,能够使一些复杂且重要的经济性状相关基因借助选择消除分析而被识别出来。本研究发现,LF组的平均测序深度略高于HF组,但HF组较LF组具有更多的SNP变异位点,表明HF组具有更多的遗传多样性。通过Fst与θπ联合筛选发现,在HF组中SP1、FGFR1、AMHR2、SMAD9、OXTR和ACTG1等基因与雌性的生殖能力、原始卵泡池形成、卵泡发生和颗粒细胞存活等过程紧密相关[18-21];LF组中TGFBR3、BMPR2、SOX11、BCL2、THBS1和ITGA11等基因与卵巢卵泡发育、胚胎发育和妊娠维持等生物过程紧密相关[22-23],且这些基因在关中奶山羊的不同繁殖力组中受到不同程度的选择,表明这些基因可能通过影响母羊的生殖过程,进而参与关中奶山羊产羔数性状的进化选择。

原始卵泡池的建立对于雌性生殖寿命的维持至关重要。在卵巢原始卵泡池的形成过程中,SP1存在于卵母细胞和体细胞中,通过调节哺乳动物卵巢中的颗粒前细胞发育来控制卵巢储备的建立[18]。Cai等[18]通过敲除小鼠前颗粒细胞的SP1表达,发现卵巢破裂、卵母细胞凋亡和原始卵泡池形成等过程被显著抑制,表明体细胞表达的SP1在原始卵泡发生过程中发挥作用。FGF信号在哺乳动物发育的许多方面都发挥着重要作用,FGFR1能够调节滋养外胚层发育并促进胚泡植入,而P65靶向FGFR1能够调节卵巢颗粒细胞的存活,进而促进卵泡生长[19-20]。在人类和小鼠中,FGFR1的失活突变会导致GnRH缺乏和一系列下游生殖障碍[21]。TGFBR3是一种辅助受体,可与TGF-β和抑制素结合并调节其活性,进而调节生殖生物学的多个方面。卵巢卵泡发育受局部产生的TGF-β超家族成员的调节。调节多种TGF-β配体(包括抑制素)作用的TGFBR3在不同的卵巢细胞类型中表达,可在小鼠垂体促性腺细胞中起抑制素A的作用[22]。Li等[23]发现,条件性敲除TGFBR3的雌性小鼠具有高繁殖力,其表现出卵泡发生增强、每个周期的排卵数增加以及产仔数增加现象。细胞间通讯和生长因子信号通路的异常可导致妊娠期间母胎相互作用的缺陷,包括胎儿/胎盘单位的免疫排斥,而胚胎植入后的子宫功能和妊娠维持需要BMPR2[24]。BMPR2是调节母胎界面的几种生理信号通路和事件的中心点,Nagashima等[24]发现,小鼠子宫蜕膜中BMPR2的缺失引发妊娠中期异常,导致血管发育异常、滋养层缺陷和子宫自然杀伤细胞缺乏,并抑制IL-15、VEGF、血管生成素和corin信号传导,这些通路的中断共同导致胎儿死亡和雌性不育。综上所述,本研究筛选到的SP1、FGFR1、TGFBR3和BMPR2等基因可能参与调节关中奶山羊母羊的生殖过程,进而影响其繁殖力。

受选择基因的功能富集能够为目标性状候选基因的筛选提供一定的参考与指导,因此本研究对关中奶山羊不同产羔数群体的受选择基因进行功能富集分析,发现一些已知的生殖相关通路被发现显著富集,包括mTOR信号通路、PI3K-Akt信号通路、Wnt信号通路、TGF-β信号通路和催产素信号通路等。在性腺中,mTOR有助于维持精原干细胞和卵泡的特性,并严格调节两个系统的分化以确保适当的配子产生,而PI3K/AKT/mTOR通路调节涉及卵泡静止、激活、发育、正常卵母细胞成熟和排卵的事件。mTOR作为颗粒细胞增殖的正向调节剂,其失调与卵泡闭锁和卵巢功能障碍有关,且通过卵母细胞-颗粒细胞通讯的mTOR活性对于维持初级卵泡处于静止状态至关重要[25]。在哺乳动物的卵巢中,未成熟的卵母细胞在被激活及排卵前储存在原始卵泡中,而原始卵泡激活的精确控制对于生殖至关重要。Habara等[26]通过wntless基因的条件性敲除,发现抑制前颗粒细胞与颗粒细胞的Wnt分泌导致卵母细胞生长受损和雌性不育,表明前颗粒细胞中的Wnt信号传导是卵巢卵泡发生和雌性生育能力所必需的。TGF-β家族是卵母细胞特异性成员之一,对各种生物的生殖功能有着深远的影响。雌性生殖系统中的TGF-β和激活素通过Smad2和Smad3发出信号,Smad2-/-胚胎不能形成原肠胚,且中胚层形成缺陷和前后轴缺陷,并在胚胎E6.5～E8.5时期死亡[27]。以上研究表明,mTOR信号通路、PI3K-Akt信号通路、Wnt信号通路和TGF-β信号通路等可能通过影响关中奶山羊的生殖过程,进而调控关中奶山羊的产羔数性状。

本研究在关中奶山羊中发现了一些影响产羔数性状的的候选基因,如SP1、FGFR1、AMHR2、SMAD9、OXTR、ACTG1、TGFBR3、BMPR2、SOX11、BCL2、THBS1和ITGA11等,这些基因有可能是应用于关中奶山羊育种实践的候选基因。此外,一些调控繁殖过程的信号通路,如TGF-β信号通路、mTOR信号通路、催产素信号通路、PI3K-Akt信号通路和Wnt信号通路等,可能在关中奶山羊产羔数性状的选择中发挥重要作用。然而,本研究筛选的候选基因是否为调控关中奶山羊产羔数性状的关键基因,仍需在群体中进行进一步验证以保证其准确性,但本研究初步筛选的候选基因可在一定程度上加快关中奶山羊的育种进程,并为实际生产中高繁个体的选育提供便利。

4 结 论

本研究通过全基因组重测序与选择消除分析筛选了与关中奶山羊产羔数性状相关的候选基因,包括SP1、FGFR1、AMHR2、SMAD9、OXTR、ACTG1、TGFBR3、BMPR2、SOX11、BCL2、THBS1和ITGA11等,为揭示关中奶山羊产羔数性状的分子遗传基础提供了相应的理论支撑。