王夢瑤,翟振翰,趙 路,王賽橋,王玉琴
(河南科技大學動物科技學院,洛陽 471023)
綿羊在我國是一種重要的家畜,能夠為人們提供羊肉、羊奶、羊毛和羊皮等優質的畜產品。隨著生活水平的提高,羊肉消費量日漸增長,但我國現有綿羊存欄量和羊肉產量雖居前列,出欄量和個體產肉量卻低于世界平均水平,目前很難滿足國內的需求[1-2],這也使得人們對綿羊繁殖能力進行深入的研究與探索。卵巢是雌性動物重要的生殖器官,其主要功能是產生并排出能夠完成受精的健康卵子,以及分泌性激素、黃體酮等類固醇激素以維持雌性動物性征和周期性繁殖活動,對動物繁殖具有關鍵性作用。
lncRNA是一種長度超過200 nt的非編碼RNA分子(non-coding RNA,ncRNA),在多種生物過程中發揮著多種作用,是繼mRNA和miRNA之后最熱門的分子生物學研究領域之一[3-4]。研究發現,哺乳動物基因組約70%被轉錄,但實際上能編碼蛋白質的只有約2%,其余被轉錄為非編碼RNA,其中lncRNA被認為是動物中表達量最高的ncRNA之一[5],最初被人們當作基因組轉錄過程中形成的“噪音”,并不具有生物學功能[6]。但隨后的研究逐漸表明lncRNA通過表觀遺傳調控、轉錄調控和轉錄后調控[7]等方面來調控基因表達,在細胞增殖分化[8]、信號轉導[9]、免疫調節[10]等多種生物學過程中發揮著重要作用。在繁殖方面,關于lncRNA的報道很多,李耀坤等[11]在高、低繁殖力川中黑山羊卵巢組織中發現ENSCHIG00000001479等lncRNA可能對山羊生殖發育過程中的卵巢發育及排卵等進程具有重要的調控作用,認為鑒定出的差異表達lncRNA與山羊產羔數有關。然后,Feng等[12]分別在高產和低產的湖羊卵巢中鑒定出5個lncRNA和76個mRNA。Su等[13]通過高通量測序篩選梅山豬和約克夏豬早期妊娠過程中差異表達lncRNA,得出XLOC-2222497及其靶標AKR1C1可以與豬子宮內膜中的孕酮相互作用以控制妊娠維持。Chen等[14]在具有不同FecB基因型的綿羊卵泡期和黃體期的下丘腦中鑒定到53個差異表達mRNAs和40個差異表達lncRNAs。這些研究證明了lncRNA在生殖組織中的存在和作用。
本研究利用高通量測序技術和生物信息分析方法對綿羊黃體期和卵泡期卵巢組織中的lncRNA進行測序,鑒定和預測不同時期差異表達基因lncRNA,通過GO和KEGG進行功能富集分析,篩選其中與繁殖相關的lncRNA及其靶基因參與的生物學過程,以期為探索lncRNA在綿羊繁育方面的作用提供理論基礎。
試驗動物為河南省洛寧農本畜牧科技開發有限公司提供的綿羊,品種為湖羊。選取年齡在1.5~2.5歲,體重43 kg左右,體型均一,發育良好的母羊陰道放置陰道栓(60 mg甲羥孕酮),12 d撤掉陰道栓。在試情兩次確認發情并撤栓后的24 h(卵泡期)和10 d(黃體期)屠宰采樣,2個時期綿羊各3只。將屠宰后采集的綿羊卵巢組織放在1.2 mL無酶凍存管內,迅速放入液氮中保存并帶回,以備提取lncRNA。
利用TRIzol法分別提取每個卵巢組織的總RNA,1.2%瓊脂糖凝膠電泳測定RNA完整性及是否存在DNA污染,Nanodorp 2000檢測RNA濃度。將測定合格的樣品送至北京諾禾致源公司進行RNA-seq。
按照要求進行lncRNA文庫構建,庫檢合格后使用Illumina上機測序。高通量測序平臺Illumina測序得到的原始數據(Raw Reads)經過質控并去除低質量、有污染、含有接頭等序列后的數據(Clean Reads)用作后續分析。使用HISAT2[15]對過濾后的reads進行參考基因組的比對分析。使用StringTie[16]對HISAT2比對的結果進行轉錄本的拼接,得到盡可能小的轉錄本集合,并對轉錄本進行定量分析。
基于轉錄組拼接結果,使用Cuffmerge軟件得到合并的轉錄本集合,根據lncRNA的結構特點以及不編碼蛋白的功能特點,設置一系列嚴格的篩選條件,通過5步篩選:外顯子數篩選、轉錄本長度篩選、轉錄本已知注釋篩選、轉錄本表達量篩選、編碼潛能篩選,取CPC2、CNCI、PFAM三種工具預測結果的交集,作為本次分析預測得到的lncRNA數據集進行后續的分析。使用StringTie-eB軟件對包括mRNA、lncRNA以及TUCP在內的轉錄本進行定量分析,得到各樣本轉錄本的FPKM信息。使用cuffdiff對不同類型轉錄本整體進行差異分析,篩選閾值設置為|log2(fold change)|≥2且q-value <0.05。
lncRNA主要通過調控mRNA來實現功能。本研究通過lncRNA與蛋白編碼基因的位置關系(co-location)和表達相關性(co-expression)預測其生物學功能。將co-location的閾值設定為lncRNA上、下游100 kb,后續再通過對lncRNA共定位的mRNA基因進行功能富集分析預測lncRNA的主功能。采用Pearson相關系數法分析樣本間lncRNA與mRNA的相關性,取相關性絕對值大于0.95的mRNA基因進行功能富集分析,預測與mRNA共表達的lncRNA的主要功能。
Gene Ontology(簡稱 GO,www.geneontology.org)是基因功能國際標準分類體系,使用GOseq[17]對差異lncRNA的靶基因進行GO富集分析。KEGG (Kyoto Encyclopedia of Genes and Genomes) 是有關生物學通路的主要公共數據庫[18],使用KOBAS(2.0)[19]進行通路富集分析,FDR≤0.05為顯著富集。
為驗證測序結果的可靠性,隨機選取5個差異表達的lncRNA進行qRT-PCR驗證。使用Primer 5.0設計引物,GAPDH作為內參基因,引物均由北京擎科生物科技有限公司合成,引物序列信息如表1所示。PCR反應體系(20 μL):10 μL lncRNA PreMix,上、下游引物各0.5 μL,1 μL cDNA和8 μL RNase-Free ddH2O。PCR反應條件:95 ℃ 3 min;95 ℃ 5 s,57 ℃ 10 s,72 ℃ 15 s,40個循環。使用2-ΔΔCt方法計算lncRNA的相對表達量,并與轉錄組結果進行對比驗證。

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

