發(fā)情周期不同階段綿羊卵巢lncRNA的鑒定和功能分析

王夢瑤，翟振翰，趙 路，王賽橋，王玉琴

(河南科技大學動物科技學院，洛陽 471023)

綿羊在我國是一種重要的家畜，能夠為人們提供羊肉、羊奶、羊毛和羊皮等優(yōu)質(zhì)的畜產(chǎn)品。隨著生活水平的提高，羊肉消費量日漸增長，但我國現(xiàn)有綿羊存欄量和羊肉產(chǎn)量雖居前列，出欄量和個體產(chǎn)肉量卻低于世界平均水平，目前很難滿足國內(nèi)的需求[1-2]，這也使得人們對綿羊繁殖能力進行深入的研究與探索。卵巢是雌性動物重要的生殖器官，其主要功能是產(chǎn)生并排出能夠完成受精的健康卵子，以及分泌性激素、黃體酮等類固醇激素以維持雌性動物性征和周期性繁殖活動，對動物繁殖具有關(guān)鍵性作用。

lncRNA是一種長度超過200 nt的非編碼RNA分子(non-coding RNA，ncRNA)，在多種生物過程中發(fā)揮著多種作用，是繼mRNA和miRNA之后最熱門的分子生物學研究領(lǐng)域之一[3-4]。研究發(fā)現(xiàn)，哺乳動物基因組約70%被轉(zhuǎn)錄，但實際上能編碼蛋白質(zhì)的只有約2%，其余被轉(zhuǎn)錄為非編碼RNA，其中l(wèi)ncRNA被認為是動物中表達量最高的ncRNA之一[5]，最初被人們當作基因組轉(zhuǎn)錄過程中形成的“噪音”，并不具有生物學功能[6]。但隨后的研究逐漸表明lncRNA通過表觀遺傳調(diào)控、轉(zhuǎn)錄調(diào)控和轉(zhuǎn)錄后調(diào)控[7]等方面來調(diào)控基因表達，在細胞增殖分化[8]、信號轉(zhuǎn)導[9]、免疫調(diào)節(jié)[10]等多種生物學過程中發(fā)揮著重要作用。在繁殖方面，關(guān)于lncRNA的報道很多，李耀坤等[11]在高、低繁殖力川中黑山羊卵巢組織中發(fā)現(xiàn)ENSCHIG00000001479等lncRNA可能對山羊生殖發(fā)育過程中的卵巢發(fā)育及排卵等進程具有重要的調(diào)控作用,認為鑒定出的差異表達lncRNA與山羊產(chǎn)羔數(shù)有關(guān)。然后，F(xiàn)eng等[12]分別在高產(chǎn)和低產(chǎn)的湖羊卵巢中鑒定出5個lncRNA和76個mRNA。Su等[13]通過高通量測序篩選梅山豬和約克夏豬早期妊娠過程中差異表達lncRNA，得出XLOC-2222497及其靶標AKR1C1可以與豬子宮內(nèi)膜中的孕酮相互作用以控制妊娠維持。Chen等[14]在具有不同F(xiàn)ecB基因型的綿羊卵泡期和黃體期的下丘腦中鑒定到53個差異表達mRNAs和40個差異表達lncRNAs。這些研究證明了lncRNA在生殖組織中的存在和作用。

本研究利用高通量測序技術(shù)和生物信息分析方法對綿羊黃體期和卵泡期卵巢組織中的lncRNA進行測序，鑒定和預測不同時期差異表達基因lncRNA，通過GO和KEGG進行功能富集分析，篩選其中與繁殖相關(guān)的lncRNA及其靶基因參與的生物學過程，以期為探索lncRNA在綿羊繁育方面的作用提供理論基礎(chǔ)。

1 材料與方法

1.1 試驗動物和樣品采集

試驗動物為河南省洛寧農(nóng)本畜牧科技開發(fā)有限公司提供的綿羊，品種為湖羊。選取年齡在1.5～2.5歲，體重43 kg左右，體型均一，發(fā)育良好的母羊陰道放置陰道栓(60 mg甲羥孕酮)，12 d撤掉陰道栓。在試情兩次確認發(fā)情并撤栓后的24 h(卵泡期)和10 d(黃體期)屠宰采樣，2個時期綿羊各3只。將屠宰后采集的綿羊卵巢組織放在1.2 mL無酶凍存管內(nèi)，迅速放入液氮中保存并帶回，以備提取lncRNA。

1.2 卵巢組織總RNA的抽提

