福清山羊與努比亞黑山羊發(fā)情期卵巢組織RNA-Seq分析

李文楊，劉遠，吳賢鋒，高承芳，黃勤樓

（福建省農業(yè)科學院畜牧獸醫(yī)研究所，福州 350013）

0 引言

【研究意義】繁殖性能是養(yǎng)羊業(yè)的重要經(jīng)濟性狀，能繁母羊產仔數(shù)的多少直接影響?zhàn)B羊業(yè)的經(jīng)濟效益。由于羊繁殖性狀受到遺傳力低、限性性狀等因素制約，常規(guī)選育存在周期長、遺傳進展慢的缺點，最經(jīng)濟有效的方法是通過選擇控制多胎性能的主效基因（分子標記）提高母羊的產羔率[1-2]。因此研究羊繁殖性狀的分子調控機理，挖掘、鑒定和充分利用控制羊繁殖性狀的主效基因，并應用于羊繁殖性狀的遺傳改良具有重要意義[3]。【前人研究進展】轉錄組是特定組織或細胞在某一環(huán)境或生理條件下所轉錄表達的所有 RNA總和，是連接基因組遺傳信息與蛋白質組生物功能的紐帶[4]。高通量轉錄組測序技術（RNA sequencing，RNA-Seq）的發(fā)展和成熟為系統(tǒng)研究基因表達及調控提供了重要的手段和方法，在畜禽重要經(jīng)濟性狀分子機制研究中得到了廣泛使用[5-7]?；诟咄繙y序技術，研究者采用不同的策略開展了羊繁殖性狀候選基因的篩選與鑒定。羅慶苗等[8]利用數(shù)字基因表達譜技術（digital gene expression tag profiling，DGE）分別檢測了處于自然發(fā)情期的 FecBBFecBB型與FecB+FecB+型小尾寒羊卵巢組織的基因表達情況，認為類固醇快速調節(jié)蛋白（steroidogenic acute regulatory protein，STAR）、卵巢濾泡激素（Folliculin，F(xiàn)LCN）等5個基因與Notch信號通路可能是參與調控小尾寒羊高低繁殖性狀差異的重要基因和調控通路。DI等[9]利用 RNA-Seq技術獲得了綿羊不同繁殖狀態(tài)下的基因表達模式。LING等[3]比較了發(fā)情期多羔和單羔安徽白山羊卵巢組織的差異表達情況，篩選出了12個可能與山羊高繁殖力有關的基因。FENG等[10]研究了高產和低產湖羊卵泡期卵巢組織 RNA表達差異，篩選出了 76個差異表達基因和 5個差異表達長鏈非編碼RNA（Long non-coding RNA, lncRNA）。付紹印等[11]對巴美肉羊性腺軸進行了轉錄組分析，功能富集分析表明，神經(jīng)活性的配體-受體相互作用信號通路在下丘腦-垂體-卵巢性腺軸調控排卵及卵泡發(fā)育方面發(fā)揮著重要調控作用?！颈狙芯壳腥朦c】目前對部分綿羊、山羊品種的繁殖力相關基因研究較多，但對福清山羊的相關研究較少。而影響不同品種繁殖力的基因也不盡相同，李文楊等[12]證實控制部分綿羊、山羊多胎性狀主效基因骨形態(tài)發(fā)生蛋白 15（bone morphogenetic protein 15，BMP15）、骨形態(tài)發(fā)生蛋白受體1B（bone morphogenetic protein receptor 1B，BMPR-1B）和生長分化因子9（growth differentiation factor 9,GDF9）并非控制福清山羊繁殖性狀的主效基因。努比亞黑山羊母羊繁殖指數(shù)（dam reproduction index，DRI）高達0.99，而福清山羊母羊DRI僅為0.45[13]。同時有研究表明利用努比亞黑山羊改良貴州黑山羊可提高后者產羔率37.2%[14]。因此，導入努比亞黑山羊基因改良福清山羊品種，對其繁殖、生產性能的提高具有重要的生產實踐意義?！緮M解決的關鍵問題】本研究擬以福建地方山羊品種——福清山羊和福建省引進的努比亞黑山羊為材料，采用RNA-Seq技術比較發(fā)情期卵巢組織的差異表達基因（differentially expressed genes,DEGs）及新轉錄本分析，并對篩選的DEGs進行GO（gene ontology）、KEGG（kyoto encyclopediaof genes and genomes）功能富集和COG（Cluster of Orthologous Groups of proteins）分類，挖掘控制福清山羊繁殖性能的相關候選基因，一方面為地方品種種質資源的遺傳機制研究提供理論依據(jù)，另一方面通過測序數(shù)據(jù)與參考基因組的比對，篩選潛在可能影響福清山羊繁殖性能的單核苷酸多態(tài)性（single nucleotide polymorphisms，SNP）和插入缺失標記（insertion-deletion，InDel）分子遺傳標記，為福清山羊繁殖性能分子遺傳改良提供新的育種素材。