表2 數據產出質量情況

表3 Reads與參考基因組比對情況
對OAF和OAL的轉錄本進行分步篩選,鑒定得到候選lncRNAs轉錄本數目為21 085個,其中包括221個已知lncRNAs和20 864個新lncRNAs。進一步對候選lncRNA進行分類,如圖1A所示,基因間型lncRNAs占46.5%,反義型lncRNAs占8.6%,內含子型的lncRNAs占44.9%。
將lncRNA與mRNA進行結構比較分析,得到lncRNA與mRNA的轉錄本長度、外顯子個數、ORF長度(圖1D)的對比情況。結果發現,lncRNA和mRNA的長度整體趨勢基本一致(圖1B);lncRNA的外顯子個數多為2個,顯著低于mRNA(圖1C),這與馬駿杰等[20]的研究結果一致;mRNA的ORF長度明顯大于lncRNA,這提示lncRNA在轉錄及轉錄后等表觀調控方面起著重要作用[21]。
利用cuffdiff進行差異分析得到初步差異分析的結果。設置篩選差異基因的標準為|log2(fold change)|≥2且q-value<0.05,得到的差異基因結果如圖2所示,在OAF和OAL中獲得1 379個差異表達的lncRNAs,其中1 158個表達上調,221個表達下調。隨機挑選5個差異表達lncRNAs進行qRT-PCR,驗證RNA測序數據的可靠性。驗證結果與測序數據結果基本一致(圖3),表明測序數據結果可靠性高。
對1 379個差異表達的lncRNA的共定位靶基因進行生物功能分析,共顯著富集到3個GO條目(P<0.05),包括1個分子功能(molecular function,MF),2個生物學過程(biological process,BP),沒有富集到細胞組分(cellular component,CC)。顯著富集到的GO條目如表4所示,在分子功能歸類中,差異基因顯著富集到谷胱甘肽轉移酶活性(glutathione transferase activity)。在生物學過程中,差異基因顯著富集到病毒過程的負調控(negative regulation of viral process)、病毒基因組復制的負調控(negative regulation of viral genome replication)。