利用TRIzol法分別提取每個卵巢組織的總RNA，1.2%瓊脂糖凝膠電泳測定RNA完整性及是否存在DNA污染，Nanodorp 2000檢測RNA濃度。將測定合格的樣品送至北京諾禾致源公司進行RNA-seq。

1.3 文庫構(gòu)建、測序及轉(zhuǎn)錄本拼裝

按照要求進行l(wèi)ncRNA文庫構(gòu)建，庫檢合格后使用Illumina上機測序。高通量測序平臺Illumina測序得到的原始數(shù)據(jù)(Raw Reads)經(jīng)過質(zhì)控并去除低質(zhì)量、有污染、含有接頭等序列后的數(shù)據(jù)(Clean Reads)用作后續(xù)分析。使用HISAT2[15]對過濾后的reads進行參考基因組的比對分析。使用StringTie[16]對HISAT2比對的結(jié)果進行轉(zhuǎn)錄本的拼接，得到盡可能小的轉(zhuǎn)錄本集合，并對轉(zhuǎn)錄本進行定量分析。

1.4 差異表達lncRNA的篩選

基于轉(zhuǎn)錄組拼接結(jié)果，使用Cuffmerge軟件得到合并的轉(zhuǎn)錄本集合，根據(jù)lncRNA的結(jié)構(gòu)特點以及不編碼蛋白的功能特點，設(shè)置一系列嚴格的篩選條件，通過5步篩選：外顯子數(shù)篩選、轉(zhuǎn)錄本長度篩選、轉(zhuǎn)錄本已知注釋篩選、轉(zhuǎn)錄本表達量篩選、編碼潛能篩選，取CPC2、CNCI、PFAM三種工具預測結(jié)果的交集，作為本次分析預測得到的lncRNA數(shù)據(jù)集進行后續(xù)的分析。使用StringTie-eB軟件對包括mRNA、lncRNA以及TUCP在內(nèi)的轉(zhuǎn)錄本進行定量分析，得到各樣本轉(zhuǎn)錄本的FPKM信息。使用cuffdiff對不同類型轉(zhuǎn)錄本整體進行差異分析，篩選閾值設(shè)置為|log2(fold change)|≥2且q-value <0.05。

1.5 lncRNA靶向mRNA預測

lncRNA主要通過調(diào)控mRNA來實現(xiàn)功能。本研究通過lncRNA與蛋白編碼基因的位置關(guān)系(co-location)和表達相關(guān)性(co-expression)預測其生物學功能。將co-location的閾值設(shè)定為lncRNA上、下游100 kb，后續(xù)再通過對lncRNA共定位的mRNA基因進行功能富集分析預測lncRNA的主功能。采用Pearson相關(guān)系數(shù)法分析樣本間lncRNA與mRNA的相關(guān)性，取相關(guān)性絕對值大于0.95的mRNA基因進行功能富集分析，預測與mRNA共表達的lncRNA的主要功能。

1.6 GO富集和KEGG通路分析

Gene Ontology(簡稱 GO，www.geneontology.org)是基因功能國際標準分類體系，使用GOseq[17]對差異lncRNA的靶基因進行GO富集分析。KEGG (Kyoto Encyclopedia of Genes and Genomes) 是有關(guān)生物學通路的主要公共數(shù)據(jù)庫[18]，使用KOBAS(2.0)[19]進行通路富集分析，F(xiàn)DR≤0.05為顯著富集。

1.7 lncRNA的熒光定量PCR驗證

為驗證測序結(jié)果的可靠性，隨機選取5個差異表達的lncRNA進行qRT-PCR驗證。使用Primer 5.0設(shè)計引物，GAPDH作為內(nèi)參基因，引物均由北京擎科生物科技有限公司合成，引物序列信息如表1所示。PCR反應(yīng)體系(20 μL)：10 μL lncRNA PreMix，上、下游引物各0.5 μL，1 μL cDNA和8 μL RNase-Free ddH2O。PCR反應(yīng)條件：95 ℃ 3 min；95 ℃ 5 s，57 ℃ 10 s，72 ℃ 15 s，40個循環(huán)。使用2-ΔΔCt方法計算lncRNA的相對表達量，并與轉(zhuǎn)錄組結(jié)果進行對比驗證。

表1 lncRNA引物序列信息

2 結(jié) 果

2.1 測序數(shù)據(jù)的質(zhì)量及比對信息

