基于高通量測序的山羊卵巢基質(zhì)、卵泡轉(zhuǎn)錄組分析

陸婷婷，鄒 嫻，楊鎮(zhèn)瑋，雷 放，陳奕業(yè)，柳廣斌，孫寶麗，劉德武，李耀坤

(1 華南農(nóng)業(yè)大學(xué) 動物科學(xué)學(xué)院，廣東 廣州 510642； 2 廣東省農(nóng)業(yè)科學(xué)院 動物科學(xué)研究所/畜禽育種國家重點實驗室/廣東省畜禽育種與營養(yǎng)研究重點實驗室，廣東 廣州 510640)

山羊Capra hircus繁殖力影響?zhàn)B羊業(yè)的經(jīng)濟效應(yīng)，受遺傳、營養(yǎng)、管理水平等因素影響，是復(fù)雜的數(shù)量性狀。母羊的產(chǎn)羔數(shù)是山羊繁殖性能之一，由排卵數(shù)決定，主要受遺傳因素影響，與卵泡發(fā)育息息相關(guān)[1]。卵泡發(fā)育是一個極其復(fù)雜的過程，在發(fā)育過程中絕大多數(shù)卵泡會閉鎖，只有少數(shù)卵泡能成熟并排卵[2-3]。

目前，關(guān)于整個卵巢的研究越來越多[4-5]，但卵巢作為一個動態(tài)發(fā)育的生殖器官，包含不同發(fā)育階段的卵泡，如原始卵泡、生長卵泡、成熟卵泡等。因此，在研究卵巢時，有必要對其細(xì)分后再開展研究。Terenina等[6]對豬卵巢健康小卵泡和閉鎖小卵泡轉(zhuǎn)錄組的研究發(fā)現(xiàn)CCDC80、GADD45A、ZNF628等基因在閉鎖小卵泡中的表達(dá)量超過健康小卵泡5倍，認(rèn)為可能是卵泡閉鎖的標(biāo)志基因。Zhang等[7]在豬健康卵泡和早期閉鎖卵泡的轉(zhuǎn)錄組研究中發(fā)現(xiàn)450個顯著差異表達(dá)基因。李鵬飛等[8-9]對牛的優(yōu)勢卵泡和從屬卵泡進(jìn)行轉(zhuǎn)錄組測序發(fā)現(xiàn)差異顯著的PRSS35、PTGRF和ARID4B等基因與卵泡發(fā)育相關(guān)。Hatzirodos等[10]在牛健康大卵泡和健康小卵泡的轉(zhuǎn)錄組研究中發(fā)現(xiàn)，在大卵泡中STAT4和XBP1基因表達(dá)上調(diào)，MGEA5基因下調(diào)；而KIT、IHH和MEST基因在小卵泡中最為活躍。吳陽升等[11]通過分析綿羊小卵泡和中卵泡轉(zhuǎn)錄組發(fā)現(xiàn)，在中卵泡表達(dá)上調(diào)的基因有IGFBP1、INHBA和TIMP1等，下調(diào)的有C-FOS、EGR1、FOSB等。綜上，豬、牛和綿羊不同發(fā)育階段卵泡的基因表達(dá)模式不同，然而山羊的卵泡轉(zhuǎn)錄組研究以及卵巢基質(zhì)和卵泡之間的基因表達(dá)模式鮮見報道。

川中黑山羊是我國優(yōu)秀的地方山羊品種[12]，經(jīng)過長期的自然和人工選擇，繁殖性能突出，遺傳性能穩(wěn)定，母羊平均產(chǎn)羔率為236.78%，羔羊成活率為91%[13]。目前對川中黑山羊繁殖性能的研究較少。轉(zhuǎn)錄組包括編碼RNA和非編碼RNA，指在某一特定條件下細(xì)胞或組織內(nèi)所有轉(zhuǎn)錄產(chǎn)物的總和，從整體水平研究基因結(jié)構(gòu)和功能[14]。轉(zhuǎn)錄組測序是通過高通量測序技術(shù)，對生物體全部基因的測序，反映生命體處于不同環(huán)境條件下的不同基因表達(dá)水平及調(diào)控模式[15]。本研究以發(fā)情期的川中黑山羊為研究對象，應(yīng)用高通量測序技術(shù)對卵巢基質(zhì)、大卵泡和小卵泡進(jìn)行轉(zhuǎn)錄組測序，比較分析三者的轉(zhuǎn)錄水平差異，探索卵巢基質(zhì)、大卵泡和小卵泡間的mRNA表達(dá)模式。旨在從轉(zhuǎn)錄組水平揭示卵巢不同部位的遺傳差異，豐富山羊卵巢基質(zhì)、大卵泡和小卵泡的轉(zhuǎn)錄組數(shù)據(jù)信息，篩選出的基因可為山羊卵泡發(fā)育提供一定的基因資源和新的研究方向。

