大麥種子萌發lncRNA-miRNA-mRNA差異表達分析及ceRNA調控網絡構建

張文英, 高浩然, 潘銳, 常浩雯, 徐樂

長江大學作物抗逆技術研究中心，長江大學農學院，湖北 荊州 434025

大麥(HordeumvulgareL.)是世界上最古老的禾本科作物之一，是全球第四大糧食作物，主要用作糧食、飼料加工、啤酒釀造等[1]。種子萌發是植物生命的起始階段，決定著幼苗的形成和生長[2]。大麥種子的萌發速度和整齊度直接影響啤酒釀造的成本和質量[3]。因此，深入了解大麥種子萌發的分子調控機理意義重大。

近些年，以RNA-Seq為代表的轉錄組測序技術極大推動了在RNA水平研究基因表達及轉錄調控機制的研究。在大麥種子萌發方面，mRNA分析揭示了大麥種子萌發過程的耐鹽分子機制[4]、赤霉素合成位置以及合成降解途徑[5]，miRNA分析發現大麥種子萌發過程中miR156等5種miRNA在種胚中高度表達[6]。

長鏈非編碼RNA(long non-coding RNA, lncRNA)是一類本身不編碼蛋白、轉錄本長度超過200nt的長鏈非編碼RNA分子。lncRNA曾被認為是“轉錄噪音”不具有任何生物學功能。近些年的研究表明，lncRNA在植物的生長發育、表觀遺傳、脅迫響應等方面起著非常重要的調節功能[7, 8]。

通常情況下，大部分的轉錄組研究僅關注某一類或兩類RNA，沒有關注多種類型RNA之間的相互作用及其在轉錄調控中的作用機制。2011年PANDOLFI 研究小組提出競爭性內源性RNA假說，認為mRNA、lncRNA、環狀RNA(circRNA)和假基因轉錄物通過與miRNA競爭性結合來調節靶基因的穩定性或翻譯活性，從而實現轉錄后功能基因調控。ceRNA機制將miRNA作為相互調節RNA的載體，使編碼和非編碼基因成為整個轉錄組內復雜調控網絡的有機整體[9]。越來越多的研究表明，ceRNA在動植物的生長發育過程中發揮著重要的調控作用。在玉米種子發育研究中發現了7種新的lncRNA具有潛在的ceRNA功能[10]。利用ceRNA調控網絡揭示circRNA在擬南芥葉片發育和水稻旗葉衰老過程發揮著重要作用[11,12]。利用ceRNA調控網絡找到6個參與玉米花藥發育的靶基因對以及2個參與油菜花藥發育的重要circRNA[13,14]。

目前，關于大麥種子萌發過程中lncRNA、mRNA及miRNA構建的ceRNA調控網絡尚未見報道。本研究利用大麥種子萌發相關的轉錄組數據，確定了種子萌發過程中miRNA、mRNA和lncRNA的差異表達譜，篩選種子萌發相關的特異性模塊進行GO富集和KEGG分析，構建lncRNA-miRNA-mRNA相互作用的ceRNA調控網絡，以期為研究作物種子萌發分子機制提供新思路。

1 材料與方法

1.1 數據來源

從BioProject數據庫下載大麥種子萌發相關的RNA-Seq原始轉錄本數據PRJNA577343和PRJNA578897，其中包括3個品種的36個萌發樣本和10個未萌發樣本；下載PRJNA389146數據集的miRNA的count數據，其中包括1個種子發育時期和2個種子萌發時期的樣本。本研究將種子發育時期視為未萌發種子。

大麥參考基因組的IBSC_v2，基因組DNA序列以及gtf注釋文件均從Ensembl Plants下載[15]。

1.2 lncRNA鑒定

采用FASTQC對RNA原始讀段進行質量評估和篩選[16]。利用Hisat 2軟件將有效讀段比對到大麥的參考基因組，計算每個基因的TMM值(Trimmed Mean of M-values)[17]。

