基于轉錄組測序挖掘槐山羊高繁關鍵候選基因

韓浩園,李世凱,2,楊瑞巧,李曼曼,李 君,哈 斯,趙金艷,魏紅芳,權 凱*

(1.河南牧業經濟學院動物科技學院,鄭州 450046;2.華南農業大學動物科學學院,廣州 510642)

槐山羊是我國優良的地方品種,核心產區位于周口市沈丘縣,群體分布覆蓋整個豫東地區,是河南養羊業的第一名片,“槐山羊肉”、“槐山羊板皮(槐皮)”為沈丘縣國家地理標志保護產品。槐山羊以皮質好、肉質鮮嫩、繁殖率高三大優勢而聞名于世[1],槐山羊性成熟早,多胎多產,繁殖性能優異,能一年兩產或兩年三產,經產母羊產羔率平均為329.20%。2003年以來由于受到國際市場的沖擊,波爾山羊等品種引入后對槐山羊大面積雜化,造成槐山羊群體血統含量參差不齊,生產力方向和生產性能也發生了巨大的變化,目前槐山羊多胎性狀的功能基因和遺傳機制尚不明晰。

目前,國內關于羊的高繁基因研究主要集中于綿羊[2-8],針對山羊品種的高繁殖力研究僅限于少數幾個品種,前期基于DNA多態性研究,在隴東絨山羊、黔北麻羊、貴州黑山羊、遼寧絨山羊、陜北白絨山等品種中發現GJB6、DNAH1、NCOA1、FSHR等基因遺傳變異與山羊繁殖力相關[9-14]。隨著測序技術的革新,生命科學研究進入了后基因組時代,誕生了各種組學研究方法,其中轉錄組學研究是從轉錄水平上反映出基因的表達情況及其調控規律,是研究基因功能和篩選新基因的有效技術之一[15-16]。目前,在山羊中已陸續開展與繁殖性能相關的轉錄組研究,Ling等[17]對安徽白山羊卵泡期單、多羔的山羊卵巢進行轉錄組測序,通過比較分析信使RNA(mRNAs)、微RNA(miRNAs)和長非編碼RNA(lncRNAs)的差異表達,發現在多羔組卵巢中miR-6404和miR-29c可能對山羊的繁殖性能有著重要調控作用。此外,大部分研究集中于卵泡不同發育階段的基因表達模式,基于單細胞轉錄組,王俊杰[18]針對濟寧青山羊卵泡不同發育階段的卵母細胞和顆粒細胞,篩選出一部分與階段發育相關聯的功能性轉錄活動。Xu等[19]通過獲取大足黑山羊不同發育階段的完整卵泡轉錄組,發現128個mRNA、4個lncRNAs、49個miRNAs和290個環狀RNA(circRNAs)在大、小濾泡中表達差異顯著,差異mRNA富集到與卵泡發育相關的信號通路中。Liu等[20]通過生物信息學分析確定安徽白山羊卵泡期和黃體期卵巢中差異表達mRNA的表達模式,鑒定出3 770個差異表達mRNA,功能富集分析表明,HSD17B7、3BHSD和SRD5A28基因可能與孕酮的合成有關,RPL12、RPS13和RPL10基因可能與卵母細胞的生長和成熟有關。

目前,對于槐山羊繁殖性能功能基因及其調控機制的研究極少,本試驗以高繁殖力品種槐山羊為研究對象,采樣個體均來自沈丘縣農牧科技研發中心,飼養標準與管理水平一致,依據經產母羊平均每胎產羔數分為單羔組和多羔組,利用轉錄組測序分析篩選與多羔性狀相關的候選基因,可以為探討山羊產羔性狀的分子機理提供理論參考。

1 材料與方法

1.1 樣本采集及建庫測序

槐山羊卵巢組織采自沈丘縣農牧科技研發中心,選取健康無病的3歲左右第3胎次母羊,試驗個體體格大小相似,飼養標準一致,依據研發中心的槐山羊繁殖性能記錄,采集3只平均每胎產羔數為1只的母羊為單羔組(S1～S3),采集3只平均每胎產羔數大于2只(分別為2.3、2.6和3.3只)的母羊為多羔組(M1～M3),屠宰后采集其卵巢組織,立即置于液氮中保存。利用TRIzol法提取卵巢組織總RNA,取RIN值>7.0,完整性良好的RNA送至派森諾進行建庫測序,利用Illumina平臺建庫測序,采用paired-end 2×150 bp測序模式。

1.2 數據處理