1 材料與方法

1.1 試驗材料

本研究選用的福清山羊和努比亞黑山羊由福建省農業(yè)科學院福清市漁溪肉羊試驗基地提供，試驗于2017年2—6月進行。選擇2—3歲營養(yǎng)良好、體態(tài)勻稱的空懷雌性福清山羊和努比亞黑山羊各5只。母羊發(fā)情鑒定以公羊試情為主，結合母羊外陰部觀察。試情公羊每日早晚各試情 1次（早 9：00、晚 21：00），母羊連續(xù)2次試情均接受試情公羊的爬跨，并在爬跨時靜立不動即鑒定為發(fā)情。鑒定母羊發(fā)情后按照國家實驗動物處理飼喂準則屠宰，取有成熟卵泡發(fā)育一側的卵巢組織，置于-80℃?zhèn)溆谩Ｃ總€品種屠宰3只，福清山羊標記為F-1、F-2、F-3，努比亞黑山羊分別標記為 N-4、N-5、N-6。

1.2 RNA提取、轉錄組文庫構建及測序

利用TaKaRa MiniBEST Universal RNA Extraction Kit （TaKaRa，大連寶生物）提取每個樣品卵巢組織的總 RNA，單個建池。為了保證數(shù)據(jù)準確性，采用Nanodrop檢測RNA的純度，Agilent 2100精確檢測RNA的完整性。RNA樣品檢測合格后，送北京百邁克生物科技有限公司進行文庫構建以及測序。文庫構建完成后，對建成文庫的質量進行檢測。文庫質量合格后，不同的文庫利用 Illumina HiSeq2500平臺（Illumina，America），基于邊合成邊測序（sequencing by synthesis，SBS）技術進行高通量測序。

1.3 轉錄組測序數(shù)據(jù)處理與分析

原始測序數(shù)據(jù)經(jīng)過去除含有接頭和低質量數(shù)據(jù)過濾后獲得高質量測序數(shù)據(jù)。利用 TopHat2[15]軟件將獲得的 Clean reads比對到山羊參考基因組（ftp：//ftp.ncbi.nlm.nih.gov/genomes/all/ GCF/001/704/415/GCF_001704415.1_ASM170441v1/GCF_001704415.1_ASM 170441v1_genomic.fna.gz）[16]，獲取序列在參考基因組或基因上的位置信息，以及樣品特有的序列特征信息。

1.4 差異表達基因的篩選、功能注釋和富集分析

參考JIANG等報道的數(shù)學模型[17]，使用Cufflinks軟件對轉錄本和基因的表達水平進行定量。采用FPKM作為衡量轉錄本或基因表達水平的指標[18]。皮爾遜相關系數(shù)r作為生物學重復相關性的評估指標[19]。

以福清山羊為對照，利用DEseq進行2個品種間的差異表達分析[20]。將差異倍數(shù)≥2且錯誤發(fā)現(xiàn)率（False Discovery Rate，F(xiàn)DR）＜0.01作為DEGs的篩選標準。繪制2個品種DEGs火山圖并進行聚類分析。分別利用GO數(shù)據(jù)庫[21]、COG數(shù)據(jù)庫[22]和KEGG數(shù)據(jù)庫[23]對DEGs進行生物信息學分析。

1.5 SNP/InDel分析

基于各樣品 reads與參考基因組序列的 TopHat2比對結果，使用 GATK軟件[24]識別測序樣品與參考基因組建的單堿基錯配和插入缺失，識別潛在的SNP及 InDel位點。進一步利用 SnpEff軟件[25]注釋篩選到的變異（SNP、InDel）和預測變異造成的影響。

1.6 熒光定量PCR驗證

