王夢(mèng)瑤,翟振翰,趙 路,王賽橋,王玉琴
(河南科技大學(xué)動(dòng)物科技學(xué)院,洛陽(yáng) 471023)
綿羊在我國(guó)是一種重要的家畜,能夠?yàn)槿藗兲峁┭蛉狻⒀蚰?、羊毛和羊皮等?yōu)質(zhì)的畜產(chǎn)品。隨著生活水平的提高,羊肉消費(fèi)量日漸增長(zhǎng),但我國(guó)現(xiàn)有綿羊存欄量和羊肉產(chǎn)量雖居前列,出欄量和個(gè)體產(chǎn)肉量卻低于世界平均水平,目前很難滿足國(guó)內(nèi)的需求[1-2],這也使得人們對(duì)綿羊繁殖能力進(jìn)行深入的研究與探索。卵巢是雌性動(dòng)物重要的生殖器官,其主要功能是產(chǎn)生并排出能夠完成受精的健康卵子,以及分泌性激素、黃體酮等類(lèi)固醇激素以維持雌性動(dòng)物性征和周期性繁殖活動(dòng),對(duì)動(dòng)物繁殖具有關(guān)鍵性作用。
lncRNA是一種長(zhǎng)度超過(guò)200 nt的非編碼RNA分子(non-coding RNA,ncRNA),在多種生物過(guò)程中發(fā)揮著多種作用,是繼mRNA和miRNA之后最熱門(mén)的分子生物學(xué)研究領(lǐng)域之一[3-4]。研究發(fā)現(xiàn),哺乳動(dòng)物基因組約70%被轉(zhuǎn)錄,但實(shí)際上能編碼蛋白質(zhì)的只有約2%,其余被轉(zhuǎn)錄為非編碼RNA,其中l(wèi)ncRNA被認(rèn)為是動(dòng)物中表達(dá)量最高的ncRNA之一[5],最初被人們當(dāng)作基因組轉(zhuǎn)錄過(guò)程中形成的“噪音”,并不具有生物學(xué)功能[6]。但隨后的研究逐漸表明lncRNA通過(guò)表觀遺傳調(diào)控、轉(zhuǎn)錄調(diào)控和轉(zhuǎn)錄后調(diào)控[7]等方面來(lái)調(diào)控基因表達(dá),在細(xì)胞增殖分化[8]、信號(hào)轉(zhuǎn)導(dǎo)[9]、免疫調(diào)節(jié)[10]等多種生物學(xué)過(guò)程中發(fā)揮著重要作用。在繁殖方面,關(guān)于lncRNA的報(bào)道很多,李耀坤等[11]在高、低繁殖力川中黑山羊卵巢組織中發(fā)現(xiàn)ENSCHIG00000001479等lncRNA可能對(duì)山羊生殖發(fā)育過(guò)程中的卵巢發(fā)育及排卵等進(jìn)程具有重要的調(diào)控作用,認(rèn)為鑒定出的差異表達(dá)lncRNA與山羊產(chǎn)羔數(shù)有關(guān)。然后,F(xiàn)eng等[12]分別在高產(chǎn)和低產(chǎn)的湖羊卵巢中鑒定出5個(gè)lncRNA和76個(gè)mRNA。Su等[13]通過(guò)高通量測(cè)序篩選梅山豬和約克夏豬早期妊娠過(guò)程中差異表達(dá)lncRNA,得出XLOC-2222497及其靶標(biāo)AKR1C1可以與豬子宮內(nèi)膜中的孕酮相互作用以控制妊娠維持。Chen等[14]在具有不同F(xiàn)ecB基因型的綿羊卵泡期和黃體期的下丘腦中鑒定到53個(gè)差異表達(dá)mRNAs和40個(gè)差異表達(dá)lncRNAs。這些研究證明了lncRNA在生殖組織中的存在和作用。
本研究利用高通量測(cè)序技術(shù)和生物信息分析方法對(duì)綿羊黃體期和卵泡期卵巢組織中的lncRNA進(jìn)行測(cè)序,鑒定和預(yù)測(cè)不同時(shí)期差異表達(dá)基因lncRNA,通過(guò)GO和KEGG進(jìn)行功能富集分析,篩選其中與繁殖相關(guān)的lncRNA及其靶基因參與的生物學(xué)過(guò)程,以期為探索lncRNA在綿羊繁育方面的作用提供理論基礎(chǔ)。
試驗(yàn)動(dòng)物為河南省洛寧農(nóng)本畜牧科技開(kāi)發(fā)有限公司提供的綿羊,品種為湖羊。選取年齡在1.5~2.5歲,體重43 kg左右,體型均一,發(fā)育良好的母羊陰道放置陰道栓(60 mg甲羥孕酮),12 d撤掉陰道栓。在試情兩次確認(rèn)發(fā)情并撤栓后的24 h(卵泡期)和10 d(黃體期)屠宰采樣,2個(gè)時(shí)期綿羊各3只。將屠宰后采集的綿羊卵巢組織放在1.2 mL無(wú)酶凍存管內(nèi),迅速放入液氮中保存并帶回,以備提取lncRNA。
利用TRIzol法分別提取每個(gè)卵巢組織的總RNA,1.2%瓊脂糖凝膠電泳測(cè)定RNA完整性及是否存在DNA污染,Nanodorp 2000檢測(cè)RNA濃度。將測(cè)定合格的樣品送至北京諾禾致源公司進(jìn)行RNA-seq。
按照要求進(jìn)行l(wèi)ncRNA文庫(kù)構(gòu)建,庫(kù)檢合格后使用Illumina上機(jī)測(cè)序。高通量測(cè)序平臺(tái)Illumina測(cè)序得到的原始數(shù)據(jù)(Raw Reads)經(jīng)過(guò)質(zhì)控并去除低質(zhì)量、有污染、含有接頭等序列后的數(shù)據(jù)(Clean Reads)用作后續(xù)分析。使用HISAT2[15]對(duì)過(guò)濾后的reads進(jìn)行參考基因組的比對(duì)分析。使用StringTie[16]對(duì)HISAT2比對(duì)的結(jié)果進(jìn)行轉(zhuǎn)錄本的拼接,得到盡可能小的轉(zhuǎn)錄本集合,并對(duì)轉(zhuǎn)錄本進(jìn)行定量分析。
基于轉(zhuǎn)錄組拼接結(jié)果,使用Cuffmerge軟件得到合并的轉(zhuǎn)錄本集合,根據(jù)lncRNA的結(jié)構(gòu)特點(diǎn)以及不編碼蛋白的功能特點(diǎn),設(shè)置一系列嚴(yán)格的篩選條件,通過(guò)5步篩選:外顯子數(shù)篩選、轉(zhuǎn)錄本長(zhǎng)度篩選、轉(zhuǎn)錄本已知注釋篩選、轉(zhuǎn)錄本表達(dá)量篩選、編碼潛能篩選,取CPC2、CNCI、PFAM三種工具預(yù)測(cè)結(jié)果的交集,作為本次分析預(yù)測(cè)得到的lncRNA數(shù)據(jù)集進(jìn)行后續(xù)的分析。使用StringTie-eB軟件對(duì)包括mRNA、lncRNA以及TUCP在內(nèi)的轉(zhuǎn)錄本進(jìn)行定量分析,得到各樣本轉(zhuǎn)錄本的FPKM信息。使用cuffdiff對(duì)不同類(lèi)型轉(zhuǎn)錄本整體進(jìn)行差異分析,篩選閾值設(shè)置為|log2(fold change)|≥2且q-value <0.05。
lncRNA主要通過(guò)調(diào)控mRNA來(lái)實(shí)現(xiàn)功能。本研究通過(guò)lncRNA與蛋白編碼基因的位置關(guān)系(co-location)和表達(dá)相關(guān)性(co-expression)預(yù)測(cè)其生物學(xué)功能。將co-location的閾值設(shè)定為lncRNA上、下游100 kb,后續(xù)再通過(guò)對(duì)lncRNA共定位的mRNA基因進(jìn)行功能富集分析預(yù)測(cè)lncRNA的主功能。采用Pearson相關(guān)系數(shù)法分析樣本間lncRNA與mRNA的相關(guān)性,取相關(guān)性絕對(duì)值大于0.95的mRNA基因進(jìn)行功能富集分析,預(yù)測(cè)與mRNA共表達(dá)的lncRNA的主要功能。
Gene Ontology(簡(jiǎn)稱(chēng) GO,www.geneontology.org)是基因功能?chē)?guó)際標(biāo)準(zhǔn)分類(lèi)體系,使用GOseq[17]對(duì)差異lncRNA的靶基因進(jìn)行GO富集分析。KEGG (Kyoto Encyclopedia of Genes and Genomes) 是有關(guān)生物學(xué)通路的主要公共數(shù)據(jù)庫(kù)[18],使用KOBAS(2.0)[19]進(jìn)行通路富集分析,F(xiàn)DR≤0.05為顯著富集。
為驗(yàn)證測(cè)序結(jié)果的可靠性,隨機(jī)選取5個(gè)差異表達(dá)的lncRNA進(jìn)行qRT-PCR驗(yàn)證。使用Primer 5.0設(shè)計(jì)引物,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個(gè)循環(huán)。使用2-ΔΔCt方法計(jì)算lncRNA的相對(duì)表達(dá)量,并與轉(zhuǎn)錄組結(jié)果進(jìn)行對(duì)比驗(yàn)證。
表1 lncRNA引物序列信息
對(duì)測(cè)序下機(jī)的Raw Reads進(jìn)行質(zhì)控過(guò)濾,如表2所示,卵泡期湖羊卵巢(OAF)文庫(kù)和黃體期湖羊卵巢(OAL)文庫(kù)的Q20和Q30比例均在93%以上。對(duì)質(zhì)控過(guò)濾后的reads進(jìn)行參考基因組的比對(duì)分析,如表3所示,測(cè)序數(shù)據(jù)與基因組比對(duì)率在87%以上,79%以上的序列在參考序列上有唯一比對(duì)位置。各項(xiàng)數(shù)據(jù)表明測(cè)序數(shù)據(jù)質(zhì)量較好,滿足后續(xù)分析的要求。
表2 數(shù)據(jù)產(chǎn)出質(zhì)量情況
表3 Reads與參考基因組比對(duì)情況
對(duì)OAF和OAL的轉(zhuǎn)錄本進(jìn)行分步篩選,鑒定得到候選lncRNAs轉(zhuǎn)錄本數(shù)目為21 085個(gè),其中包括221個(gè)已知lncRNAs和20 864個(gè)新lncRNAs。進(jìn)一步對(duì)候選lncRNA進(jìn)行分類(lèi),如圖1A所示,基因間型lncRNAs占46.5%,反義型lncRNAs占8.6%,內(nèi)含子型的lncRNAs占44.9%。
將lncRNA與mRNA進(jìn)行結(jié)構(gòu)比較分析,得到lncRNA與mRNA的轉(zhuǎn)錄本長(zhǎng)度、外顯子個(gè)數(shù)、ORF長(zhǎng)度(圖1D)的對(duì)比情況。結(jié)果發(fā)現(xiàn),lncRNA和mRNA的長(zhǎng)度整體趨勢(shì)基本一致(圖1B);lncRNA的外顯子個(gè)數(shù)多為2個(gè),顯著低于mRNA(圖1C),這與馬駿杰等[20]的研究結(jié)果一致;mRNA的ORF長(zhǎng)度明顯大于lncRNA,這提示lncRNA在轉(zhuǎn)錄及轉(zhuǎn)錄后等表觀調(diào)控方面起著重要作用[21]。
利用cuffdiff進(jìn)行差異分析得到初步差異分析的結(jié)果。設(shè)置篩選差異基因的標(biāo)準(zhǔn)為|log2(fold change)|≥2且q-value<0.05,得到的差異基因結(jié)果如圖2所示,在OAF和OAL中獲得1 379個(gè)差異表達(dá)的lncRNAs,其中1 158個(gè)表達(dá)上調(diào),221個(gè)表達(dá)下調(diào)。隨機(jī)挑選5個(gè)差異表達(dá)lncRNAs進(jìn)行qRT-PCR,驗(yàn)證RNA測(cè)序數(shù)據(jù)的可靠性。驗(yàn)證結(jié)果與測(cè)序數(shù)據(jù)結(jié)果基本一致(圖3),表明測(cè)序數(shù)據(jù)結(jié)果可靠性高。
對(duì)1 379個(gè)差異表達(dá)的lncRNA的共定位靶基因進(jìn)行生物功能分析,共顯著富集到3個(gè)GO條目(P<0.05),包括1個(gè)分子功能(molecular function,MF),2個(gè)生物學(xué)過(guò)程(biological process,BP),沒(méi)有富集到細(xì)胞組分(cellular component,CC)。顯著富集到的GO條目如表4所示,在分子功能歸類(lèi)中,差異基因顯著富集到谷胱甘肽轉(zhuǎn)移酶活性(glutathione transferase activity)。在生物學(xué)過(guò)程中,差異基因顯著富集到病毒過(guò)程的負(fù)調(diào)控(negative regulation of viral process)、病毒基因組復(fù)制的負(fù)調(diào)控(negative regulation of viral genome replication)。
A. lncRNA類(lèi)型分布; B. lncRNA與mRNA外顯子數(shù)目;C. lncRNA與mRNA的轉(zhuǎn)錄本長(zhǎng)度;D. lncRNA與mRNA ORF長(zhǎng)度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類(lèi)型分布及其與mRNA的結(jié)構(gòu)特征比較分析Fig.1 Distribution of lncRNA types and structural characteristics comparison with mRNA
圖2 差異表達(dá)lncRNAs的火山圖Fig.2 Volcano plot of differentially expressed lncRNAs
圖3 差異表達(dá)lncRNAs的qRT-PCR驗(yàn)證Fig.3 Validation of the differentially expressed lncRNAs by qRT-PCR
表4 lncRNA的共定位基因前3個(gè)GO條目
在KEGG富集分析中,差異表達(dá)lncRNA的共定位靶基因涉及265條信號(hào)通路,其中共顯著富集到10條通路(P<0.05),包括胞吞作用(endocytosis)、谷胱甘肽代謝(glutathione metabolism)、藥物代謝-細(xì)胞色素P450(drug metabolism-cytochrome P450)、精氨酸和脯氨酸代謝(arginine and proline metabolism)、細(xì)胞色素P450代謝外源物質(zhì)(metabo-lism of xenobiotics by cytochrome P450)等。P值前20的通路見(jiàn)圖4。
對(duì)1 379個(gè)差異表達(dá)lncRNAs的共表達(dá)靶基因進(jìn)行生物功能分析,共顯著富集到482個(gè)GO條目(P<0.05),包括63個(gè)分子功能,378個(gè)生物學(xué)過(guò)程,41個(gè)細(xì)胞組分。前10條顯著富集到的GO條目如表5所示,在分子功能歸類(lèi)中,差異基因顯著富集到結(jié)合(binding)、蛋白結(jié)合(protein binding)等。在生物學(xué)過(guò)程中,差異基因顯著富集到卵泡發(fā)育(ovarian follicle development)、排卵周期過(guò)程(ovulation cycle process)等。在細(xì)胞組分中,差異基因顯著富集到質(zhì)膜(plasma membrane part)、神經(jīng)元投射(neuron projection)等。
在KEGG富集分析中,差異表達(dá)lncRNA的共表達(dá)靶基因涉及275條信號(hào)通路,其中共顯著富集到28條通路(P<0.05),包括鈣離子信號(hào)通路(calcium signaling pathway)、卵母細(xì)胞減數(shù)分裂(oocyte meiosis)、cAMP信號(hào)通路(cAMP signaling pathway)、催產(chǎn)素信號(hào)通路(oxytocin signaling pathway)、MAPK信號(hào)通路(MAPK signaling pathway)、甲狀腺激素合成通路(thyroid hormone synthesis)、雌激素信號(hào)通路(estrogen signaling pathway)等。P值前20的通路見(jiàn)圖5。
本研究KEGG富集分析得到了可能與繁殖有關(guān)的通路基因(表6)。這些通路包括鈣離子信號(hào)通路、卵母細(xì)胞減數(shù)分裂、cAMP信號(hào)通路、催產(chǎn)素信號(hào)通路、MAPK信號(hào)通路、甲狀腺激素合成通路、雌激素信號(hào)通路。
通過(guò)篩選KGEE富集到信號(hào)通路,篩選出5個(gè)顯著差異基因,與多個(gè)不同lncRNA存在靶標(biāo)關(guān)系(表7)。其中多個(gè)lncRNAs靶向同一個(gè)差異基因,如LNC_003902、LNC_003906、LNC_003907均靶向PPP3CA等;也存在一個(gè)lncRNA同時(shí)具有多個(gè)靶基因,如ADCY5、ADCY1、CDC25A等都是LNC_012847的潛在目標(biāo),這些lncRNAs呈現(xiàn)出了多種方向的調(diào)控作用。
圖4 差異表達(dá)lncRNA共定位基因的KEGG分析Fig.4 KEGG analysis of differentially expressed lncRNA co-location genes
表5 lncRNA的共表達(dá)基因前10個(gè)GO條目
圖5 差異表達(dá)lncRNA共表達(dá)基因的KEGG分析Fig.5 KEGG analysis of differentially expressed lncRNA co-expression genes
表6 繁殖相關(guān)的候選信號(hào)通路及基因
表7 差異基因與lncRNAs的靶標(biāo)關(guān)系
繁殖能力對(duì)綿羊的經(jīng)濟(jì)效益有重要影響。lncRNA因其在基因調(diào)控網(wǎng)絡(luò)和廣泛的生物學(xué)過(guò)程中的作用而受到廣泛關(guān)注。研究發(fā)現(xiàn),lncRNA涉及多種動(dòng)物繁殖相關(guān)的過(guò)程,如妊娠[22]、性腺激素應(yīng)答[23]、卵子成熟[24]和胎盤(pán)形成[25]等。本研究對(duì)綿羊不同發(fā)情周期卵巢轉(zhuǎn)錄組分析,選擇不同發(fā)情周期的卵巢篩選出的差異基因顯著富集到的通路及通路中的靶基因進(jìn)行有關(guān)繁殖調(diào)控的討論。
本研究在KEGG富集分析得到的多個(gè)與繁殖調(diào)控相關(guān)的候選通路中選擇卵母細(xì)胞減數(shù)分裂、MAPK、甲狀腺激素合成、雌激素信號(hào)通路,并對(duì)通路及其中的靶基因進(jìn)行繁殖相關(guān)調(diào)控的討論。
卵母細(xì)胞減數(shù)分裂被認(rèn)為是形成單倍體配子所需的特殊細(xì)胞分裂形式,卵母細(xì)胞減數(shù)分裂的重新啟動(dòng)在卵泡期中也起到關(guān)鍵作用[26]。卵母細(xì)胞減數(shù)分裂期間染色體分離的精確調(diào)節(jié)對(duì)哺乳動(dòng)物的繁殖至關(guān)重要。PPP3CA(protein phosphatase 3 catalytic subunit alpha)是蛋白磷酸酶3的催化亞基A,廣泛存在于真核生物的多種細(xì)胞中[27]。研究發(fā)現(xiàn),PPP3CA與山羊繁殖力性狀有關(guān)[28],基因內(nèi)缺失突變顯著影響山羊產(chǎn)仔數(shù)[29]。在卵泡期顯著上調(diào)的LNC_003902、LNC_003906、LNC_003907靶向的PPP3CA在此通路中顯著富集。推測(cè)上述lncRNA可能通過(guò)調(diào)控PPP3CA影響綿羊產(chǎn)羔數(shù)。CDC23(cell division cycle 23)是小鼠卵母細(xì)胞減數(shù)分裂期間紡錘體組裝和染色體分離的關(guān)鍵調(diào)節(jié)因子[30],這表明CDC23在減數(shù)分裂卵母細(xì)胞中具有獨(dú)特的作用。LNC_020278、LNC_020280、LNC_020281、LNC_020279靶向的CDC23在此通路中顯著富集。根據(jù)已知CDC23功能的研究,推測(cè)上述lncRNA可能在卵母顆粒細(xì)胞減數(shù)分裂中有重要作用。Liu等[31]在高、低繁殖力寒泊羊篩選的卵泡期差異表達(dá)lncRNA也顯著富集到此通路,與本研究結(jié)果相似。
MAPK信號(hào)通路是是真核生物信號(hào)傳遞網(wǎng)絡(luò)中的重要途徑之一,MAPK被一系列細(xì)胞外的刺激信號(hào)激活并介導(dǎo)信號(hào)從細(xì)胞膜向細(xì)胞核傳導(dǎo),是調(diào)節(jié)細(xì)胞增殖、分化、凋亡和死亡等多種生理過(guò)程的關(guān)鍵信號(hào)通路[32-33]。MAPK信號(hào)通路被催乳素抑制,該通路的功能對(duì)正常卵泡發(fā)育至關(guān)重要[34]。MAPK1(mitogen-activated protein kinase 1)是MAP激酶家族的成員,在顆粒細(xì)胞中的磷酸化由LH受體激活引起并介導(dǎo)卵母細(xì)胞成熟[35-36]。本研究發(fā)現(xiàn),在黃體期顯著上調(diào)的LNC_011239靶向的MAPK1在此通路顯著富集,根據(jù)靶基因的功能分析,推測(cè)LNC_011239可能在黃體期通過(guò)調(diào)控MAPK1影響卵母細(xì)胞成熟。
甲狀腺激素合成通路已被確定對(duì)脊椎動(dòng)物胚胎發(fā)生和胎兒成熟至關(guān)重要。此外,甲狀腺素(T4)和三碘甲狀腺原氨酸(T3)對(duì)正常發(fā)育、生長(zhǎng)和代謝穩(wěn)態(tài)至關(guān)重要[37]。有研究報(bào)道,大鼠在懷孕期間甲狀腺激素缺乏,導(dǎo)致產(chǎn)仔數(shù)減少[38]。雌激素是由卵巢和胎盤(pán)分泌的一種類(lèi)固醇類(lèi)激素,在動(dòng)物發(fā)情、分娩等生殖功能方面有重要作用[39],雌激素信號(hào)通路是指雌激素與細(xì)胞核或細(xì)胞膜中的雌激素受體(ER)結(jié)合后,作用于雌激素受體反應(yīng)元件調(diào)節(jié)基因表達(dá)。當(dāng)雌激素與受體結(jié)合后,可以維持雌性個(gè)體第二性征,可以刺激和維持黃體功能,在雌性動(dòng)物繁殖周期中對(duì)卵泡的生長(zhǎng)發(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在這兩條信號(hào)通路均有顯著富集,提示LNC_012847可能與綿羊繁殖力存在一定相關(guān)性。
本研究采用高通量測(cè)序技術(shù)和生物信息學(xué)工具,對(duì)發(fā)情周期不同階段綿羊卵巢的lncRNAs進(jìn)行了鑒定,在OAF和OAL中獲得1 379個(gè)差異表達(dá)的lncRNAs,其中1 158個(gè)表達(dá)上調(diào),221個(gè)表達(dá)下調(diào)。共表達(dá)和共定位的GO注釋和KEGG富集分析差異lncRNAs的靶基因富集到多個(gè)繁殖相關(guān)信號(hào)通路中,其中LNC_011239、LNC_012847、LNC_003902、LNC_003906、LNC_003907等靶向的MAPK1、ADCY1、ADCY5、PPP3CA和CDC23可能發(fā)揮關(guān)鍵的調(diào)控作用。此外,這些差異lncRNAs表達(dá)譜為揭示綿羊繁殖能力的分子機(jī)制提供了有價(jià)值的資源。綿羊繁殖相關(guān)的候選lncRNAs及鑒定的關(guān)鍵靶基因可能需要通過(guò)敲除、過(guò)表達(dá)等方法進(jìn)一步研究,以驗(yàn)證lncRNA在綿羊卵巢中的特異功能。