不同生理階段邊雞卵巢轉錄譜的構建及卵泡發育相關基因的分析

張俊珍，李彩娥，劉 博，李建慧，張 蒙，李亞莉，常強強

(山西農業大學動物科學學院，太谷 030801)

家禽的就巢性影響著其產蛋性能，就巢性越強的品種產蛋性能越低。卵泡的發育與環境、營養和激素調控有關系，決定著家禽就巢行為的發生和維持，對產蛋性能具有重要影響。卵泡發育取決于細胞的增殖、分化和卵黃質大量沉積，進而促進了卵泡質量的快速增長[1]。在家禽生產過程中，一些與卵泡發育相關的基因參與調控家禽的就巢性。隨著家禽選育程度的不斷提高，多基因調控研究在家禽遺傳改良中得到廣泛的應用。

邊雞是肉蛋兼用的中國地方品種，主要分布在內蒙古自治區與山西省北部相毗連的長城內外一帶，具有蛋重大、肉質好、適應性強、耐粗飼和抗寒等優點。由于邊雞具有很強的就巢性，產蛋量很少，限制了該品種資源的開發和利用。因此，提高邊雞產蛋性能已受到相關家禽育種研究者的廣泛關注。家禽育種研究者通常以卵巢作為研究家禽就巢性的試驗材料，從分子水平上研究與就巢性相關基因的差異表達。江紅霞等[2]對克氏原螯蝦的卵巢和肝胰腺進行轉錄組測序，共獲得53 006個單基因(unigene)，卵巢與肝、胰對比分析，發現卵巢中有20 382個差異性表達基因(differentially expressed genes, DEGs)，通過GO功能分類分析發現，部分DEGs被注釋到繁殖(reproduction)、繁殖過程(reproduction process)、免疫系統過程(immune system process)和生長(growth)GO條目。孟金柱等[3]對貴州白山羊的優勢卵泡和從屬卵泡進行轉錄組測序，篩選出233個上調表達基因，并對這些基因進行GO、KEGG分析，最終篩選出4個與卵泡發育有關的基因。郝慶玲等[4]對山西省文水縣海福特母牛的不同生理階段發育卵泡進行轉錄組分析，共獲得608個差異表達基因，其中有9個參與牛卵泡發育的信號通路調控，并篩選出11個直接參與信號通路調節及細胞增殖分化的基因。朱志明等[5]對山麻鴨開產期和產蛋高峰期卵巢組織轉錄組測序，發現存在1 929 個 差異表達基因，其中上調基因 989 個，下調基因 940 個，并將這些基因在Nr、GO、KEGG數據庫進行分析比較，獲得了差異表達基因的功能、分類和代謝通路。Shen等[6]使用轉錄組測序來分析靖黃雞三種類型的卵泡細胞中circRNA和mRNA，利用生物信息學分析確定了14 502個circRNA，其中5 622個廣泛分布在卵泡細胞發育的所有階段，差異表達分析表明，在與生殖相關的途徑中富集了不同表達的circRNA和mRNA，包括TGF-β信號通路、卵母細胞減數分裂和血管平滑肌收縮，為卵泡發育提供了理論依據。董新龍等[7]利用轉錄組對京海黃雞卵巢組織進行分析，發現了4 431個新轉錄本，這些新轉錄本比對到 GO、Nr、KEGG數據庫的數量分別為 917、1 795、422。與這些數據庫的成功比對，為進一步證實新基因的存在提供了基礎；此外，經過與KEGG數據庫比對，發現絲裂原活化蛋白激酶信號通路(MAPK signaling pathway)排在富集通路的首位。雖然轉錄組測序已經在畜禽不同生理階段卵巢組織差異表達基因的研究中廣泛應用，但是目前尚無對于邊雞卵巢組織轉錄組差異表達基因的研究。本研究為了挖掘邊雞卵巢卵泡發育過程中起調節作用的遺傳因子，揭示卵巢發育期、開產期、產蛋高峰期和就巢期卵巢組織的差異表達基因，篩選與卵泡發育相關的候選基因，利用轉錄組測序技術構建不同日齡的邊雞卵巢的轉錄組文庫，對邊雞在不同生理狀態下卵巢組織做轉錄組差異分析、差異表達基因GO功能富集分析和KEGG分析，將差異表達基因進行功能注釋，旨在研究與邊雞卵巢發育相關的基因，以揭示邊雞就巢性強的遺傳機制。