1 材料與方法

1.1 試驗動物和樣品采集

選用南方某種羊場10頭健康的川中黑山羊母羊，年齡3.5～4.5歲。采用前列腺素1次注射法對每只母羊注射氯前列烯醇0.1 mg，24 h后第1次發(fā)情，18 d后確認(rèn)第2次發(fā)情，采用公羊試情法確認(rèn)母羊是否發(fā)情，以母羊搖尾、站立并接受爬跨視為發(fā)情[16-19]。發(fā)情后24 h內(nèi)進(jìn)行屠宰并采集卵巢。

卵巢用PBS緩沖液洗3次，分離小卵泡(Small follicle，SF，d＜3 mm)、大卵泡 (Large follicle，LF，d＞1 cm)和卵巢基質(zhì)，各6個重復(fù)組，其中1個小卵泡樣本包含8～10個小卵泡，1個大卵泡樣本僅包含1個大卵泡，連同分離后的卵巢基質(zhì)放入凍存管，立即置于液氮中，帶回實驗室于?80 ℃冰箱保存，用于后續(xù)試驗。試驗所得的卵巢基質(zhì)樣品中有少量d＜3 mm 的小卵泡，所有d＞1 cm 的大卵泡已分離出來。

1.2 總RNA抽提、cDNA文庫構(gòu)建和轉(zhuǎn)錄組測序

將卵巢基質(zhì)、大卵泡和小卵泡完全研磨，用Trizol試劑提取總RNA(Invitrogen，美國)。通過DanoDrop ND-2000 分光度計 (Thermo Science，美國)、10 g/L瓊脂糖凝膠分別檢測RNA的濃度和純度、完整性。每個樣品取3 μL總RNA，構(gòu)建cDNA文庫，利用Illumina HiSeq測序平臺進(jìn)行雙末端測序。整個流程委托上海派森諾生物科技股份公司完成。

1.3 轉(zhuǎn)錄組數(shù)據(jù)分析

對測序數(shù)據(jù)采用Cutadapt軟件去除3′端的接頭，去除的部分與接頭有至少10 bp重疊(AGATC GGAAG)，允許20%堿基錯配；去除平均占比低于Q20(堿基識別準(zhǔn)確率在99%以上的堿基所占百分比)的Reads。通過Bowtie2軟件(GCF_001704415.1_ARS1_genomic；http://www.ensembl.org/)建立參考基因組索引，使用Tophat2軟件的microexon-search和library-type=fr-firststrand參數(shù)將過濾后的Reads比對到參考基因組上，默認(rèn)Reads和參考基因組序列的堿基錯配在2個之內(nèi)即為比對成功。采用R語言DESeq程序包對差異基因進(jìn)行篩選，條件為：|log2表達(dá)差異倍數(shù)|＞1，顯著性P＜0.05，至少有4個樣品的FPKM(每百萬Reads中來自某一基因每千堿基長度的 Reads數(shù)目，F(xiàn)ragments per kilobase per million of reads)大于 1。通過 R 語言 Pheatmap軟件包對所有基因和樣品進(jìn)行雙向聚類分析。

1.4 差異基因功能注釋

利用Gene ontology(GO)數(shù)據(jù)庫對所有差異基因進(jìn)行注釋，通過 Kyoto encyclopedia of genes and genomes(KEGG)數(shù)據(jù)庫對差異基因進(jìn)行通路分析，采用超幾何分布計算差異基因顯著富集的GO條目和KEGG通路，P＜0.05為顯著富集。

1.5 qRT-PCR驗證分析

從差異基因中隨機選擇5個基因，利用qRTPCR進(jìn)行轉(zhuǎn)錄組表達(dá)譜驗證分析，PCR反應(yīng)體系10 μL：5 μL SYBR?Green Master Mix，上下游引物各 0.5 μL，1 μL cDNA 和 3 μL ddH2O。反應(yīng)條件：95 ℃，1 min；95 ℃，30 s，58 ℃，30 s，72 ℃，1 min，34個循環(huán)。引物設(shè)計參考NCBI的山羊相關(guān)基因序列，β-actin作為內(nèi)參基因，用Primer 5.0軟件設(shè)計，并送 Sangon Biotech 合成引物 (表1)。

表1 用于qRT-PCR檢測的基因引物Table 1 Gene primers for qRT-PCR detection

1.6 數(shù)據(jù)統(tǒng)計分析

試驗數(shù)據(jù)用SPSS 22.0統(tǒng)計軟件分析，通過獨立樣本t檢驗來進(jìn)行差異顯著性分析。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)控