檢驗合格樣品經過上機測序,生成 FASTQ 的原始數據(raw data)。采用cutadapt去除3′端帶接頭的序列,去除平均質量分數低于Q20的reads,生成clean data。使用HISAT2(http://ccb.jhu.edu/software/hisat2/index.shtml )軟件BWT算法將過濾后的reads比對到山羊參考基因組(Capra_hircus.ARS1.dna.toplevel.fa)上。

1.3 差異表達基因及功能富集分析

采用DESeq對基因表達進行差異分析,表達差異倍數|log2FoldChange|>1,顯著性P<0.05的基因為差異表達基因(differentially expressed genes,DEGs)。使用R語言pheatmap軟件包對差異基因的并集和樣品進行雙向聚類分析,根據同一基因在不同樣品中的表達水平和同一樣品中不同基因的表達模式進行聚類,采用euclidean方法計算距離,采用層次聚類最長距離法(complete linkage)進行聚類。本研究使用topGO進行GO富集分析,使用KAAS(http://www.genome.jp/tools/kaas/)中bi-directional best hit(BBH)參數進行KEGG功能注釋,P<0.05為顯著富集。

1.4 熒光定量PCR驗證

通過實時熒光定量PCR(qRT-PCR)技術,以GAPDH為內參基因,挑選6個基因(ADCY8、FOS、PAK1、NR4A2、ADAMTS、NR4A1)進行驗證。熒光定量PCR反應體系共10 μL,包括:2×SYBR qPCR master mix 5 μL,上、下游引物(表1)各0.2 μL,cDNA 0.1 μL,ddH2O 4.5 μL;反應程序為:95 ℃預變性3 min;95 ℃ 10 s,60 ℃ 30 s,40個循環。利用2-ΔΔCt方法計算基因表達量,以RNA-Seq方法的基因表達log2(Fold Change)值為橫坐標,以實時熒光定量方法的log2(2-ΔΔCt)值為縱坐標,進行相關分析,并計算皮爾遜相關系數。

1.5 遺傳變異分析

采用varscan程序獲取SNP和InDel位點,過濾標準為:SNP位點堿基Q>20;覆蓋該位點的reads數目>8;支持突變位點的reads數目>2;SNP位點的P<0.01。

2 結 果

2.1 下機及過濾數據

對每個樣品的下機數據(raw data)分別進行統計,并對帶接頭、低質量的reads進行過濾,統計結果見表2。平均reads總數為41 163 239條,平均堿基總數為6 174 485 900 bp,平均clean reads數為37 908 182條,平均clean data堿基數為5 686 227 300 bp,clean data的reads和堿基占raw data的百分比為92.08%。

2.2 參考基因組數據比對

將過濾后的高質量序列與參考基因組序列進行比對,平均95.60%的clean data比對到參考基因組,其中96.21%的序列為特異比對,可用于后續分析,基本信息統計見表3。將比對到基因組上的reads分布情況進行統計,73.24%的reads定位到基因區域,26.76%的reads定位到基因間區,84.56%的reads定位到外顯子區域,比對結果見表4。

表3 卵巢組織轉錄組數據與參考序列比對結果Table 3 The mapping results between transcriptome data from ovarian tissue and reference sequences

表4 轉錄組數據位置比對Table 4 RNA-Seq mapped events

2.3 差異表達基因分析

本研究在多羔與單羔組間共發現121個差異表達基因(|log2(Fold Change)|>1,P<0.05),包括30個上調基因和91個下調基因,差異表達基因火山圖見圖1,藍色為多羔組相比于單羔組下調基因,紅色為多羔組相比于單羔組上調基因。下調基因中顯著性前5的基因為RXFP1、DUSP1、THSD7B、CCN1和MRAP2,上調基因中顯著性前5的基因為ENSCHIG00000021945、SLF1、ENSCHIG00000001740、ND2和BRINP3(表5)。

紅點表示多羔組上調基因,藍點表示多羔組下調基因Red dots represent up-regulated genes, blue dots represent down-regulated genes圖1 差異表達基因火山圖Fig.1 Volcano map of DEGs

表5 多羔組與單羔組間差異表達基因(前10)Table 5 The top 10 DEGs between multi-lamb group and single-lamb group

2.4 聚類分析

基于差異表達基因的表達量,雙向聚類分析結果表明單羔組(Single)S1～S3個體聚為一類,多羔組(Multiple)M1～M3個體聚為一類,由圖2可見,單羔組和多羔組間的差異表達基因的表達模式差異明顯,推測本研究篩選的差異基因表達量可能與槐山羊的產羔數相關。

紅色表示高表達基因,綠色表示低表達基因Red represents up-regulated genes, green represents down-regulated genes圖2 差異表達基因雙向聚類結果Fig.2 Two-way clustering results of DEGs

2.5 功能富集分析

對差異表達基因進行功能富集分析,GO分析共顯著富集到440個功能條目(P<0.05),極顯著富集到100個條目(P<0.01),圖3顯示了顯著性最高的前20個GO條目,其中較多的差異基因富集到發育過程(developmental process)、解剖結構發育(anatomical structure development)、多細胞生物發育(multicellular organism development)和動物器官發育(animal organ development)等功能條目中(表6),PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1和SLC4A10等基因在單羔與多羔組間表達差異顯著,且參與多條組織器官結構發育功能條目,可能影響母羊繁殖系統的發育過程。本研究發現,NR4A1和NR4A2基因富集到促腎上腺皮質激素釋放激素應答(response to corticotropin-releasing hormone)和促腎上腺皮質激素釋放激素刺激的細胞應答(cellular response to corticotropin-releasing hormone stimulus)2個功能條目中(表6),促腎上腺皮質激素釋放激素在外周與生殖功能相關,推測NR4A1和NR4A2基因可能通過參與促腎上腺皮質激素釋放激素應答,進而調控槐山羊的產羔數。

圖3 GO富集分析氣泡圖Fig.3 The bubble diagram of GO enrichment analysis

表6 GO富集分析Table 6 GO enrichment analysis

GO有向無環圖結果表明(圖4),生物過程中的發育過程(developmental process)、多細胞生物發育(multicellular organismal process)和行為(behavior)條目顯著富集,且其下屬功能條目均被顯著富集,其中,發育過程(developmental process)、多細胞生物過程(multicellular organismal process)、解剖結構發育(anatomical structure development)、多細胞生物發育(multicellular organism development)、系統發育(system development)功能條目均包括超過30個差異表達基因,本結果說明槐山羊產羔數高低可能與其細胞生物過程、解剖結構發育及系統發育有關。

圖4 生物過程中部分條目的GO有向無環圖Fig.4 The GO directed acyclic graph based on a part of terms in biology process

KEGG分析將差異表達基因共富集到15條通路中(圖5),其中6條信號通路與內分泌系統有關(表7),包括甲狀旁腺激素的合成、分泌和作用(parathyroid hormone synthesis, secretion and action)、松弛素信號通路(relaxin signaling pathway)、醛固酮合成與分泌(aldosterone synthesis and secretion)、脂肪細胞中脂肪分解的調節(regulation of lipolysis in adipocytes)、皮質醇的合成與分泌(cortisol synthesis and secretion)、催產素信號通路(oxytocin signaling pathway),在多羔組與單羔組比較中發現FOS、MMP13、NR4A1、NR4A2、ADCY8基因差異表達,且富集在多條內分泌相關的KEGG通路中(表7),本結果說明槐山多羔性狀可能與激素分泌及上述5個基因的表達情況有關。

圖5 KEGG富集分析氣泡圖Fig.5 The bubble diagram of KEGG enrichment analysis

表7 KEGG功能富集分析Table 7 KEGG functional enrichment analysis

2.6 實時熒光定量驗證

以GAPDH基因為內參基因,用實時熒光定量法檢測單羔組和多羔組槐山羊卵巢組織中的ADCY8、FOS、PAK1、NR4A2、ADAMTS、NR4A1基因表達量,結果顯示6個基因在槐山羊多羔組中顯著下調,熒光定量PCR結果與RNA-Seq結果一致,分別以RNA-Seq及qRT-PCR的單羔、多羔組間的表達量差異倍數為橫、縱坐標,兩種方法得到的6個基因表達量的皮爾遜相關系數R2為0.971 7,顯著性P<0.01(圖6),證明RNA-Seq的結果準確可靠。

圖6 實時熒光定量PCR驗證結果Fig.6 The validation results of quantitative real-time PCR

2.7 遺傳變異篩選

SNP/Indel篩選結果表明6個個體的SNP數為82 358～106 639個,Indel數為8 856～10 844個(表8)。在上述結果中發現的與槐山羊產羔數相關的11個基因(PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1、SLC4A10、FOS、NR4A1、NR4A2、ADCY8)中共發現236個SNPs,其中39個SNPs位于外顯子區域,7個SNPs為錯義突變,4個位于ADAMTS1基因,1個位于PTX3基因,2個位于CCN1基因(表9)。

表8 SNP/Indel篩選結果Table 8 The screen results of SNP/Indel

表9 與槐山羊產羔數相關基因的非同義突變Table 9 The nonsynonymous mutations related to lambing numbers of Huai goats

3 討 論

繁殖性能包括生育力、生產力和繁殖力3個方面,它們分別指的是母畜每年的妊娠率、每胎的產羔數和每年的產羔數,排卵率和窩產羔數是繁殖性能的直接體現。國內外關于山羊的研究報道相對匱乏,相關重視程度和技術發展遠遠比不上綿羊,與產羔數相關的研究主要集中于黔北麻羊、濟寧青山羊和絨山羊等品種中,對于槐山羊尚無相關報道。基于轉錄組測序技術,本研究構建了6個單、多羔槐山羊卵巢組織mRNA文庫,共發現了121個差異表達基因,其中多羔組上調基因30個,下調基因91個,富集到多條組織器官結構發育功能條目及多條內分泌系統相關激素分泌信號通路。在單、雙羔燕山絨山羊卵巢組織轉錄組測序中發現的差異表達基因富集于繁殖過程相關通路中,包括松弛素信號通路、cAMP 信號通路、TGF-beta信號通路、MAPK信號通路和雌激素信號通路等[21]。在濟寧青山羊的轉錄組分析中也發現類固醇代謝和胰島素通路可能參與了產羔數性狀的進化選擇[18],證明這些信號通路與山羊的產羔數相關。本研究篩選出PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1、SLC4A10、ADCY8、NR4A1和NR4A2等主要候選基因,而燕山絨山羊中篩選的基因包括ADAMTS16、GABRA1、TIMP3、TRPC1、SOS2、WNT2、SLC6A3、PRLHR和SERPINB11等,說明不同品種間控制產羔數的調控基因可能存在差異。在Pelibuey羊的卵巢組織研究中同樣篩選出NR4A1基因可能與羊繁殖相關[22],與本研究結果一致。

本研究在單羔與多羔組間發現的差異表達基因包括PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1、SLC4A10、ADCY8、NR4A1和NR4A2等。已有研究表明其中PAK1、PTX3、MMP13、ADAMTS1、ADCY8、NR4A1和NR4A2基因均在卵巢組織中表達,且與卵泡發育、排卵及類固醇激素合成相關。如在對濟寧青山羊研究中發現PAK1基因與山羊產羔數相關[18],與本研究結果一致。PTX3基因可通過抑制PI3K/AKT/mTOR和Erk1/2信號通路,從而調控奶山羊卵巢顆粒細胞的增殖、凋亡,并調控類固醇合成關鍵酶CYP19A1和CYP11A1的表達,進而作用于顆粒細胞類固醇激素的生成[23]。MMP13膠原酶則參與了卵泡壁的降解,在卵泡形成和排卵中具有重要功能[24-26]。ADAMTS1基因mRNA和蛋白在優勢卵泡周圍富集,參與子宮和卵巢的形態發生、子宮內膜蛻膜化、卵子的發生和排卵過程,促進卵泡的發育和成熟[27-28]。也有研究發現,ADCY8基因參與卵巢類固醇生成、G蛋白和雌激素信號通路,與卵母細胞減數分裂、孕激素介導的卵母細胞成熟有關[29]。

此外,本研究發現NR4A1和NR4A2基因參與促腎上腺皮質激素釋放激素應答,促腎上腺皮質激素釋放激素在外周與生殖功能相關,推測NR4A1和NR4A2基因可能通過促腎上腺皮質激素釋放激素應答,進而調控槐山羊的產羔數。前人研究發現NR4A1基因作為一個孤核受體的轉錄因子,編碼類固醇-甲狀腺激素超家族,同樣調節關鍵類固醇激素合成酶基因的轉錄[30-31],而孤兒核受體NR4A2與細胞生長或凋亡密切相關,可能參與卵巢上皮性癌的發生與發展[32]。結合本研究結果,推測PAK1、PTX3、MMP13、ADAMTS1、ADCY8、NR4A1和NR4A2基因表達于山羊卵巢組織,且調控山羊卵泡發育及產羔數。

本研究在與槐山羊產羔數相關的11個基因(PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1、SLC4A10、FOS、NR4A1、NR4A2、ADCY8)中共發現236個SNPs,其中編碼蛋白發生改變的SNPs為7個,位于ADAMTS1、PTX3和CCN1基因中,本研究所發現的遺傳變異位點與產羔數的關系需要進一步驗證。目前尚無這些基因的SNP與山羊產羔數關系的研究。已發現的與中國本地山羊高繁性能有關的SNPs研究主要集中在FSHR、NCOA1、RBP4和DNAH1等基因[9]。在長白豬群體的NR4A1基因第5內含子發現了A/G突變,多態性分型結果表明胎次產仔數呈現GG>AG>AA的趨勢[33],該基因SNP在山羊產羔數中的作用有待進一步研究及驗證。

4 結 論

本研究構建了多羔及單羔共6個個體的mRNA文庫,發現PTX3、MMP13、PAK1、ADAMTS1、COL1A2、CCN1、SLC4A10、FOS、NR4A1、NR4A2、ADCY8等11個基因顯著差異表達,且參與到組織器官的結構發育及內分泌系統相關激素分泌,可能通過調控山羊卵泡發育、排卵及類固醇激素分泌,進而影響山羊產羔數。本研究為解析山羊多羔性狀遺傳機制提供了理論依據。