采用如下方法鑒定大麥種子萌發相關的lncRNA：利用StringTie軟件將每個文庫比對結果組裝起來，獲得轉錄本結構，利用gffcompare軟件將所有文庫組裝結果合并，剔除與蛋白編碼基因重疊的轉錄本；剔除序列長度小于200nt或含有大于300nt開放閱讀框的轉錄本；將剩余的轉錄本在Pfam數據庫中進行比對，剔除包含保守結構域或者模體的轉錄本；將剩余的轉錄本在NCBI-NR數據庫進行BLASTx比對，剔除與數據庫中蛋白序列同源的轉錄本；采用CPC2軟件評估剩余轉錄本編碼潛能，剔除具有較高編碼潛能的轉錄本；利用Rfam數據庫剔除rRNA、tRNA、snoRNA和snRNA等持家非編碼RNA[18]。

1.3 差異表達基因篩選

使用DESeq2軟件包對mRNA和lncRNA數據進行整理，篩選并確定差異表達的mRNA和lncRNA；利用edgeR函數包對miRNA整理分析，篩選并確定差異表達的miRNA。計算3種不同RNA差異表達倍數，篩選P<0.05且差異表達倍數不低于2倍的轉錄本為差異表達轉錄本。

1.4 WGCNA分析

利用R軟件WGCNA包分別對差異表達的mRNA和lncRNA構建基因共表達網絡分析。根據scalefree拓撲標準，確定合適的β值作為構建網絡的軟閾值，通過計算所有基因之間的皮爾遜相關性來構建鄰接矩陣，并轉化為拓撲矩陣，并根據相似性對基因進行層次聚類，挑選連接度高于0.6的特定模塊進行分析[19]。

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