將本研究所得的測序數(shù)據(jù)經(jīng)過質(zhì)量控制，對Clean reads進(jìn)行堿基含量分布和堿基質(zhì)量分布統(tǒng)計，Q30(堿基識別準(zhǔn)確率在99.9%以上的堿基所占百分比)在91%以上，測序數(shù)據(jù)與山羊的數(shù)據(jù)庫序列相似性在83%以上(表2)。綜合表明測序數(shù)據(jù)質(zhì)量好，滿足后續(xù)分析的要求。

2.2 差異基因篩選

測序結(jié)果分析表明，卵巢基質(zhì)和大卵泡之間顯著差異表達(dá)的基因總數(shù)為524個，其中顯著上調(diào)的有238個，顯著下調(diào)的有286個；卵巢基質(zhì)和小卵泡之間顯著差異表達(dá)的基因總數(shù)為180個，其中顯著上調(diào)的有25個，顯著下調(diào)的有155個；大卵泡和小卵泡之間顯著差異表達(dá)的基因總數(shù)為403個，其中顯著上調(diào)的有281個，顯著下調(diào)的有122個。結(jié)合近年哺乳動物中與卵泡發(fā)育相關(guān)的基因的報道[6,20-33]，我們發(fā)現(xiàn)本試驗中15個差異表達(dá)的基因可能與卵泡發(fā)育相關(guān)，包括INHA、INHBA、CYP19A1、MGST1、BEX2、APOA1、APOA2、MMP15、CITED1、PRSS23、LRP8、ADAMTS1、SPARCL1、TNFRSF1B和TNFRSF19。按照P值大小篩選排名前10的差異表達(dá)基因，再根據(jù)差異倍數(shù)從大到小進(jìn)行排序[34-35]，結(jié)果如表3所示。通過聚類圖(圖1)發(fā)現(xiàn)，卵巢基質(zhì)、大卵泡和小卵泡各樣本的基因分別聚類在一起，且同一組幾個不同樣本聚類在同一個簇中。此外，卵巢基質(zhì)與小卵泡間差異基因的數(shù)量最少，且大部分基因表達(dá)量相似(圖1)，推測小卵泡和卵巢基質(zhì)間的表達(dá)模式較為相似。大卵泡與卵巢基質(zhì)、小卵泡的差異表達(dá)基因數(shù)量較多，且表達(dá)量相差較大，推測大卵泡與卵巢基質(zhì)、小卵泡的基因表達(dá)模式差異較大。

表3 不同卵巢部位間排名前10位的差異表達(dá)基因Table 3 Top ten differentially expressed genes among different parts of ovary

圖1 不同卵巢部位間差異表達(dá)基因的聚類圖Fig. 1 Cluster analysis of differentially expressed genes among different parts of ovary

2.3 共有的差異表達(dá)基因