為了進一步驗證高通量測序結果以及發(fā)現(xiàn)新轉錄本序列的準確性，以18S rRNA為內參基因，隨機挑選6個 DEGs和2個新轉錄本，采用實時熒光定量PCR（qRT-PCR）驗證基因表達量。內參基因引物參考文獻設計[21]，其余基因和新轉錄本引物根據(jù)測序序列信息由ABI公司的Primer Express Software v2.0設計，鉑尚生物技術（上海）有限公司合成，基因名稱及引物信息見表1。以1.2中各樣品RNA池為模板，采用TAKARA PrimeScript? RT reagent Kit with gDNA Eraser（TaKaRa，大連寶生物）試劑盒逆轉錄后進行qRT-PCR擴增，擴增體系20μL：cDNA模板1μL，上、下游引物各 1μL，2×SYBR? Select Master Mix 10μL，ddH2O7μL。qRT-PCR擴增程序：95℃預變性2 min；95℃變性10s，60℃退火10s，72℃延伸40s，50個循環(huán)；最后 72℃延伸 5min；4℃保存。采用 2-ΔΔCt法計算基因的相對表達量。

表1 實時熒光定量PCR引物信息Table 1 The qRT-PCR primers

2 結果

2.1 測序數(shù)據(jù)評估及對比分析

經(jīng)過文庫質量檢測和測序質量控制，各樣品均符合測序要求。6個樣品共得到46.68Gb Clean Data，各樣品Clean Data 均達到7.24Gb以上，獲得Clean reads 24 283 922—27 945 197條。測序數(shù)據(jù)堿基質量（Quality Score）Q30堿基百分比為91.15%—91.99%，堿基GC含量為50.98%—52.19%，各堿基質量穩(wěn)定在20%—30%之間，低質量堿基比例小，測序質量好。表明文庫構建質量和測序質量高，測序數(shù)據(jù)準確可靠，滿足后續(xù)分析條件。

將獲得的 Clean reads與參考山羊基因組進行對比分析，結果表明各樣品的 Clean reads與參考基因組的比對效率在 84.62%—87.14%之間，比對到參考基因組唯一位置的Reads占Clean Reads的百分比為 81.28%—84.35%，多處位置的百分比為2.30%—3.71%（表2），表明6個樣本總RNA不存在污染。

表2 測序數(shù)據(jù)與參考基因組的序列比對結果統(tǒng)計Table 2 RNA-Seq and mapping to the reference genome

2.2 DEGs的篩選、功能注釋和富集分析

2.2.1 DEGs的篩選 通過定位到參考基因組區(qū)域或基因外顯子區(qū)的Reads的計數(shù)來估計基因（轉錄本）的表達水平，基因（轉錄本）的真實表達水平與Reads計數(shù)、基因（轉錄本）的長度和測序深度成正相關。采用FPKM作為衡量轉錄本或基因表達水平的指標，計算得到了6個樣本18 996個表達基因的FPKM值。

為了保證測序數(shù)據(jù)的可重復性以及可靠性，分析了6個樣品之間轉錄本表達水平的相關性。皮爾遜相關系數(shù)的平方（r2）越接近于1，說明2個樣品的表達模式相似度越高。Encode計劃建議r2大于0.92，福清山羊3個樣品的r2為0.926—0.935，努比亞黑山羊3個樣品的r2為0.928—0.947，表明2個品種內部樣品之間表達模式的相似度較高，重復性好，試驗結果可信度高。

根據(jù)差異基因篩選標準，以福清山羊為對照，獲得了努比亞黑山羊和福清山羊發(fā)情期卵巢組織的DEGs 149個，其中上調基因53個，下調基因96個，圖2為DEGs火山圖。進一步根據(jù)DEGs在2個品種卵巢組織的差異表達值和文獻檢索，初步認為輸卵管素基因（oviductin，OVN）、類固醇合成快速調節(jié)蛋白基因(steroidogenic acute regulatory protein，STAR)、早期生長應答1基因(early growth response 1，EGR1)可作為福清山羊繁殖性能的候選基因。

圖1 樣品的表達量相關性熱圖Fig. 1 Correlation heat map of gene expression level in 6 samples

2.2.2 DEGs的GO功能富集 GO注釋系統(tǒng)是一個有向無環(huán)圖，包含生物學過程（biological process,BP），細胞組分（cellular component, CC）和分子功能（molecular function, MF）3個主要分支。篩選出的149個DEGs中的108個基因能夠被GO數(shù)據(jù)庫注釋（圖3）。被注釋的DEGs分別參與了BP、CC和MF的22、19和20個功能亞分類。在BP分類中，細胞過程、單一的生物過程和生物調節(jié)參與的DEGs最多，與繁殖相關聯(lián)的繁殖過程和繁殖的DEGs也被注釋。在CC分類中，參與細胞部分；2細胞；3細胞器的DEGs最多。在MF分類中，DEGs在結合中所占比例最高，其次為催化活性及分子傳感器活性。