對測序下機的Raw Reads進行質(zhì)控過濾，如表2所示，卵泡期湖羊卵巢(OAF)文庫和黃體期湖羊卵巢(OAL)文庫的Q20和Q30比例均在93%以上。對質(zhì)控過濾后的reads進行參考基因組的比對分析，如表3所示，測序數(shù)據(jù)與基因組比對率在87%以上，79%以上的序列在參考序列上有唯一比對位置。各項數(shù)據(jù)表明測序數(shù)據(jù)質(zhì)量較好，滿足后續(xù)分析的要求。

表2 數(shù)據(jù)產(chǎn)出質(zhì)量情況

表3 Reads與參考基因組比對情況

2.2 lncRNAs的篩選與鑒定

對OAF和OAL的轉(zhuǎn)錄本進行分步篩選，鑒定得到候選lncRNAs轉(zhuǎn)錄本數(shù)目為21 085個，其中包括221個已知lncRNAs和20 864個新lncRNAs。進一步對候選lncRNA進行分類，如圖1A所示，基因間型lncRNAs占46.5%，反義型lncRNAs占8.6%，內(nèi)含子型的lncRNAs占44.9%。

將lncRNA與mRNA進行結(jié)構(gòu)比較分析，得到lncRNA與mRNA的轉(zhuǎn)錄本長度、外顯子個數(shù)、ORF長度(圖1D)的對比情況。結(jié)果發(fā)現(xiàn)，lncRNA和mRNA的長度整體趨勢基本一致(圖1B)；lncRNA的外顯子個數(shù)多為2個，顯著低于mRNA(圖1C)，這與馬駿杰等[20]的研究結(jié)果一致；mRNA的ORF長度明顯大于lncRNA，這提示lncRNA在轉(zhuǎn)錄及轉(zhuǎn)錄后等表觀調(diào)控方面起著重要作用[21]。

2.3 lncRNAs的差異表達分析及驗證

利用cuffdiff進行差異分析得到初步差異分析的結(jié)果。設(shè)置篩選差異基因的標準為|log2(fold change)|≥2且q-value<0.05，得到的差異基因結(jié)果如圖2所示，在OAF和OAL中獲得1 379個差異表達的lncRNAs，其中1 158個表達上調(diào)，221個表達下調(diào)。隨機挑選5個差異表達lncRNAs進行qRT-PCR，驗證RNA測序數(shù)據(jù)的可靠性。驗證結(jié)果與測序數(shù)據(jù)結(jié)果基本一致(圖3)，表明測序數(shù)據(jù)結(jié)果可靠性高。

2.4 差異lncRNAs共定位基因的GO和KEGG分析

對1 379個差異表達的lncRNA的共定位靶基因進行生物功能分析，共顯著富集到3個GO條目(P<0.05)，包括1個分子功能(molecular function，MF)，2個生物學過程(biological process，BP)，沒有富集到細胞組分(cellular component，CC)。顯著富集到的GO條目如表4所示，在分子功能歸類中，差異基因顯著富集到谷胱甘肽轉(zhuǎn)移酶活性(glutathione transferase activity)。在生物學過程中，差異基因顯著富集到病毒過程的負調(diào)控(negative regulation of viral process)、病毒基因組復制的負調(diào)控(negative regulation of viral genome replication)。

A. lncRNA類型分布; B. lncRNA與mRNA外顯子數(shù)目；C. lncRNA與mRNA的轉(zhuǎn)錄本長度；D. lncRNA與mRNA ORF長度A. Distribution of lncRNA types; B. lncRNA and mRNA exon number; C. lncRNA and mRNA transcripts length; D. lncRNA and mRNA ORF length圖1 lncRNA類型分布及其與mRNA的結(jié)構(gòu)特征比較分析Fig.1 Distribution of lncRNA types and structural characteristics comparison with mRNA

圖2 差異表達lncRNAs的火山圖Fig.2 Volcano plot of differentially expressed lncRNAs

圖3 差異表達lncRNAs的qRT-PCR驗證Fig.3 Validation of the differentially expressed lncRNAs by qRT-PCR

表4 lncRNA的共定位基因前3個GO條目

在KEGG富集分析中，差異表達lncRNA的共定位靶基因涉及265條信號通路，其中共顯著富集到10條通路(P<0.05)，包括胞吞作用(endocytosis)、谷胱甘肽代謝(glutathione metabolism)、藥物代謝-細胞色素P450(drug metabolism-cytochrome P450)、精氨酸和脯氨酸代謝(arginine and proline metabolism)、細胞色素P450代謝外源物質(zhì)(metabo-lism of xenobiotics by cytochrome P450)等。P值前20的通路見圖4。