在卵巢基質(zhì)vs大卵泡、卵巢基質(zhì)vs小卵泡、小卵泡vs大卵泡這3組中有4個共有的差異表達(dá)基因；卵巢基質(zhì)vs大卵泡、卵巢基質(zhì)vs小卵泡之間有73個共有的差異表達(dá)基因；卵巢基質(zhì)vs大卵泡、小卵泡vs大卵泡之間有186個共有的差異表達(dá)基因；卵巢基質(zhì)vs小卵泡、小卵泡vs大卵泡之間有29個共有的差異表達(dá)基因。

2.4 差異表達(dá)基因的生物學(xué)功能分析

基于GO數(shù)據(jù)庫，對差異基因進(jìn)行功能分類，各功能分類下差異前10位的GO條目如表4、表5和表6所示。在卵巢基質(zhì)vs大卵泡(表4)中有15個GO條目顯著富集(P＜0.05)，包括催化活性、細(xì)胞質(zhì)、細(xì)胞膜、細(xì)胞內(nèi)膜結(jié)合細(xì)胞器、轉(zhuǎn)移酶活性等。在卵巢基質(zhì)vs小卵泡(表5)中有25個GO條目顯著富集(P＜0.05)，包括小分子結(jié)合、碳水化合物衍生物結(jié)合、核苷酸結(jié)合、磷酸核苷結(jié)合、嘌呤核苷三磷酸結(jié)合等。在小卵泡vs大卵泡(表6)中有64個GO條目顯著富集(P＜0.05)，包括細(xì)胞過程、細(xì)胞內(nèi)部、線粒體部分、嘌呤核苷酸結(jié)合等。

表4 卵巢基質(zhì)vs大卵泡部分差異表達(dá)基因的GO功能分類Table 4 GO functional classification of partial differentially expressed genes in ovarian stroma vs large follicle

表5 卵巢基質(zhì)vs小卵泡部分差異表達(dá)基因的GO功能分類Table 5 GO functional classification of partial differentially expressed genes in ovarian stroma vs small follicle

基于KEGG數(shù)據(jù)庫的通路富集分析(各組對比項目中前10位的信號通路見表7)，在卵巢基質(zhì)vs大卵泡中差異表達(dá)基因涉及55條通路，顯著富集的通路有46個(P＜0.05)，包括氧化磷酸化、代謝途徑、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工、檸檬酸循環(huán)、類固醇生物合成等。在卵巢基質(zhì)vs小卵泡中差異表達(dá)基因涉及27條通路，顯著富集的通路有17個(P＜0.05)，包括細(xì)胞周期、錯配修復(fù)、DNA復(fù)制、卵母細(xì)胞減數(shù)分裂、RNA轉(zhuǎn)運、孕酮介導(dǎo)的卵母細(xì)胞成熟等。在小卵泡vs大卵泡中差異表達(dá)基因涉及42條通路，顯著富集的通路有30個(P＜0.05)，包括氧化磷酸化、蛋白酶體、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工等。

表6 小卵泡vs大卵泡部分差異表達(dá)基因的GO功能分類Table 6 GO functional classification of partial differentially expressed genes in small follicle vs large follicle

表7 卵巢基質(zhì)、大卵泡和小卵泡部分差異表達(dá)基因的KEGG代謝通路Table 7 KEGG metabolic pathway of partial differentially expressed genes in ovarian stroma, large follicle and small follicle

2.5 基因篩選及qRT-PCR驗證分析

從3組差異基因中篩選5個基因，分別是CYP19A1和INHBA(卵巢基質(zhì)vs大卵泡)、GPNMB(卵巢基質(zhì)vs小卵泡)、KITLG和SFRP1(小卵泡vs大卵泡)，進(jìn)行qRT-PCR驗證(圖2)。這些基因的表達(dá)變化趨勢與測序結(jié)果基本一致，說明測序結(jié)果可靠。

圖2 不同卵巢部位間差異表達(dá)基因的qRT-PCR驗證Fig. 2 Quantitative identification of differentially expressed genes among different parts of ovary by qRT-PCR

3 討論