1 材料與方法

1.1 試驗動物和卵巢樣本采集

試驗用邊雞在標準化飼養條件下進行飼養管理，每次取材前進行連續7 d的產蛋觀察：產下第一枚蛋起，表示性成熟；產蛋率達到70%以上，代表產蛋高峰；320日齡后產蛋率下降為30%以下，觀察記錄一些健康的停產并表現出就巢行為的個體。隨機選取70日齡(卵巢發育期)、165日齡(開產期)、220日齡(產蛋高峰期)和330日齡(就巢期)健康、體重相近且符合取樣條件的邊雞個體，每個日齡5只，頸靜脈放血法處死，采集卵巢組織，迅速投入液氮中，轉運后置于-80 ℃超低溫冰箱中保存備用。

1.2 cDNA文庫構建和測序

采用TRIzol試劑(美國Invitrogen)提取不同日齡邊雞的卵巢組織總RNA，用Agilent 2100生物分析儀進行總RNA品質檢測。用帶有poly-T oligo附著的磁珠富集含有poly A的mRNA，使用逆轉錄酶和隨機引物將mRNA反轉錄合成第一條cDNA鏈，使用DNA polymerase I和RNase H合成第二條cDNA鏈。通過PCR擴增方法純化和富集cDNA，構建邊雞卵巢cDNA文庫。在BGISEQ-500平臺進行測序(華大基因公司，深圳，中國)。

1.3 生物信息學分析

過濾低質量片段，去除帶接頭污染和未知堿基含量比例過高的片段獲得clean reads。Clean reads與參考基因組比對預測新基因和轉錄本，檢測基因剪接事件，識別卵巢樣本間的差異表達基因(DEGs)，進行聚類分析和功能注釋。

1.4 新轉錄本預測

利用HISAT軟件將clean reads定位到雞參考基因組。基因組定位后，用StringTie[8]重建轉錄本，基于基因組注釋信息，用Cuffcompare[9]識別新的轉錄本，并使用CPC[10]預測新轉錄本的編碼能力。

1.5 基因表達、DEGs檢測、GO功能富集和KEGG分析

利用Bowtie2[11]將新的編碼轉錄本與轉錄本合并，將clean reads映射成完整的參考基因組后，利用RSEM[12]計算基因表達水平。使用Cor計算卵巢樣本間的Pearson相關性，使用R軟件里的 hclust 函數進行層次聚類分析，圖形的繪制采用了R軟件中的ggplot2包。

采用PossionDis[13]算法檢測不同日齡邊雞卵巢間的差異表達基因。根據GO注釋結果，對DEGs進行功能分類和生物通路分類，同時使用R軟件中的phyper函數進行GO的功能和KEGG富集，用超幾何分布函數計算出P值，計算每個P值的錯誤發現率(FDR)，當FDR≤0.01時則被認為顯著富集。

1.6 實時熒光定量PCR(qPCR)驗證

為了驗證轉錄組測序的結果，從已知參與卵泡發育的差異表達基因中隨機選擇11個基因，通過qPCR進行驗證。將上述從測序同樣的個體樣本中提取的不同日齡邊雞卵巢的總RNA用MMLV逆轉錄酶(TaKaRa，大連，中國)將其逆轉錄為cDNA。通過SYBR方法進行qPCR驗證(引物見表1)。擴增條件為：預變性95 ℃ 30 s；(變性95 ℃ 5 s、退火55 ℃ 30 s、延伸72 ℃ 10 s)×40循環；熔解過程(95 ℃ 1 min、55 ℃ 30 s、55 ℃ 30 s)。內參基因β-actin與目的基因的樣本要在同樣的條件下進行qPCR，通過2-△△CT法計算目的基因的相對表達量，結果以“平均值±標準差(Means±SD)”表示，用Graph-Pad Prism 8.0軟件分析數據，采用one-way ANOVA進行單因素方差分析。每個樣本做4次技術重復。

表1 引物序列

2 結 果

2.1 轉錄組測序原始數據及轉錄本分析

