大麥種子萌發(fā)lncRNA-miRNA-mRNA差異表達(dá)分析及ceRNA調(diào)控網(wǎng)絡(luò)構(gòu)建

張文英, 高浩然, 潘銳, 常浩雯, 徐樂(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ī)制提供新思路。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

從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]。

1.2 lncRNA鑒定

采用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]。

1.3 差異表達(dá)基因篩選

使用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)錄本。

1.4 WGCNA分析

利用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]。

1.5 富集分析及代謝通路分析

利用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)。

1.6 ceRNA網(wǎng)絡(luò)構(gòu)建

利用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ò)。

1.7 實(shí)時(shí)定量聚合酶鏈反應(yīng)(qRT-PCR)分析

采用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

2 結(jié)果與分析

2.1 lncRNA鑒定結(jié)果

對(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

2.2 差異表達(dá)基因篩選

根據(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

2.3 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建及特異性模塊鑒定

利用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

2.4 mRNA特異性模塊富集分析和通路分析

為進(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

2.5 ceRNA網(wǎng)絡(luò)構(gòu)建結(jié)果

利用上述篩選到的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

2.6 qRT-PCR分析結(jié)果

利用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

3 討論與結(jié)論

種子萌發(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ǔ)。