山羊性成熟時，以卵泡波的形式募集卵泡，并使卵泡發(fā)育至成熟排卵[36]。在形成有腔卵泡之前，卵泡的發(fā)育不受生殖激素調(diào)控，發(fā)育比較緩慢；在形成腔體時，性腺軸會分泌激素，促進(jìn)卵泡發(fā)育或閉鎖[37-38]。山羊是單胎或雙胎動物[3]，與人類相似，研究卵泡不同發(fā)育階段的作用機制，山羊是較好的動物模型。

本研究采用高通量測序技術(shù)，首次對川中黑山羊的卵巢基質(zhì)、大卵泡和小卵泡的mRNA表達(dá)量進(jìn)行研究。結(jié)果表明卵巢基質(zhì)、大卵泡和小卵泡的表達(dá)模式不同，其中卵巢基質(zhì)和小卵泡的表達(dá)模式更相近，可能是因為在分離卵泡時，每個卵巢基質(zhì)中d＞1 cm的大卵泡已分離出來，但不可避免的，卵巢基質(zhì)中含有少量的d＜3 mm的小卵泡，從而造成卵巢基質(zhì)和小卵泡的基因表達(dá)模式更相近。

通過差異基因分析，篩選出與其他哺乳動物卵泡發(fā)育或閉鎖有關(guān)的基因，包括在大卵泡上調(diào)的抑制素基因家族、細(xì)胞色素P450家族19亞家族成員1(CYP19A1)、MGST1、BEX家族、APO家族、MMP家族、CITED家族、PRSS23、LRP8；在大卵泡下調(diào)的ADAMTS1、SPARCL1、TNFRSF家族。在豬[6]和牛[8,11]的研究中，抑制素基因家族在閉鎖卵泡或小卵泡中低表達(dá)，通過負(fù)反饋調(diào)節(jié)抑制垂體分泌促卵泡激素，影響卵泡閉鎖。INHA和INHBA在牛大卵泡上調(diào)[10]，有助于卵泡膜產(chǎn)生雄激素并抑制垂體產(chǎn)生卵泡刺激素。CYP19A1的過表達(dá)增加大竇狀卵泡類固醇合成[6-8]，促進(jìn)卵泡發(fā)育。谷胱甘肽代謝MGST1與細(xì)胞的代謝和解毒能力有關(guān)，其高表達(dá)可能表示健康卵泡的解毒能力更好[39]。BEX家族是X染色體連鎖基因家族，BEX2在小卵泡低表達(dá)，引起線粒體凋亡，使細(xì)胞周期阻滯于G1期，從而調(diào)控顆粒細(xì)胞增殖和凋亡[9,40-41]。載脂蛋白APOA1在中卵泡中高表達(dá)，可能參與了卵泡的選擇及發(fā)育過程[11]?；|(zhì)金屬蛋白酶 (Matrix metallopeptidase，MMP)家族是一類依賴金屬鋅離子的蛋白水解酶家族,能有效地降解細(xì)胞外基質(zhì)，可能與卵泡體積的快速增大有關(guān)[11]。ADAMTS1基因家族與膜外結(jié)構(gòu)域的分裂有關(guān)，在小卵泡中表達(dá)上調(diào)[9]，在經(jīng)產(chǎn)奶牛排卵前的卵泡中上調(diào)[42]。SPARCL1是SPARC家族一員，在小卵泡中高表達(dá)，阻滯細(xì)胞停留在G1期，抑制細(xì)胞增殖[9,43]。但I(xiàn)GFBP5、IGF1、MMP9、PGRMC1等基因的表達(dá)量與其他研究報道的不一致。值得注意的是，在本研究中經(jīng)典顆粒細(xì)胞凋亡標(biāo)志物TNF受體超家族(TNFRSF1B、TNFRSF19)在大卵泡中下調(diào)[8]，與一些研究報道的無顯著差異不同[6-7]。