2.5 差異lncRNAs共表達基因的GO和KEGG分析

對1 379個差異表達lncRNAs的共表達靶基因進行生物功能分析，共顯著富集到482個GO條目(P<0.05)，包括63個分子功能，378個生物學過程，41個細胞組分。前10條顯著富集到的GO條目如表5所示，在分子功能歸類中，差異基因顯著富集到結(jié)合(binding)、蛋白結(jié)合(protein binding)等。在生物學過程中，差異基因顯著富集到卵泡發(fā)育(ovarian follicle development)、排卵周期過程(ovulation cycle process)等。在細胞組分中，差異基因顯著富集到質(zhì)膜(plasma membrane part)、神經(jīng)元投射(neuron projection)等。

在KEGG富集分析中，差異表達lncRNA的共表達靶基因涉及275條信號通路，其中共顯著富集到28條通路(P<0.05)，包括鈣離子信號通路(calcium signaling pathway)、卵母細胞減數(shù)分裂(oocyte meiosis)、cAMP信號通路(cAMP signaling pathway)、催產(chǎn)素信號通路(oxytocin signaling pathway)、MAPK信號通路(MAPK signaling pathway)、甲狀腺激素合成通路(thyroid hormone synthesis)、雌激素信號通路(estrogen signaling pathway)等。P值前20的通路見圖5。

2.6 繁殖相關(guān)的候選lncRNAs及其靶基因篩選

本研究KEGG富集分析得到了可能與繁殖有關(guān)的通路基因(表6)。這些通路包括鈣離子信號通路、卵母細胞減數(shù)分裂、cAMP信號通路、催產(chǎn)素信號通路、MAPK信號通路、甲狀腺激素合成通路、雌激素信號通路。

通過篩選KGEE富集到信號通路，篩選出5個顯著差異基因，與多個不同lncRNA存在靶標關(guān)系(表7)。其中多個lncRNAs靶向同一個差異基因，如LNC_003902、LNC_003906、LNC_003907均靶向PPP3CA等；也存在一個lncRNA同時具有多個靶基因，如ADCY5、ADCY1、CDC25A等都是LNC_012847的潛在目標，這些lncRNAs呈現(xiàn)出了多種方向的調(diào)控作用。

圖4 差異表達lncRNA共定位基因的KEGG分析Fig.4 KEGG analysis of differentially expressed lncRNA co-location genes

表5 lncRNA的共表達基因前10個GO條目

圖5 差異表達lncRNA共表達基因的KEGG分析Fig.5 KEGG analysis of differentially expressed lncRNA co-expression genes

表6 繁殖相關(guān)的候選信號通路及基因

表7 差異基因與lncRNAs的靶標關(guān)系

3 討 論

繁殖能力對綿羊的經(jīng)濟效益有重要影響。lncRNA因其在基因調(diào)控網(wǎng)絡(luò)和廣泛的生物學過程中的作用而受到廣泛關(guān)注。研究發(fā)現(xiàn)，lncRNA涉及多種動物繁殖相關(guān)的過程，如妊娠[22]、性腺激素應(yīng)答[23]、卵子成熟[24]和胎盤形成[25]等。本研究對綿羊不同發(fā)情周期卵巢轉(zhuǎn)錄組分析，選擇不同發(fā)情周期的卵巢篩選出的差異基因顯著富集到的通路及通路中的靶基因進行有關(guān)繁殖調(diào)控的討論。

本研究在KEGG富集分析得到的多個與繁殖調(diào)控相關(guān)的候選通路中選擇卵母細胞減數(shù)分裂、MAPK、甲狀腺激素合成、雌激素信號通路，并對通路及其中的靶基因進行繁殖相關(guān)調(diào)控的討論。