通過在BGISEQ-500平臺上對不同日齡邊雞卵巢樣本進行測序，構建基因表達譜，每個樣品約產生11.04 Gb的序列。共鑒定出18 891個基因，其中已知基因18 063個，預測的新基因828個。此外，共檢測出21 299個新的轉錄本，其中14 376個屬于已知基因且具有蛋白編碼潛力的新轉錄本，925個屬于新基因的具有蛋白編碼潛力的新轉錄本，其余5 998個為長鏈非編碼RNA。

2.2 不同生理階段邊雞卵巢基因表達及相關性分析

通過RSEM分析每組樣本的基因表達水平，計算每兩個樣本之間的Pearson相關性并對所有樣本進行層次聚類分析。結果表明，在165與220 d邊雞卵巢基因表達聚類最近，二者間的Pearson值為0.983，說明有顯著相關性；而70和330 d與165和220 d邊雞卵巢基因表達均未聚類，且相互間相關性相對較低(圖1)。

樣品間彼此靠的越近代表樣品的表達譜越相似The closer the samples are to each other, the more similar the expression profiles of the samples are圖1 不同日齡邊雞卵巢表達基因的層次聚類和Pearson相關性Fig.1 Hierarchical clustering and Pearson correlation of genes expressed among the ovaries of Bian chicken at different ages

2.3 不同日齡邊雞卵巢基因表達譜中的差異表達基因

通過PossionDis計算差異表達基因水平，70、165、220和330 d 邊雞4個發育階段的卵巢轉錄組數據兩兩對比，形成6組差異表達的mRNA(圖2)。通過比較發現，70 d邊雞卵巢轉錄組與其他階段差異最大，獲得的差異表達基因也相對較多，而且下調基因的數量大于上調基因；而165與330 d 相比、220與 330 d相比，差異表達基因的上調數量大于下調基因的數量；165與220 d相比，差異表達基因的總量最少。

以差異倍數>1.5，FDR<0.01為標準，使用DESeq2R包對轉錄組數據進行篩選，找出了在不同生理階段與發育卵泡數變化相關的差異表達基因40個(表2)。結果表明，有18個基因參與卵泡發育的整個過程，有22個基因僅在220 d前表達。其中35個基因的表達量與發育卵泡數成正相關，FDX1L在70、220、165和330 d邊雞卵巢中表達量依次遞減，NDRG4在70、165、220和330 d邊雞卵巢中表達量隨日齡增加呈遞減趨勢，MC3R和PRLHR在70和330 d邊雞卵巢中表達量高于165和220 d的表達量，VIP僅在70、165和220 d邊雞卵巢中表達。

圖2 不同日齡邊雞卵巢中的差異表達基因Fig.2 Differentially expressed genes in the ovary of Bian chickens at different ages

表2 不同日齡邊雞卵巢中可能與卵泡發育相關的差異表達基因及其表達量(RPKM值)

(轉下頁 Carried forward)

(轉下頁 Carried forward)

2.4 不同日齡邊雞卵巢差異表達基因GO和KEGG富集分析

GO分為三大類，即生物過程、細胞組分、分子功能。將4個混合樣本相互組合形成6個對比組，對篩選獲得的差異表達基因進行GO功能富集分析并繪制差異表達基因GO富集柱狀圖(圖3)。結果表明：70與165 d相比，共有1 748個差異基因，其中有615個 得到功能注釋；70與220 d相比，有1 704個差異基因，其中596個得到功能注釋；70 與330 d相比，有2 013個差異基因，其中有622個得到功能注釋； 165與220 d相比，有656個差異基因，其中169個 得到功能注釋；165與330 d相比，有1 541個差異基因，其中有480個得到功能注釋；220與 330 d相比，有1 593個差異基因，其中有474個得到功能注釋。差異基因主要集中在“生物過程”，其次是“細胞組分”，最后是“分子功能”。對比6個組合的差異基因功能富集GO term，其差異基因在“生物過程”分類中，細胞過程(cellular process)所占比例最多，其次是生物調控(biological regulation)和生物過程調控(regulation of biological process)，在“細胞組分”分類中，細胞(cell)占的比例最大，其次是細胞部分(cell part)和細胞膜等；在“分子功能”分類中，結合(binding)占的比例最多，其次是催化活性(catalytic activity)。此外，可以看到“生物過程”分類中與發育繁殖相關的有發育過程 (developmental process)、繁殖(reproduction)、繁殖過程(reproductive process)。在70與165 d相比，共涉及272個與發育繁殖相關的差異表達基因；70與220 d相比，共涉及269個與發育繁殖相關的差異表達基因；70與330 d相比，共涉及262個與發育繁殖相關的差異表達基因；165與 220 d相比，共涉及57個與發育繁殖相關的差異表達基因；165與330 d相比，共涉及194個與發育繁殖相關的差異表達基因；220 d 與330 d相比，共涉及188個與發育繁殖相關的差異表達基因。