A. lncRNA類型分布; B. lncRNA與mRNA外顯子數目;C. lncRNA與mRNA的轉錄本長度;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的結構特征比較分析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代謝外源物質(metabo-lism of xenobiotics by cytochrome P450)等。P值前20的通路見圖4。
對1 379個差異表達lncRNAs的共表達靶基因進行生物功能分析,共顯著富集到482個GO條目(P<0.05),包括63個分子功能,378個生物學過程,41個細胞組分。前10條顯著富集到的GO條目如表5所示,在分子功能歸類中,差異基因顯著富集到結合(binding)、蛋白結合(protein binding)等。在生物學過程中,差異基因顯著富集到卵泡發育(ovarian follicle development)、排卵周期過程(ovulation cycle process)等。在細胞組分中,差異基因顯著富集到質膜(plasma membrane part)、神經元投射(neuron projection)等。
在KEGG富集分析中,差異表達lncRNA的共表達靶基因涉及275條信號通路,其中共顯著富集到28條通路(P<0.05),包括鈣離子信號通路(calcium signaling pathway)、卵母細胞減數分裂(oocyte meiosis)、cAMP信號通路(cAMP signaling pathway)、催產素信號通路(oxytocin signaling pathway)、MAPK信號通路(MAPK signaling pathway)、甲狀腺激素合成通路(thyroid hormone synthesis)、雌激素信號通路(estrogen signaling pathway)等。P值前20的通路見圖5。
本研究KEGG富集分析得到了可能與繁殖有關的通路基因(表6)。這些通路包括鈣離子信號通路、卵母細胞減數分裂、cAMP信號通路、催產素信號通路、MAPK信號通路、甲狀腺激素合成通路、雌激素信號通路。
通過篩選KGEE富集到信號通路,篩選出5個顯著差異基因,與多個不同lncRNA存在靶標關系(表7)。其中多個lncRNAs靶向同一個差異基因,如LNC_003902、LNC_003906、LNC_003907均靶向PPP3CA等;也存在一個lncRNA同時具有多個靶基因,如ADCY5、ADCY1、CDC25A等都是LNC_012847的潛在目標,這些lncRNAs呈現出了多種方向的調控作用。

圖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 繁殖相關的候選信號通路及基因