卵母細胞減數(shù)分裂被認為是形成單倍體配子所需的特殊細胞分裂形式，卵母細胞減數(shù)分裂的重新啟動在卵泡期中也起到關(guān)鍵作用[26]。卵母細胞減數(shù)分裂期間染色體分離的精確調(diào)節(jié)對哺乳動物的繁殖至關(guān)重要。PPP3CA(protein phosphatase 3 catalytic subunit alpha)是蛋白磷酸酶3的催化亞基A，廣泛存在于真核生物的多種細胞中[27]。研究發(fā)現(xiàn)，PPP3CA與山羊繁殖力性狀有關(guān)[28]，基因內(nèi)缺失突變顯著影響山羊產(chǎn)仔數(shù)[29]。在卵泡期顯著上調(diào)的LNC_003902、LNC_003906、LNC_003907靶向的PPP3CA在此通路中顯著富集。推測上述lncRNA可能通過調(diào)控PPP3CA影響綿羊產(chǎn)羔數(shù)。CDC23(cell division cycle 23)是小鼠卵母細胞減數(shù)分裂期間紡錘體組裝和染色體分離的關(guān)鍵調(diào)節(jié)因子[30]，這表明CDC23在減數(shù)分裂卵母細胞中具有獨特的作用。LNC_020278、LNC_020280、LNC_020281、LNC_020279靶向的CDC23在此通路中顯著富集。根據(jù)已知CDC23功能的研究，推測上述lncRNA可能在卵母顆粒細胞減數(shù)分裂中有重要作用。Liu等[31]在高、低繁殖力寒泊羊篩選的卵泡期差異表達lncRNA也顯著富集到此通路，與本研究結(jié)果相似。

MAPK信號通路是是真核生物信號傳遞網(wǎng)絡(luò)中的重要途徑之一，MAPK被一系列細胞外的刺激信號激活并介導信號從細胞膜向細胞核傳導，是調(diào)節(jié)細胞增殖、分化、凋亡和死亡等多種生理過程的關(guān)鍵信號通路[32-33]。MAPK信號通路被催乳素抑制，該通路的功能對正常卵泡發(fā)育至關(guān)重要[34]。MAPK1(mitogen-activated protein kinase 1)是MAP激酶家族的成員，在顆粒細胞中的磷酸化由LH受體激活引起并介導卵母細胞成熟[35-36]。本研究發(fā)現(xiàn)，在黃體期顯著上調(diào)的LNC_011239靶向的MAPK1在此通路顯著富集，根據(jù)靶基因的功能分析，推測LNC_011239可能在黃體期通過調(diào)控MAPK1影響卵母細胞成熟。

甲狀腺激素合成通路已被確定對脊椎動物胚胎發(fā)生和胎兒成熟至關(guān)重要。此外，甲狀腺素(T4)和三碘甲狀腺原氨酸(T3)對正常發(fā)育、生長和代謝穩(wěn)態(tài)至關(guān)重要[37]。有研究報道，大鼠在懷孕期間甲狀腺激素缺乏，導致產(chǎn)仔數(shù)減少[38]。雌激素是由卵巢和胎盤分泌的一種類固醇類激素，在動物發(fā)情、分娩等生殖功能方面有重要作用[39]，雌激素信號通路是指雌激素與細胞核或細胞膜中的雌激素受體(ER)結(jié)合后，作用于雌激素受體反應(yīng)元件調(diào)節(jié)基因表達。當雌激素與受體結(jié)合后，可以維持雌性個體第二性征，可以刺激和維持黃體功能，在雌性動物繁殖周期中對卵泡的生長發(fā)育起重要作用。ADCY1(adenylate cyclase 1)和ADCY5(adenylate cyclase 5)是腺苷酸環(huán)化酶家族的成員。ADCY1和ADCY5被發(fā)現(xiàn)是影響山羊和奶牛繁殖力的候選基因[40-41]。在卵泡期顯著上調(diào)的LNC_003902、LNC_003906、LNC_003907、LNC_012843、LNC_012847、LNC_021474、LNC_021475和LNC_012847靶向的ADCY1和ADCY5在這兩條信號通路均有顯著富集，提示LNC_012847可能與綿羊繁殖力存在一定相關(guān)性。

4 結(jié) 論

本研究采用高通量測序技術(shù)和生物信息學工具，對發(fā)情周期不同階段綿羊卵巢的lncRNAs進行了鑒定，在OAF和OAL中獲得1 379個差異表達的lncRNAs，其中1 158個表達上調(diào)，221個表達下調(diào)。共表達和共定位的GO注釋和KEGG富集分析差異lncRNAs的靶基因富集到多個繁殖相關(guān)信號通路中，其中LNC_011239、LNC_012847、LNC_003902、LNC_003906、LNC_003907等靶向的MAPK1、ADCY1、ADCY5、PPP3CA和CDC23可能發(fā)揮關(guān)鍵的調(diào)控作用。此外，這些差異lncRNAs表達譜為揭示綿羊繁殖能力的分子機制提供了有價值的資源。綿羊繁殖相關(guān)的候選lncRNAs及鑒定的關(guān)鍵靶基因可能需要通過敲除、過表達等方法進一步研究，以驗證lncRNA在綿羊卵巢中的特異功能。