通過對卵泡發育的不同生理階段的邊雞卵巢差異表達基因KEGG通路比較均可發現信號轉導、免疫系統以及信號分子與相互作用路徑上富集最多(表3)。

2.5 差異表達基因mRNA的實時定量PCR(qPCR)驗證

為進一步驗證轉錄組測序的可靠性，本研究通過qPCR隨機驗證了11個差異基因的表達。結果表明，所選基因的mRNA表達趨勢與轉錄組測序數據一致(圖4)，可以證明轉錄組測序的結果可靠。

3 討 論

家禽的就巢行為與卵巢發育有直接聯系，就巢行為的發生會引起產蛋性能的下降。轉錄組測序以它檢測范圍廣、測序通量高、檢測成本低、分辨率好和靈敏度高等特點，已經被廣泛應用[14]，在金定鴨的卵巢[15]、京海黃雞卵巢[16]中都有應用。本研究對邊雞卵巢組織的轉錄組測序，共鑒定出18 891個基因，其中已知基因18 063個，預測的新基因828個。此外，共檢測出21 299個新的轉錄本。本研究挖掘了邊雞卵巢組織轉錄組信息，為開展邊雞就巢性狀相關基因的研究及分子調控機制研究奠定基礎。

對所有樣本間的相關性進行層次聚類分析，用pearson法計算每兩個樣本間所有基因表達量的相關系數，相關系數越接近1，表明這兩個樣本間的基因表達的相似度越高。本研究發現，在165和220 d之間的相關系數為0.983，基因表達相似度最高，說明卵巢發育成熟后，開始排卵和排卵高峰兩個樣本卵巢組織基因表達相關性最強，與層次聚類分析表達譜結果一致。從邊雞不同生理階段卵巢組織差異表達基因數量顯示，存在上調基因和下調基因，165與220 d相比，差異表達基因的總量最少(656個)，與樣本間的相關性進行層次聚類分析結果一致，這些差異表達基因可能促進卵泡發育。

本研究使用DESeq2R包對轉錄組數據進行篩選，以樣本間的相關性和不同生理階段卵泡數為篩選基礎，找出了40個可能與卵泡發育相關的差異表達基因，qPCR隨機驗證所選基因的mRNA表達趨勢與轉錄組測序數據一致，證明了轉錄組測序的結果可靠。其中35個基因的表達量與發育卵泡數成正相關，5個基因表達量與發育卵泡數不成正相關(鐵氧還蛋白-1L(FDX1L)、N-Myc下游調控基因-4(NDRG4)、黑素皮質素受體(MC3R)、催乳素釋放激素受體(PRLHR)和血管活性腸肽(VIP))。FDX1L基因編碼一種183個氨基酸的線粒體蛋白Fdx2，它在鐵-硫簇生物發生過程中起重要作用[17]。它們通常通過半胱氨酸的巰基側鏈與假體蛋白結合，并作為電子受體或供體，在關鍵的細胞活動中發揮重要作用[18]。本試驗中發現，FDX1L在70、220、165和330 d邊雞卵巢中表達量依次遞減，表明該基因可能參與卵巢發育和調控排卵過程，作用機

圖4 qPCR驗證不同日齡邊雞卵巢中的差異表達基因Fig.4 Differentially expressed genes in the ovary of Bian chicken at different ages verified by the qPCR