圖2 差異表達基因火山圖Fig. 2 Volcano plot of differentially expressed genes

圖3 差異表達基因GO注釋Fig. 3 GO annotation of DEGs

2.2.3 差異表達基因的COG分類 利用COG數(shù)據(jù)庫可以對基因產物進行直系同源分類。篩選出的149個DEGs中，有30個DEGs能夠被COG數(shù)據(jù)庫注釋，分類統(tǒng)計結果見圖 4。其中一般的功能預測基因數(shù)目最多，其次為脂質運輸與代謝。

2.2.4 差異表達基因的KEGG注釋 經(jīng)過分析，可被注釋到KEGG數(shù)據(jù)庫的91個DEGs參與了多個代謝通路，將差異表達基因KEGG的注釋結果按照KEGG中通路類型進行分類，分類圖見圖 5。差異表達基因參與最多的是人類疾病相關的通路，其次是生物系統(tǒng)相關通路和環(huán)境信息處理相關通路。

以KEGG數(shù)據(jù)庫總通路為單位，應用超幾何檢驗，找出與整個基因組背景相比，在DEGs中顯著性富集的通路，KEGG通路分析表明DEGs共富集到21條信號通路中，6條通路顯著富集（P＜0.05）。圖 6為DEGs的KEGG通路富集分析結果，圖中呈現(xiàn)了顯著性Q值最小的20個通路，其中自身免疫性甲狀腺病、同種異體移植排斥、吞噬體、病毒性心肌炎、哮喘 5個代謝通路的富集顯著性最可靠，推測其相關的DGEs可能與山羊繁殖性能相關。

圖4 差異表達基因COG注釋分類統(tǒng)計圖Fig. 4 COG annotation classification statistics of DEGs

2.3 新基因分析

使用Cufflinks軟件對Mapped Reads進行拼接，并與參考基因組注釋信息進行比較，尋找原來未被注釋的轉錄區(qū)，發(fā)掘該物種的新轉錄本和新基因，補充和完善參考基因組注釋信息。過濾掉編碼的肽鏈過短（少于 50個氨基酸殘基）或只包含單個外顯子的序列，共發(fā)掘1 506個新基因（轉錄本），其中25個在2個品種中差異表達。

2.4 SNP/InDel分析

經(jīng)過與參考基因組序列的對比，已知的 124個DEGs中發(fā)現(xiàn)了3 689個SNP/InDel位點，其中6個DEGs未發(fā)現(xiàn)突變（圖7）。位于內含子區(qū)域（INTRON）的 SNP/InDel位點最多，占基因序列突變總量的43.2%，其次為下游序列（DOWNSTREAM），占比16.8%。而引起基因編碼序列非同義突變（NONSYNONYMOUS-CODING）的 SNP/InDel位點占比7.8%。

圖5 差異表達基因顯著富集的KEGG通路Fig. 5 List of KEGG pathway for DEGs

2.5 qRT-PCR驗證測序

6個基因（轉錄本）的 qRT-PCR結果見圖8，選擇的目的基因在2個品種個體間表達變化模式與轉錄組測序結果一致，表明本研究利用轉錄組測序獲得差異表達基因以及新基因（轉錄本）序列數(shù)據(jù)較為準確。

圖6 差異表達基因KEGG通路富集散點圖Fig. 6 Enriched scatter map of DEGs KEGG pathway

圖8 測序結果的qRT-PCR驗證Fig. 8 Validation of sequencing results by qRT-PCR

3 討論

隨著高通量測序技術的發(fā)展和廣泛應用，動物遺傳學從以前候選基因、單一性狀為主的微觀研究轉向以物種進化構架內的全基因組水平、多性狀、組學化的宏觀研究，利用高通量測序技術研究畜禽主要經(jīng)濟性狀取得了較大的研究進展[26]。卵巢是雌性動物重要的繁殖器官，具有提供卵細胞和產生類固醇激素的功能，常用于發(fā)情和排卵等繁殖生理過程的研究。羅慶苗等[8]利用 DGE技術研究了小尾寒羊高繁殖性狀相關基因，但篩選的DEGs直接與已知動物繁殖性狀相關基因極少，DI等[9-10]在不同品種研究中也得到了類似結果，與本研究結果一致，并且不同品種篩選出的與繁殖性狀相關的DEGs也不盡相同。蘭道亮等[27]認為部分DEGs可能涉及卵巢的相關功能及活動，但缺乏現(xiàn)有試驗驗證；同時可能也與器官功能的復雜性和多樣性有關[8，28]。