表7 差異基因與lncRNAs的靶標關系
繁殖能力對綿羊的經濟效益有重要影響。lncRNA因其在基因調控網絡和廣泛的生物學過程中的作用而受到廣泛關注。研究發現,lncRNA涉及多種動物繁殖相關的過程,如妊娠[22]、性腺激素應答[23]、卵子成熟[24]和胎盤形成[25]等。本研究對綿羊不同發情周期卵巢轉錄組分析,選擇不同發情周期的卵巢篩選出的差異基因顯著富集到的通路及通路中的靶基因進行有關繁殖調控的討論。
本研究在KEGG富集分析得到的多個與繁殖調控相關的候選通路中選擇卵母細胞減數分裂、MAPK、甲狀腺激素合成、雌激素信號通路,并對通路及其中的靶基因進行繁殖相關調控的討論。
卵母細胞減數分裂被認為是形成單倍體配子所需的特殊細胞分裂形式,卵母細胞減數分裂的重新啟動在卵泡期中也起到關鍵作用[26]。卵母細胞減數分裂期間染色體分離的精確調節對哺乳動物的繁殖至關重要。PPP3CA(protein phosphatase 3 catalytic subunit alpha)是蛋白磷酸酶3的催化亞基A,廣泛存在于真核生物的多種細胞中[27]。研究發現,PPP3CA與山羊繁殖力性狀有關[28],基因內缺失突變顯著影響山羊產仔數[29]。在卵泡期顯著上調的LNC_003902、LNC_003906、LNC_003907靶向的PPP3CA在此通路中顯著富集。推測上述lncRNA可能通過調控PPP3CA影響綿羊產羔數。CDC23(cell division cycle 23)是小鼠卵母細胞減數分裂期間紡錘體組裝和染色體分離的關鍵調節因子[30],這表明CDC23在減數分裂卵母細胞中具有獨特的作用。LNC_020278、LNC_020280、LNC_020281、LNC_020279靶向的CDC23在此通路中顯著富集。根據已知CDC23功能的研究,推測上述lncRNA可能在卵母顆粒細胞減數分裂中有重要作用。Liu等[31]在高、低繁殖力寒泊羊篩選的卵泡期差異表達lncRNA也顯著富集到此通路,與本研究結果相似。
MAPK信號通路是是真核生物信號傳遞網絡中的重要途徑之一,MAPK被一系列細胞外的刺激信號激活并介導信號從細胞膜向細胞核傳導,是調節細胞增殖、分化、凋亡和死亡等多種生理過程的關鍵信號通路[32-33]。MAPK信號通路被催乳素抑制,該通路的功能對正常卵泡發育至關重要[34]。MAPK1(mitogen-activated protein kinase 1)是MAP激酶家族的成員,在顆粒細胞中的磷酸化由LH受體激活引起并介導卵母細胞成熟[35-36]。本研究發現,在黃體期顯著上調的LNC_011239靶向的MAPK1在此通路顯著富集,根據靶基因的功能分析,推測LNC_011239可能在黃體期通過調控MAPK1影響卵母細胞成熟。
甲狀腺激素合成通路已被確定對脊椎動物胚胎發生和胎兒成熟至關重要。此外,甲狀腺素(T4)和三碘甲狀腺原氨酸(T3)對正常發育、生長和代謝穩態至關重要[37]。有研究報道,大鼠在懷孕期間甲狀腺激素缺乏,導致產仔數減少[38]。雌激素是由卵巢和胎盤分泌的一種類固醇類激素,在動物發情、分娩等生殖功能方面有重要作用[39],雌激素信號通路是指雌激素與細胞核或細胞膜中的雌激素受體(ER)結合后,作用于雌激素受體反應元件調節基因表達。當雌激素與受體結合后,可以維持雌性個體第二性征,可以刺激和維持黃體功能,在雌性動物繁殖周期中對卵泡的生長發育起重要作用。ADCY1(adenylate cyclase 1)和ADCY5(adenylate cyclase 5)是腺苷酸環化酶家族的成員。ADCY1和ADCY5被發現是影響山羊和奶牛繁殖力的候選基因[40-41]。在卵泡期顯著上調的LNC_003902、LNC_003906、LNC_003907、LNC_012843、LNC_012847、LNC_021474、LNC_021475和LNC_012847靶向的ADCY1和ADCY5在這兩條信號通路均有顯著富集,提示LNC_012847可能與綿羊繁殖力存在一定相關性。
本研究采用高通量測序技術和生物信息學工具,對發情周期不同階段綿羊卵巢的lncRNAs進行了鑒定,在OAF和OAL中獲得1 379個差異表達的lncRNAs,其中1 158個表達上調,221個表達下調。共表達和共定位的GO注釋和KEGG富集分析差異lncRNAs的靶基因富集到多個繁殖相關信號通路中,其中LNC_011239、LNC_012847、LNC_003902、LNC_003906、LNC_003907等靶向的MAPK1、ADCY1、ADCY5、PPP3CA和CDC23可能發揮關鍵的調控作用。此外,這些差異lncRNAs表達譜為揭示綿羊繁殖能力的分子機制提供了有價值的資源。綿羊繁殖相關的候選lncRNAs及鑒定的關鍵靶基因可能需要通過敲除、過表達等方法進一步研究,以驗證lncRNA在綿羊卵巢中的特異功能。