理需進一步研究探討。NDRG4在多種癌癥中被認為與腫瘤細胞增殖能力相關[19-20]，NDRG4 表達下調后卵巢癌細胞增殖能力增高、凋亡減少，其機制可能與調節 Bcl-2/Bax 表達及 P38 /MAPK 信號通路活性有關[21-22]，NDRG4通過Akt/CREB激活促進肌原性分化[23]，支持激活細胞外信號調節激酶(ERK)促進神經元分化，并減少參與神經元生長的Elk-1的激活[24]。現已有報道稱Bcl-2與鵝就巢相關[25]，雖在本試驗中未發現Bcl-2在不同卵泡發育階段的邊雞卵巢中的表達量呈差異性，但發現NDRG4在邊雞卵巢中的表達量隨日齡增加呈遞減趨勢，NDRG4基因下調可能促進卵巢機能衰老，與上述研究中NDRG4基因下調后卵巢癌細胞增殖能力增高、凋亡減少的信號通路可能相關，卵巢機能衰老與卵巢癌細胞增殖能力的相關性需要進一步研究探討。黑素皮質素受體(melanocortin receptors，MCRs)是G-蛋白偶聯受體(G protein-coupled receptors，GPCRs)超家族成員之一[26]。人、鼠和雞MCRs由5個成員組成，分別是MC1R、MC2R、MC3R、MC4R和MC5R，其中，MC3R基因不僅在外周組織中有一定表達，而且在大腦中呈現高表達，因此，又被稱為神經性黑素皮質素受體基因[27-28]。現有研究發現，MC3R主要參與調控動物的體重和能量代謝的過程[29]，也有研究發現，MC3R在人的心、乳腺、骨骼肌和腎以及睪丸和卵巢中表達[30]，在這些組織中表達很可能與免疫功能相關[31]。有關MC3R表達對禽類繁殖性能的研究尚未見報道，本試驗中發現，MC3R基因在邊雞卵巢70和330 d表達量高于165和220 d的表達量，推測MC3R可能起到抑制卵泡發育的作用，能夠促進邊雞卵巢就巢性，這一發現為進一步研究卵泡發育及排卵調控機制奠定基礎。李國輝等[32]發現，催乳素(PRL)和神經肽Y(NPY)基因作為雞繁殖性能的候選基因，為多基因聚合育種技術在家禽遺傳改良方面的應用。在本試驗中雖未發現PRL，但發現催乳素釋放激素受體(PRLHR)在卵泡發育期即165和220 d低表達，分別為1.14和1.16，而在未發育期70 d和就巢期330 d表達量相對較高，分別為3.13和2.72，PRLHR與NPY受體(NPYR)結合，從而影響PRL釋放，進而調控邊雞的就巢性。Hosseinpour 等[33]發現，血管活性腸肽(VIP)及其受體(VIPR-1)的SNP與火雞產蛋性能相關，但未報道與就巢性相關。本試驗發現，在70、165和220 d的邊雞卵巢中VIP呈差異性表達，但在就巢期即330 d的邊雞卵巢中未發現VIP表達，該基因表達可能參與卵巢發育，與Hosseinpour 等[33]的結果一致。

本研究通過GO功能分析，發現邊雞不同日齡卵巢組織差異表達基因功能主要集中在“生物過程”。在“生物過程”中，與發育繁殖相關的通路有發育過程、繁殖和繁殖過程，富集在這些通路的差異表達基因可能參與卵泡的發育和排卵調控機制；而邊雞卵巢成熟開始排卵期間，即165與220 d相比，與發育繁殖相關的差異表達基因僅57個且數量最少，與本研究中卵巢基因差異表達相關性分析結果一致，這些少量的差異表達基因可能參與發育成熟卵泡數量的調控；220與330 d相比，共涉及188個與發育繁殖相關的差異表達基因，這些差異表達基因可能與邊雞就巢性的調控相關。對不同生理階段的邊雞卵巢差異表達基因KEGG通路比較，均可發現這些差異基因在信號轉導、免疫系統以及信號分子與相互作用的路徑上富集最多。通過揭示這些差異表達基因的功能，為多基因聚合育種改良邊雞就巢性提供了可能。

4 結 論

本研究通過構建基因表達譜，獲取邊雞卵巢組織轉錄組信息，篩選出40個與卵巢發育相關的差異表達基因，通過GO功能富集和KEGG分析，獲得了差異表達基因的功能分類和通路，為開展邊雞卵泡發育相關基因的研究及分子調控機制研究奠定基礎。