利用agriGOv2(http://systemsbiology.cau.edu.cn/agriGOv2/)工具進行GO富集分析，使用基迪奧在線網站平臺(https://www.omicshare.com/tools/home/report/koenrich.html)進行KEGG通路分析，根據篩選條件(P<0.05)對分析結果進行匯總，結果以氣泡圖呈現。

1.6 ceRNA網絡構建

利用psRNATarget在線工具對差異表達的miRNA和篩選出的mRNA及lncRNA進行靶基因預測，獲得lncRNA-miRNA-mRNA的ceRNA網絡數據，利用Cytescape V3.6.1軟件進行可視化分析，根據“高表達lncRNA-低表達miRNA-高表達mRNA”和“低表達lncRNA-高表達miRNA-低表達mRNA”這2種調控軸構建ceRNA調控網絡。

1.7 實時定量聚合酶鏈反應(qRT-PCR)分析

采用qRT-PCR分析lncRNA、miRNA和mRNA的相對表達水平，利用Trizol法[20]提取大麥Barke干種子和萌發24h的RNA樣本，用大麥的actin基因作為lncRNA和mRNA的內參基因，以U6作為miRNA的內參基因，使用莖環法對miRNA進行反轉錄，引物序列見表1。

表1 熒光定量PCR引物信息Table 1 Primers for fluorescent quantitative PCR

2 結果與分析

2.1 lncRNA鑒定結果

對3個參試品種種子萌發的RNA-seq數據進行統計分析，發現46個文庫總測序序列數量為11.4×108，質控處理后有效序列總數為11.2×108，其中有992426882個讀數能比對到大麥參考基因組上。

為了識別大麥種子萌發的lncRNA，按照圖1所示流程，使用Stringtie軟件進行轉錄本重組，通過剔除已知的編碼基因，獲得了35245條候選lncRNA轉錄本，最后通過區分lncRNA和蛋白編碼轉錄本，獲得15053條可靠的lncRNA轉錄本。

圖1 大麥種子萌發 lncRNA鑒定流程 Fig.1 Identification processes of lncRNA inbarley seed germination

2.2 差異表達基因篩選

根據設定的篩選條件(P<0.05；|log2Foldchange|>2)，從PRJNA577343供試大麥品種中篩選到17369個差異表達的mRNA和3523個差異表達的lncRNA，從PRJNA578897的品種IK621篩選到8976個差異表達的mRNA和1317個差異表達的lncRNA，品種IK573篩選到10307個差異表達的mRNA和1172個差異表達的lncRNA，其中3個品種共同的差異表達mRNA有6447個，共同的差異表達lncRNA有661個(見圖2)。同時，通過miRNA Count數據，鑒定到310個差異表達的miRNA。

圖2 不同品種大麥種子萌發差異表達的mRNA和lncRNA比較Fig.2 Comparison of differentially expressed mRNA and lncRNA in seed germination of different barley varieties

2.3 加權基因共表達網絡構建及特異性模塊鑒定

利用WGCNA對3個品種共有的差異表達的mRNA和lncRNA分別進行共表達分析，基于模塊與種子萌發的相關性(R2>0.6)篩選特異性共表達模塊。

在差異表達的mRNA分析中，設置權重值β為24，結果如圖3所示。不同模塊用不同顏色表示，共獲得10個模塊(Grey 模塊表示未分配到任何模塊的基因)。其中，brown模塊包含的基因個數最多，有2246個；grey模塊包含的基因個數最少，有8個，通過共表達網絡與種子萌發進行關聯，鑒定到與種子萌發極顯著正相關的模塊greenyellow(R2=0.68)(見圖4)，該模塊包括1433個mRNA。

圖3 mRNA聚類樹及模塊構建 圖4 mRNA模塊及其與種子萌發相關性 Fig.3 mRNA clustering tree and module construction Fig.4 mRNA module and its correlation with seed germination

在差異表達的lncRNA分析中，設置權重值β為12，結果如圖5所示。通過lncRNA獲得6個模塊。其中，blue模塊包含的基因個數最多，有325個；grey模塊包含的基因個數最少，有10個，通過共表達網絡與種子萌發進行關聯分析，鑒定到與種子萌發極顯著正相關的模塊brown(R2=0.61)(見圖6)，該模塊含有116個lncRNA。

圖5 lncRNA聚類樹及模塊構建 圖6 lncRNA模塊及其與種子萌發相關性 Fig.5 lncRNA clustering tree and module construction Fig.6 lncRNA module and its correlation with seed germination

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

為進一步解析greenyellow模塊的生物學功能，利用agriGO工具對模塊內所有基因進行GO富集分析，這些GO通路(見圖7(a))涉及到生物學過程、分子功能以及細胞組分。分析結果表明，該模塊共富集到26個生物過程，其主要涉及水分響應(GO:0009415)、對酸性化學物質的響應(GO:0001101)、對無機物的響應(GO:0010035)、對含氧化合物的反應(GO:1901700)、胚發育(GO:0009790)等通路以及一些酶活通路，如肽酶調節活性(GO:0030414)、肽酶抑制劑活性(GO:0004866)、內肽酶調節活性(GO:0061135)、內肽酶抑制劑活性(GO:0061134)等。KEGG富集分析表明，greenyellow模塊共有67個基因富集到7條通路，主要涉及MAPK信號轉導通路、植物激素信號轉導通路等(見圖7(b))。

圖7 greenyellow模塊的GO和KEGG富集分析Fig.7 GO and KEGG enrichment analysis of greenyellow module

2.5 ceRNA網絡構建結果

利用上述篩選到的mRNA、lncRNA與miRNA進行靶基因預測，獲得168個miRNA-mRNA靶基因對以及310個miRNA-lncRNA靶基因對，最終根據“高表達lncRNA-低表達miRNA-高表達mRNA”和“低表達lncRNA-高表達miRNA-低表達mRNA”這2種調控軸構建ceRNA調控網絡(見圖8)，該網絡包含38個lncRNA、16個miRNA以及18個mRNA。其中“高表達lncRNA-低表達miRNA-高表達mRNA”的子網絡包含了4個lncRNA、2個miRNA以及1個mRNA，顯示基因HORVU0Hr1G039470的表達可能受到2個miRNA和4個lncRNA調控。“低表達lncRNA-高表達miRNA-低表達mRNA”子網絡包含34個lncRNA、14個miRNA以及17個mRNA，其中miR854的表達受到19個靶向基因的調控，包括3個mRNA和16個lncRNA；而miR5054、miR5059、miR5060、miR5304、miR6432以及miR916的靶向lncRNA和mRNA均只有1個。此外，對ceRNA調控網絡中的mRNA在大麥中進行了注釋和在擬南芥中進行了同源基因的注釋(見表2)。

表2 ceRNA調控網絡中mRNA功能注釋Table 2 mRNA function annotations in the ceRNA regulatory network

圖8 大麥種子萌發ceRNA網絡調控圖Fig.8 ceRNA network regulatory map of barley seed germination

2.6 qRT-PCR分析結果

利用qRT-PCR檢測基因HORVU0Hr1G039470、HORVU6Hr1G031480、MSTRG.23227、MSTRG.21252、miR1858以及miR2611的相對表達水平(見圖9)。結果顯示，HORVU0Hr1G039470、MSTRG.21252以及miR1858上調，基因HORVU6Hr1G031480、MSTRG.23227以及miR2611下調，這些基因的上下調模式與預測的ceRNA調控網絡中的上下調模式一致。由此，初步驗證了ceRNA調控網絡的可行性，同時初步驗證了HORVU0Hr1G039470-miR2611-MSTRG.21252和HORVU6Hr1G031480-miR1858-MSTRG.23227這2個關系對具有相應的ceRNA功能。

圖9 大麥種子萌發過程中不同RNA的表達水平Fig.9 Expression levels of different RNAs during barley seed germination

3 討論與結論

種子萌發是一個復雜的生物學過程，對植物的生存和繁衍至關重要。本研究通過對大麥萌發轉錄組數據進行分析，鑒定到了15053個lncRNA。通過差異表達基因分析，得到6447個差異表達的mRNA以及661個差異表達的lncRNA以及差異表達310個miRNA。對差異表達的mRNA和lncRNA進行WGCNA分析，分別鑒定到1個與大麥種子萌發高度相關的共表達模塊，在這2個模塊中，分別包含了1433個mRNA和116個lncRNA。高度相關的共表達模塊中有67個mRNA被富集到MAPK信號通路，植物激素信號轉導通路等KEGG通路上，這些通路均已證明與種子萌發有關[4]。通過對篩選到的基因進行靶基因預測，成功構建了一個包含2個子網絡的ceRNA網絡圖，包含16個miRNA節點，38個lncRNA節點以及18個mRNA節點。

在“高表達lncRNA-低表達miRNA-高表達mRNA”子網絡中，作為糖基轉移酶基因HORVU0Hr1G039470的表達可能受到miR5052、miR1858以及4個lncRNA的調控。前人研究表明，糖基轉移酶與種子的萌發有關[21]，miR1858與籽粒發育有關[22]，籽粒良好的發育是種子萌發的關鍵。MSTRG.2322或MSTRG.10816或MSTRG.10592可能通過競爭性地結合miR1858來調節HORVU0Hr1G039470的表達，從而調節大麥種子萌發過程。在“低表達lncRNA-高表達miRNA-低表達mRNA”子網絡中，miR5054、miR5671及miR854與調控MAPK信號通路或植物激素信號轉導通路相關[23-25]，miR946參與種子的萌發以及休眠[26]。MSTRG.14519與HORVU3Hr1G089250可能競爭性地結合miR5054，從而影響種子萌發過程中激素的表達；MSTRG.33901或MSTRG.5241或MSTRG.15791能夠與miR5671進行結合，從而調控谷氨酸代謝相關的基因HORVU4Hr1G007610、HORVU3Hr1G003050；與種子萌發相關的miR946以及HORVU2Hr1G110230通過與MSTRG.42080或MSTRG.11925或MSTRG.11079或MSTRG.525行使ceRNA功能，從而對大麥種子萌發產生影響。此外，miR854受到的靶向調控基因最多，說明該基因可能是大麥種子萌發過程中一個極為重要的節點miRNA。

本研究通過qRT-PCR驗證了ceRNA調控網絡中的2個關系對，其中HORVU0Hr1G039470和MSTRG.21252上調，而miR2611在大麥種子萌發中下調，基因HORVU6Hr1G031480、MSTRG.23227下調，而miR1858在大麥種子萌發中上調。這表明MSTRG.21252上調和miR2611下調可能促進miR2611對靶基因HORVU0Hr1G039470的促進作用，而MSTRG.23227下調和miR1858上調可能促進miR1858對靶基因HORVU6Hr1G031480的抑制作用。

綜上所述，本研究鑒定出大麥種子萌發相關的差異表達基因，構建了大麥種子萌發相關的lncRNA-miRNA-mRNA互作ceRNA調控網絡，初步揭示了3種RNA在大麥種子萌發中的調控模式。這些結果將有助于揭示種子萌發的分子調控機制，為種子萌發lncRNA、mRNA以及miRNA功能的進一步研究奠定基礎。