張文英, 高浩然, 潘銳, 常浩雯, 徐樂(lè)
長(zhǎng)江大學(xué)作物抗逆技術(shù)研究中心,長(zhǎng)江大學(xué)農(nóng)學(xué)院,湖北 荊州 434025
大麥(HordeumvulgareL.)是世界上最古老的禾本科作物之一,是全球第四大糧食作物,主要用作糧食、飼料加工、啤酒釀造等[1]。種子萌發(fā)是植物生命的起始階段,決定著幼苗的形成和生長(zhǎng)[2]。大麥種子的萌發(fā)速度和整齊度直接影響啤酒釀造的成本和質(zhì)量[3]。因此,深入了解大麥種子萌發(fā)的分子調(diào)控機(jī)理意義重大。
近些年,以RNA-Seq為代表的轉(zhuǎn)錄組測(cè)序技術(shù)極大推動(dòng)了在RNA水平研究基因表達(dá)及轉(zhuǎn)錄調(diào)控機(jī)制的研究。在大麥種子萌發(fā)方面,mRNA分析揭示了大麥種子萌發(fā)過(guò)程的耐鹽分子機(jī)制[4]、赤霉素合成位置以及合成降解途徑[5],miRNA分析發(fā)現(xiàn)大麥種子萌發(fā)過(guò)程中miR156等5種miRNA在種胚中高度表達(dá)[6]。
長(zhǎng)鏈非編碼RNA(long non-coding RNA, lncRNA)是一類本身不編碼蛋白、轉(zhuǎn)錄本長(zhǎng)度超過(guò)200nt的長(zhǎng)鏈非編碼RNA分子。lncRNA曾被認(rèn)為是“轉(zhuǎn)錄噪音”不具有任何生物學(xué)功能。近些年的研究表明,lncRNA在植物的生長(zhǎng)發(fā)育、表觀遺傳、脅迫響應(yīng)等方面起著非常重要的調(diào)節(jié)功能[7, 8]。
通常情況下,大部分的轉(zhuǎn)錄組研究?jī)H關(guān)注某一類或兩類RNA,沒有關(guān)注多種類型RNA之間的相互作用及其在轉(zhuǎn)錄調(diào)控中的作用機(jī)制。2011年P(guān)ANDOLFI 研究小組提出競(jìng)爭(zhēng)性內(nèi)源性RNA假說(shuō),認(rèn)為mRNA、lncRNA、環(huán)狀RNA(circRNA)和假基因轉(zhuǎn)錄物通過(guò)與miRNA競(jìng)爭(zhēng)性結(jié)合來(lái)調(diào)節(jié)靶基因的穩(wěn)定性或翻譯活性,從而實(shí)現(xiàn)轉(zhuǎn)錄后功能基因調(diào)控。ceRNA機(jī)制將miRNA作為相互調(diào)節(jié)RNA的載體,使編碼和非編碼基因成為整個(gè)轉(zhuǎn)錄組內(nèi)復(fù)雜調(diào)控網(wǎng)絡(luò)的有機(jī)整體[9]。越來(lái)越多的研究表明,ceRNA在動(dòng)植物的生長(zhǎng)發(fā)育過(guò)程中發(fā)揮著重要的調(diào)控作用。在玉米種子發(fā)育研究中發(fā)現(xiàn)了7種新的lncRNA具有潛在的ceRNA功能[10]。利用ceRNA調(diào)控網(wǎng)絡(luò)揭示circRNA在擬南芥葉片發(fā)育和水稻旗葉衰老過(guò)程發(fā)揮著重要作用[11,12]。利用ceRNA調(diào)控網(wǎng)絡(luò)找到6個(gè)參與玉米花藥發(fā)育的靶基因?qū)σ约?個(gè)參與油菜花藥發(fā)育的重要circRNA[13,14]。
目前,關(guān)于大麥種子萌發(fā)過(guò)程中l(wèi)ncRNA、mRNA及miRNA構(gòu)建的ceRNA調(diào)控網(wǎng)絡(luò)尚未見報(bào)道。本研究利用大麥種子萌發(fā)相關(guān)的轉(zhuǎn)錄組數(shù)據(jù),確定了種子萌發(fā)過(guò)程中miRNA、mRNA和lncRNA的差異表達(dá)譜,篩選種子萌發(fā)相關(guān)的特異性模塊進(jìn)行GO富集和KEGG分析,構(gòu)建lncRNA-miRNA-mRNA相互作用的ceRNA調(diào)控網(wǎng)絡(luò),以期為研究作物種子萌發(fā)分子機(jī)制提供新思路。
從BioProject數(shù)據(jù)庫(kù)下載大麥種子萌發(fā)相關(guān)的RNA-Seq原始轉(zhuǎn)錄本數(shù)據(jù)PRJNA577343和PRJNA578897,其中包括3個(gè)品種的36個(gè)萌發(fā)樣本和10個(gè)未萌發(fā)樣本;下載PRJNA389146數(shù)據(jù)集的miRNA的count數(shù)據(jù),其中包括1個(gè)種子發(fā)育時(shí)期和2個(gè)種子萌發(fā)時(shí)期的樣本。本研究將種子發(fā)育時(shí)期視為未萌發(fā)種子。
大麥參考基因組的IBSC_v2,基因組DNA序列以及gtf注釋文件均從Ensembl Plants下載[15]。
采用FASTQC對(duì)RNA原始讀段進(jìn)行質(zhì)量評(píng)估和篩選[16]。利用Hisat 2軟件將有效讀段比對(duì)到大麥的參考基因組,計(jì)算每個(gè)基因的TMM值(Trimmed Mean of M-values)[17]。
采用如下方法鑒定大麥種子萌發(fā)相關(guān)的lncRNA:利用StringTie軟件將每個(gè)文庫(kù)比對(duì)結(jié)果組裝起來(lái),獲得轉(zhuǎn)錄本結(jié)構(gòu),利用gffcompare軟件將所有文庫(kù)組裝結(jié)果合并,剔除與蛋白編碼基因重疊的轉(zhuǎn)錄本;剔除序列長(zhǎng)度小于200nt或含有大于300nt開放閱讀框的轉(zhuǎn)錄本;將剩余的轉(zhuǎn)錄本在Pfam數(shù)據(jù)庫(kù)中進(jìn)行比對(duì),剔除包含保守結(jié)構(gòu)域或者模體的轉(zhuǎn)錄本;將剩余的轉(zhuǎn)錄本在NCBI-NR數(shù)據(jù)庫(kù)進(jìn)行BLASTx比對(duì),剔除與數(shù)據(jù)庫(kù)中蛋白序列同源的轉(zhuǎn)錄本;采用CPC2軟件評(píng)估剩余轉(zhuǎn)錄本編碼潛能,剔除具有較高編碼潛能的轉(zhuǎn)錄本;利用Rfam數(shù)據(jù)庫(kù)剔除rRNA、tRNA、snoRNA和snRNA等持家非編碼RNA[18]。
使用DESeq2軟件包對(duì)mRNA和lncRNA數(shù)據(jù)進(jìn)行整理,篩選并確定差異表達(dá)的mRNA和lncRNA;利用edgeR函數(shù)包對(duì)miRNA整理分析,篩選并確定差異表達(dá)的miRNA。計(jì)算3種不同RNA差異表達(dá)倍數(shù),篩選P<0.05且差異表達(dá)倍數(shù)不低于2倍的轉(zhuǎn)錄本為差異表達(dá)轉(zhuǎn)錄本。
利用R軟件WGCNA包分別對(duì)差異表達(dá)的mRNA和lncRNA構(gòu)建基因共表達(dá)網(wǎng)絡(luò)分析。根據(jù)scalefree拓?fù)錁?biāo)準(zhǔn),確定合適的β值作為構(gòu)建網(wǎng)絡(luò)的軟閾值,通過(guò)計(jì)算所有基因之間的皮爾遜相關(guān)性來(lái)構(gòu)建鄰接矩陣,并轉(zhuǎn)化為拓?fù)渚仃?,并根?jù)相似性對(duì)基因進(jìn)行層次聚類,挑選連接度高于0.6的特定模塊進(jìn)行分析[19]。
利用agriGOv2(http://systemsbiology.cau.edu.cn/agriGOv2/)工具進(jìn)行GO富集分析,使用基迪奧在線網(wǎng)站平臺(tái)(https://www.omicshare.com/tools/home/report/koenrich.html)進(jìn)行KEGG通路分析,根據(jù)篩選條件(P<0.05)對(duì)分析結(jié)果進(jìn)行匯總,結(jié)果以氣泡圖呈現(xiàn)。
利用psRNATarget在線工具對(duì)差異表達(dá)的miRNA和篩選出的mRNA及l(fā)ncRNA進(jìn)行靶基因預(yù)測(cè),獲得lncRNA-miRNA-mRNA的ceRNA網(wǎng)絡(luò)數(shù)據(jù),利用Cytescape V3.6.1軟件進(jìn)行可視化分析,根據(jù)“高表達(dá)lncRNA-低表達(dá)miRNA-高表達(dá)mRNA”和“低表達(dá)lncRNA-高表達(dá)miRNA-低表達(dá)mRNA”這2種調(diào)控軸構(gòu)建ceRNA調(diào)控網(wǎng)絡(luò)。
采用qRT-PCR分析lncRNA、miRNA和mRNA的相對(duì)表達(dá)水平,利用Trizol法[20]提取大麥Barke干種子和萌發(fā)24h的RNA樣本,用大麥的actin基因作為lncRNA和mRNA的內(nèi)參基因,以U6作為miRNA的內(nèi)參基因,使用莖環(huán)法對(duì)miRNA進(jìn)行反轉(zhuǎn)錄,引物序列見表1。
表1 熒光定量PCR引物信息Table 1 Primers for fluorescent quantitative PCR
對(duì)3個(gè)參試品種種子萌發(fā)的RNA-seq數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,發(fā)現(xiàn)46個(gè)文庫(kù)總測(cè)序序列數(shù)量為11.4×108,質(zhì)控處理后有效序列總數(shù)為11.2×108,其中有992426882個(gè)讀數(shù)能比對(duì)到大麥參考基因組上。
為了識(shí)別大麥種子萌發(fā)的lncRNA,按照?qǐng)D1所示流程,使用Stringtie軟件進(jìn)行轉(zhuǎn)錄本重組,通過(guò)剔除已知的編碼基因,獲得了35245條候選lncRNA轉(zhuǎn)錄本,最后通過(guò)區(qū)分lncRNA和蛋白編碼轉(zhuǎn)錄本,獲得15053條可靠的lncRNA轉(zhuǎn)錄本。
圖1 大麥種子萌發(fā) lncRNA鑒定流程 Fig.1 Identification processes of lncRNA inbarley seed germination
根據(jù)設(shè)定的篩選條件(P<0.05;|log2Foldchange|>2),從PRJNA577343供試大麥品種中篩選到17369個(gè)差異表達(dá)的mRNA和3523個(gè)差異表達(dá)的lncRNA,從PRJNA578897的品種IK621篩選到8976個(gè)差異表達(dá)的mRNA和1317個(gè)差異表達(dá)的lncRNA,品種IK573篩選到10307個(gè)差異表達(dá)的mRNA和1172個(gè)差異表達(dá)的lncRNA,其中3個(gè)品種共同的差異表達(dá)mRNA有6447個(gè),共同的差異表達(dá)lncRNA有661個(gè)(見圖2)。同時(shí),通過(guò)miRNA Count數(shù)據(jù),鑒定到310個(gè)差異表達(dá)的miRNA。
圖2 不同品種大麥種子萌發(fā)差異表達(dá)的mRNA和lncRNA比較Fig.2 Comparison of differentially expressed mRNA and lncRNA in seed germination of different barley varieties
利用WGCNA對(duì)3個(gè)品種共有的差異表達(dá)的mRNA和lncRNA分別進(jìn)行共表達(dá)分析,基于模塊與種子萌發(fā)的相關(guān)性(R2>0.6)篩選特異性共表達(dá)模塊。
在差異表達(dá)的mRNA分析中,設(shè)置權(quán)重值β為24,結(jié)果如圖3所示。不同模塊用不同顏色表示,共獲得10個(gè)模塊(Grey 模塊表示未分配到任何模塊的基因)。其中,brown模塊包含的基因個(gè)數(shù)最多,有2246個(gè);grey模塊包含的基因個(gè)數(shù)最少,有8個(gè),通過(guò)共表達(dá)網(wǎng)絡(luò)與種子萌發(fā)進(jìn)行關(guān)聯(lián),鑒定到與種子萌發(fā)極顯著正相關(guān)的模塊greenyellow(R2=0.68)(見圖4),該模塊包括1433個(gè)mRNA。
圖3 mRNA聚類樹及模塊構(gòu)建 圖4 mRNA模塊及其與種子萌發(fā)相關(guān)性 Fig.3 mRNA clustering tree and module construction Fig.4 mRNA module and its correlation with seed germination
在差異表達(dá)的lncRNA分析中,設(shè)置權(quán)重值β為12,結(jié)果如圖5所示。通過(guò)lncRNA獲得6個(gè)模塊。其中,blue模塊包含的基因個(gè)數(shù)最多,有325個(gè);grey模塊包含的基因個(gè)數(shù)最少,有10個(gè),通過(guò)共表達(dá)網(wǎng)絡(luò)與種子萌發(fā)進(jìn)行關(guān)聯(lián)分析,鑒定到與種子萌發(fā)極顯著正相關(guān)的模塊brown(R2=0.61)(見圖6),該模塊含有116個(gè)lncRNA。
圖5 lncRNA聚類樹及模塊構(gòu)建 圖6 lncRNA模塊及其與種子萌發(fā)相關(guān)性 Fig.5 lncRNA clustering tree and module construction Fig.6 lncRNA module and its correlation with seed germination
為進(jìn)一步解析greenyellow模塊的生物學(xué)功能,利用agriGO工具對(duì)模塊內(nèi)所有基因進(jìn)行GO富集分析,這些GO通路(見圖7(a))涉及到生物學(xué)過(guò)程、分子功能以及細(xì)胞組分。分析結(jié)果表明,該模塊共富集到26個(gè)生物過(guò)程,其主要涉及水分響應(yīng)(GO:0009415)、對(duì)酸性化學(xué)物質(zhì)的響應(yīng)(GO:0001101)、對(duì)無(wú)機(jī)物的響應(yīng)(GO:0010035)、對(duì)含氧化合物的反應(yīng)(GO:1901700)、胚發(fā)育(GO:0009790)等通路以及一些酶活通路,如肽酶調(diào)節(jié)活性(GO:0030414)、肽酶抑制劑活性(GO:0004866)、內(nèi)肽酶調(diào)節(jié)活性(GO:0061135)、內(nèi)肽酶抑制劑活性(GO:0061134)等。KEGG富集分析表明,greenyellow模塊共有67個(gè)基因富集到7條通路,主要涉及MAPK信號(hào)轉(zhuǎn)導(dǎo)通路、植物激素信號(hào)轉(zhuǎn)導(dǎo)通路等(見圖7(b))。
圖7 greenyellow模塊的GO和KEGG富集分析Fig.7 GO and KEGG enrichment analysis of greenyellow module
利用上述篩選到的mRNA、lncRNA與miRNA進(jìn)行靶基因預(yù)測(cè),獲得168個(gè)miRNA-mRNA靶基因?qū)σ约?10個(gè)miRNA-lncRNA靶基因?qū)?,最終根據(jù)“高表達(dá)lncRNA-低表達(dá)miRNA-高表達(dá)mRNA”和“低表達(dá)lncRNA-高表達(dá)miRNA-低表達(dá)mRNA”這2種調(diào)控軸構(gòu)建ceRNA調(diào)控網(wǎng)絡(luò)(見圖8),該網(wǎng)絡(luò)包含38個(gè)lncRNA、16個(gè)miRNA以及18個(gè)mRNA。其中“高表達(dá)lncRNA-低表達(dá)miRNA-高表達(dá)mRNA”的子網(wǎng)絡(luò)包含了4個(gè)lncRNA、2個(gè)miRNA以及1個(gè)mRNA,顯示基因HORVU0Hr1G039470的表達(dá)可能受到2個(gè)miRNA和4個(gè)lncRNA調(diào)控?!暗捅磉_(dá)lncRNA-高表達(dá)miRNA-低表達(dá)mRNA”子網(wǎng)絡(luò)包含34個(gè)lncRNA、14個(gè)miRNA以及17個(gè)mRNA,其中miR854的表達(dá)受到19個(gè)靶向基因的調(diào)控,包括3個(gè)mRNA和16個(gè)lncRNA;而miR5054、miR5059、miR5060、miR5304、miR6432以及miR916的靶向lncRNA和mRNA均只有1個(gè)。此外,對(duì)ceRNA調(diào)控網(wǎng)絡(luò)中的mRNA在大麥中進(jìn)行了注釋和在擬南芥中進(jìn)行了同源基因的注釋(見表2)。
表2 ceRNA調(diào)控網(wǎng)絡(luò)中mRNA功能注釋Table 2 mRNA function annotations in the ceRNA regulatory network
圖8 大麥種子萌發(fā)ceRNA網(wǎng)絡(luò)調(diào)控圖Fig.8 ceRNA network regulatory map of barley seed germination
利用qRT-PCR檢測(cè)基因HORVU0Hr1G039470、HORVU6Hr1G031480、MSTRG.23227、MSTRG.21252、miR1858以及miR2611的相對(duì)表達(dá)水平(見圖9)。結(jié)果顯示,HORVU0Hr1G039470、MSTRG.21252以及miR1858上調(diào),基因HORVU6Hr1G031480、MSTRG.23227以及miR2611下調(diào),這些基因的上下調(diào)模式與預(yù)測(cè)的ceRNA調(diào)控網(wǎng)絡(luò)中的上下調(diào)模式一致。由此,初步驗(yàn)證了ceRNA調(diào)控網(wǎng)絡(luò)的可行性,同時(shí)初步驗(yàn)證了HORVU0Hr1G039470-miR2611-MSTRG.21252和HORVU6Hr1G031480-miR1858-MSTRG.23227這2個(gè)關(guān)系對(duì)具有相應(yīng)的ceRNA功能。
圖9 大麥種子萌發(fā)過(guò)程中不同RNA的表達(dá)水平Fig.9 Expression levels of different RNAs during barley seed germination
種子萌發(fā)是一個(gè)復(fù)雜的生物學(xué)過(guò)程,對(duì)植物的生存和繁衍至關(guān)重要。本研究通過(guò)對(duì)大麥萌發(fā)轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行分析,鑒定到了15053個(gè)lncRNA。通過(guò)差異表達(dá)基因分析,得到6447個(gè)差異表達(dá)的mRNA以及661個(gè)差異表達(dá)的lncRNA以及差異表達(dá)310個(gè)miRNA。對(duì)差異表達(dá)的mRNA和lncRNA進(jìn)行WGCNA分析,分別鑒定到1個(gè)與大麥種子萌發(fā)高度相關(guān)的共表達(dá)模塊,在這2個(gè)模塊中,分別包含了1433個(gè)mRNA和116個(gè)lncRNA。高度相關(guān)的共表達(dá)模塊中有67個(gè)mRNA被富集到MAPK信號(hào)通路,植物激素信號(hào)轉(zhuǎn)導(dǎo)通路等KEGG通路上,這些通路均已證明與種子萌發(fā)有關(guān)[4]。通過(guò)對(duì)篩選到的基因進(jìn)行靶基因預(yù)測(cè),成功構(gòu)建了一個(gè)包含2個(gè)子網(wǎng)絡(luò)的ceRNA網(wǎng)絡(luò)圖,包含16個(gè)miRNA節(jié)點(diǎn),38個(gè)lncRNA節(jié)點(diǎn)以及18個(gè)mRNA節(jié)點(diǎn)。
在“高表達(dá)lncRNA-低表達(dá)miRNA-高表達(dá)mRNA”子網(wǎng)絡(luò)中,作為糖基轉(zhuǎn)移酶基因HORVU0Hr1G039470的表達(dá)可能受到miR5052、miR1858以及4個(gè)lncRNA的調(diào)控。前人研究表明,糖基轉(zhuǎn)移酶與種子的萌發(fā)有關(guān)[21],miR1858與籽粒發(fā)育有關(guān)[22],籽粒良好的發(fā)育是種子萌發(fā)的關(guān)鍵。MSTRG.2322或MSTRG.10816或MSTRG.10592可能通過(guò)競(jìng)爭(zhēng)性地結(jié)合miR1858來(lái)調(diào)節(jié)HORVU0Hr1G039470的表達(dá),從而調(diào)節(jié)大麥種子萌發(fā)過(guò)程。在“低表達(dá)lncRNA-高表達(dá)miRNA-低表達(dá)mRNA”子網(wǎng)絡(luò)中,miR5054、miR5671及miR854與調(diào)控MAPK信號(hào)通路或植物激素信號(hào)轉(zhuǎn)導(dǎo)通路相關(guān)[23-25],miR946參與種子的萌發(fā)以及休眠[26]。MSTRG.14519與HORVU3Hr1G089250可能競(jìng)爭(zhēng)性地結(jié)合miR5054,從而影響種子萌發(fā)過(guò)程中激素的表達(dá);MSTRG.33901或MSTRG.5241或MSTRG.15791能夠與miR5671進(jìn)行結(jié)合,從而調(diào)控谷氨酸代謝相關(guān)的基因HORVU4Hr1G007610、HORVU3Hr1G003050;與種子萌發(fā)相關(guān)的miR946以及HORVU2Hr1G110230通過(guò)與MSTRG.42080或MSTRG.11925或MSTRG.11079或MSTRG.525行使ceRNA功能,從而對(duì)大麥種子萌發(fā)產(chǎn)生影響。此外,miR854受到的靶向調(diào)控基因最多,說(shuō)明該基因可能是大麥種子萌發(fā)過(guò)程中一個(gè)極為重要的節(jié)點(diǎn)miRNA。
本研究通過(guò)qRT-PCR驗(yàn)證了ceRNA調(diào)控網(wǎng)絡(luò)中的2個(gè)關(guān)系對(duì),其中HORVU0Hr1G039470和MSTRG.21252上調(diào),而miR2611在大麥種子萌發(fā)中下調(diào),基因HORVU6Hr1G031480、MSTRG.23227下調(diào),而miR1858在大麥種子萌發(fā)中上調(diào)。這表明MSTRG.21252上調(diào)和miR2611下調(diào)可能促進(jìn)miR2611對(duì)靶基因HORVU0Hr1G039470的促進(jìn)作用,而MSTRG.23227下調(diào)和miR1858上調(diào)可能促進(jìn)miR1858對(duì)靶基因HORVU6Hr1G031480的抑制作用。
綜上所述,本研究鑒定出大麥種子萌發(fā)相關(guān)的差異表達(dá)基因,構(gòu)建了大麥種子萌發(fā)相關(guān)的lncRNA-miRNA-mRNA互作ceRNA調(diào)控網(wǎng)絡(luò),初步揭示了3種RNA在大麥種子萌發(fā)中的調(diào)控模式。這些結(jié)果將有助于揭示種子萌發(fā)的分子調(diào)控機(jī)制,為種子萌發(fā)lncRNA、mRNA以及miRNA功能的進(jìn)一步研究奠定基礎(chǔ)。