Ovn具有促進受精和早期胚胎發(fā)育的作用[29]，張譯夫等[30]研究表明其還在卵母細胞成熟和排卵過程中具有重要的生物學作用。Ovn的合成受卵巢類固醇激素的調控[31]，STAR編碼的蛋白通過增強膽固醇轉化為孕烯醇酮在類固醇激素合成的急性調節(jié)中起著關鍵作用[32]，羅慶苗等[8]研究認為STAR的顯著下調可導致 FecBBFecBB小尾寒羊性激素合成能力下降而提高繁殖力。本研究發(fā)現(xiàn)，與福清山羊對比，努比亞黑山羊發(fā)情期卵巢組織的Ovn表達量顯著上調，而STAR表達量顯著下調，很可能是造成2品種繁殖性能差異的關鍵因素之一。EGR1涉及不同組織多種基因的反式激活及功能調控，在牛的排卵前卵泡中EGR1具有促性腺激素依賴性的誘導表達[33]，吳陽升等[34]研究表明EGR1參與綿羊卵泡早期發(fā)育功能，本研究中該基因顯著上調表達，可作為福清山羊繁殖性能候選基因進行進一步驗證研究。此外，其他 DEGs（包括新轉錄本）與繁殖性能相關的確切功能及作用尚不明確，有待進一步的實驗驗證。

由于針對不同羊品種 RNA-Seq篩選到的 DEGs存在差異，對應差異表達基因KEGG顯著富集到的代謝通路也不盡相同[9-11]。羅慶苗等[8]在小尾寒羊中富集到細胞內吞過程、Notch信號通路和嘧啶代謝3個調控通路，LING等[3]在安徽白山羊繁殖力差異研究中富集到細胞外基質受體相互作用、粘著斑和肌動蛋白細胞骨架調節(jié)等11條信號通路，本研究富集到自身免疫性甲狀腺病、同種異體移植排斥、吞噬體、病毒性心肌炎、哮喘等免疫和人類疾病的代謝通路，與蘭道亮等[27,35-36]的研究結果相似，MIAO 等[35]認為疾病相關通路的富集表明代謝、免疫等途徑可能調控家畜繁殖力，蘭道亮等[27]則認為DEGs富集的疾病通路雖然源于疾病發(fā)生的生理過程，但這些通路中的部分基因可能涉及卵巢的相關功能及活動。表明畜禽繁殖性能形成的分子機理其調控機制十分復雜，造成不同品種間高繁殖力差異的調控基因可能存在差異。

繁殖性狀是羊的重要經(jīng)濟性狀，而繁殖性狀是受多基因控制的，遺傳力低，傳統(tǒng)的選擇方法對其遺傳進展的改良非常緩慢，利用顯著影響羊繁殖性能的分子標記進行輔助選擇能提高選育效果[29]。但研究表明分子標記具有品種特異性，因此在不同品種中篩選出更豐富、更適用的分子標記是當前的重要工作之一 。RNA-Seq高通量轉錄組技術在完善基因結構及挖掘新轉錄本及新基因方面具有強大的優(yōu)勢[27]。本研究發(fā)掘了1 506個新基因（轉錄本），其中25個為2個品種的差異表達基因。經(jīng)過進一步與參考基因組序列的對比，在已知的 124個 DEGs中發(fā)現(xiàn)了 3689個SNP/InDel位點，引起基因編碼序列非同義突變（NON-SYNONYMOUS-CODING）的 SNP/InDel位點占比7.8%，驗證這些潛在影響福清山羊繁殖性能的分子標記、評價其在福清山羊群體中的遺傳效應，篩選出控制福清山羊繁殖性狀的主效基因（分子標記）是下一步研究的關鍵。

4 結論

本研究通過對DRI指數(shù)存在明顯差異的福清山羊和努比亞黑山羊的發(fā)情期卵巢組織進行 RNA-Seq分析，在轉錄組水平上篩選出了可能與羊繁殖性能的149個DGEs，并發(fā)掘了1 506個新基因（轉錄本）。初步認為OVN、STAR、EGR1在福清山羊和努比亞黑山羊發(fā)情期卵巢組織中的差異表達可能是造成2品種繁殖力差異的關鍵因素之一，可作為福清山羊繁殖性能的候選基因進行深入驗證研究。該研究為揭示福清山羊繁殖性能形成的分子機制提供了線索，為開發(fā)和利用畜禽優(yōu)良性狀分子標記提供了參考依據(jù)。