進(jìn)一步對差異基因功能分析發(fā)現(xiàn)，卵巢基質(zhì)vs大卵泡、卵巢基質(zhì)vs小卵泡、小卵泡vs大卵泡三者富集的通路大不相同，其中，卵巢基質(zhì)vs大卵泡、小卵泡vs大卵泡富集的通路較多，且有大部分是相同的；而卵巢基質(zhì)vs小卵泡富集的通路較少，與其他2組差異較大。說明可能是大卵泡的生物學(xué)功能較為豐富。功能富集表明不同組織表現(xiàn)出不同的富集特征，核苷酸結(jié)合和磷酸核苷結(jié)合過程在3組中都共同富集，但在卵巢基質(zhì)vs大卵泡中，特有的GO條目包括RNA加工過程、細(xì)胞質(zhì)、細(xì)胞內(nèi)膜、酶的催化活性等，可能是因為大卵泡階段代謝過程比較旺盛，大卵泡與卵巢基質(zhì)之間物質(zhì)運輸比較頻繁。在卵巢基質(zhì)vs小卵泡中特有的GO條目有細(xì)胞應(yīng)激反應(yīng)、初級代謝過程等，可能是因為小卵泡正處于初級發(fā)育階段。在小卵泡vs大卵泡中特有的GO條目有線粒體和線粒體膜、細(xì)胞器內(nèi)膜等，可能是因為卵泡發(fā)育過程主要受到線粒體膜的調(diào)控。

與卵巢基質(zhì)相比，檸檬酸循環(huán)和類固醇生物合成信號通路在大卵泡中上調(diào)。檸檬酸循環(huán)保證機體能量供給和促進(jìn)物質(zhì)代謝[44]。類固醇生物合成在卵泡發(fā)育中起著重要的作用[45]，例如類固醇合成快速調(diào)節(jié)蛋白 (Steroidogenic acute regulatory protein，STAR)可介導(dǎo)膽固醇進(jìn)入線粒體內(nèi)膜，調(diào)控孕酮和雌激素的合成[46-48]，影響卵泡發(fā)育；健康卵泡中CYP19A1的上調(diào)，促進(jìn)孕烯醇酮形成，孕烯醇酮是黃體酮(Progesterone，P4)和雌二醇(Estradiol，E2)的底物，可能是促進(jìn)卵泡發(fā)育的主要因素[7]。

在卵巢基質(zhì)vs小卵泡中，差異基因富集的信號通路與細(xì)胞分裂周期有關(guān)。有研究[49-50]報道，在小鼠卵泡發(fā)育階段，細(xì)胞周期、DNA復(fù)制和減數(shù)分裂等多個信號通路被激活，與本試驗結(jié)果一致，可能在小卵泡階段，卵母細(xì)胞發(fā)育上調(diào)了與細(xì)胞分裂相關(guān)的信號通路。

在小卵泡vs大卵泡中，差異基因富集的甲狀腺激素信號通路、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工、AMPK信號通路、類固醇生物合成和Notch信號通路在豬[7]和牛[10]的卵泡轉(zhuǎn)錄組研究中均有報道。內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工和甲狀腺激素信號通路的顯著富集，可能是因為隨著卵泡發(fā)育，顆粒細(xì)胞分化成熟，甾體激素合成逐漸增強[51-52]，物質(zhì)代謝也因此增強。在人和大鼠的研究中發(fā)現(xiàn)[53]，類固醇激素(如孕酮)與不同膜受體結(jié)合，抑制顆粒細(xì)胞的有絲分裂或凋亡，從而抑制卵泡發(fā)育或閉鎖，維持卵泡的平衡狀態(tài)。關(guān)于綿羊的研究發(fā)現(xiàn)，Notch信號通路調(diào)節(jié)細(xì)胞增殖、分化和凋亡，參與卵巢基質(zhì)卵泡的發(fā)育[54]。

4 結(jié)論

卵巢基質(zhì)、大卵泡、小卵泡的差異mRNA表達(dá)模式不相同，其中卵巢基質(zhì)與小卵泡更相近。就卵泡發(fā)育而言，深入研究這些特異性mRNA在卵泡不同發(fā)育階段的作用有著重要意義。本文從生物信息學(xué)角度分析卵巢不同組織表達(dá)譜，為下一步差異基因的功能驗證及其在卵泡發(fā)育中的作用提供一定的理論基礎(chǔ)。

致謝：感謝練志全幫助分離小卵泡和大卵泡，感謝趙智鋒幫助分析轉(zhuǎn)錄組測序結(jié)果！