microRNA調控膜系統及離子轉運介導青稞鹽堿脅迫的應答

王玉林，原紅軍，扎桑，楊春葆，徐齊君

（省部共建青稞和牦牛種質資源與遺傳改良國家重點實驗室/西藏自治區農牧科學院農業研究所，西藏拉薩850032）

青稞（Hordeum vulgare var.nudum）是我國西部青藏高原常年種植的主要糧食作物和重要牲畜飼料來源[1]，約占該地區糧食播種面積的80%左右。但是，由于對青稞食品需求量的增加以及氣候變化的嚴重影響，青稞產量的提高對滿足青藏高原地區糧食消費需求有著重要的意義[2]。青稞生長在極端環境，因此對非生物脅迫有著獨特的適應機理。

土壤鹽堿化是一種嚴重的非生物脅迫，影響作物生長并造成世界范圍內作物產量損失[3]。據估計，在未來幾十年中，鹽堿脅迫可能會影響越來越多的耕地[4]，并使植物遭到滲透脅迫、高pH脅迫和離子毒性。為了應對高鹽堿脅迫，植物通過復雜的基因調控網絡作出反應，改變基因表達水平并參與膜運輸、信號轉導和能量代謝相關的翻譯后修飾[5-6]。但是，目前對植物鹽堿反應的分子機制了解甚少，提高農作物的鹽堿耐受性是一個重要的全球性問題。

MicroRNA（miRNA）屬于一類小的單鏈RNA，它們在轉錄后水平上負調控基因表達[7]。許多研究表明，植物miRNA在應對非生物脅迫中起重要作用[8-9]，如鹽堿脅迫[10-12]。近年來，miRNA在模式植物和幾種谷物中介導與脅迫適應相關基因的調控通路[13]。研究表明，不同植物物種在鹽堿脅迫下可能具有不同的miRNA介導的調控策略。高通量測序技術已廣泛用于miRNA鑒定[13-15]，為青稞鹽堿響應的miRNA檢測及研究提供了很好的技術支持。

本研究的目的是在青稞中鑒定鹽堿脅迫響應相關的miRNA，預測其靶標以及探索檢測到的miRNA的功能，并揭示miRNA的表達和調控模式。研究結果可為進一步驗證青稞在鹽堿脅迫下的miRNA-mRNA調控模式提供重要的信息，并為解析青稞和其他農作物中鹽堿脅迫響應的分子機理提供基礎。

1 材料與方法

1.1 植物材料和鹽堿處理

本研究以青稞地方品種0119和0226作為材料，其中0119對鹽堿脅迫具有很高的抗性，而0226對鹽堿脅迫敏感[16]。用2%H2O2將種子表面滅菌30 min，再用無菌水沖洗，放在25℃的黑暗環境下潮濕的薄紙中2 d。然后將發芽的幼苗轉移至溫室中的Hoagland營養液中，光溫條件分別為14 h（光照）/10 h（黑暗）和22℃（光照）/20℃（黑暗）。鹽堿脅迫處理，將2個材料的3周齡幼苗中的1/2移至添加了200 mmol/L Na+（NaHCO3和Na2CO3物質的量之比為1∶1）的培養液中72 h。其余1/2植物作為對照生長在Hoagland營養液中。

1.2 總RNA提取和miRNA文庫的構建

分別取2個材料在鹽堿脅迫（2、8、24、48、72 h）和對照（0、2、8、24、48、72 h）處理后的植株新鮮葉片，使用Trizol試劑（Invitrogen，Carlsbad，CA，美國）從中提取總RNA，并使用RNeasy Plant Mini Kit（Qiagen，Valencia，CA，美國）純化。RNA樣品的濃度和質量使用Nano Drop 2000分光光度計（Thermo Fisher Scientific，Inc.）通過測量A260/A280和A260/A230的比例來確定。將純化的高質量RNA溶于無RNase的水中。取10μg總RNA樣品，通過PAGE凝膠選擇8～30個核苷酸的RNA片段。將RNA片段連接至5'接頭和3'接頭，通過PAGE凝膠再次過濾片段大小。使用One Step Primer Script miRNA cDNA合成試劑盒（Trans Gen）對連接的RNA進行RT-PCR，產生cDNA產物。cDNA文庫通過凝膠純化，并用安捷倫（美國，加利福尼亞州圣克拉拉）生物分析儀2100進行驗證。最終使用HiSeq2000（Illumina）對文庫進行深度測序。每個處理的每個時間點取3個生物學重復。本研究中使用的原始序列數據已提交GSA（Genome Sequence Archive）數據庫[17]及中國科學院北京基因組研究所（BIG）生命與健康大數據中心（2019年會員），注冊號為CRA001345，可在http://bigd.big.ac.cn/gsa上公開訪問。

1.3 miRNA的鑒定和基因表達譜分析

對原始數據去接頭[19]、過濾低質量讀長序列之后，將miRNA讀長檢測序列與青稞參考基因組進行比對[18]，使用miRDeep2（https://github.com/rajewskylab/mirdeep2）軟件鑒定青稞miRNA。在所有樣品中，miRNA讀數均大于66，以此識別候選miRNA。使用BLAST將這些潛在的miRNA與Rfam數據庫進行比對[20]（版本：10.1），鑒定其他ncRNA（非編碼RNA）。為了鑒定已知的miRNA，使用BLAST將潛在的miRNA標簽與miRBase中存在的成熟miRNA進行比對[21]，不超過2個錯配的miRNA被認為是已知的miRNA。DESeq2 R包[22]用于鑒定差異表達的基因。顯著性變化的閾值為fold-change(倍數變化)>2，校正后的P值<0.05。

1.4 靶基因預測和靶功能富集分析

使用psRobot（v1.2,http://omicslab.genetics.ac.cn/psRobot/downloads.php）軟件在默認參數下預測miRNA的靶基因。使用Cytoscape 3.6.1可視化miRNA-mRNA調控網絡[23]。使用Top GO R包[24]對顯著差異表達的miRNA的靶基因進行功能分類和富集分析。

2 結果與分析

2.1 青稞miRNA的鑒定

對于所有66個miRNA文庫，在將低質量讀長序列和接頭序列進行過濾后，共剩余897 112 430條凈讀長序列（clean reads），每個文庫的讀長總數范圍為9 908 031～20 151 566條，唯一讀長數目范圍為1 059 453～5 097 449條。最豐富的miRNA的長度范圍為20～24個核苷酸，最常見的miRNA的長度為21和24個核苷酸。讀長比對率為58.01%～92.76%。在66個樣品中，共13 757個潛在的miRNA有表達。在本分析中，與其他非編碼RNA類型匹配的miRNA標簽被去除，但是沒有潛在的miRNA標簽映射到Rfam數據庫（http://rfam.xfam.org/）。分布在單個染色體上的miRNA數量在1 513（第4染色體）與1 996（第2染色體）之間，并且這些miRNA大致均勻地分布在整個染色體上。在本研究所有文庫中，miRNA的相對表達水平分布為66～7 821 211個讀長。

2.2 已知和新miRNA的BLAST比對

使用BLAST將潛在的miRNA序列與miRBase中公開的miRNA序列進行比較[13]，鑒定已知的miRNA，結果顯示：有225個miRNA在miRBase中是保守的，主要分布在hvu-miR5049和hvu-miR1436家族；其余的13 532個miRNA是推定的新型miRNA。長度在20～24個核苷酸的miRNA數量為11 759，其中28.3%長度為21個核苷酸，53.7%長度為24個核苷酸。miRNA前體的長度為34～113個核苷酸，平均長度為72個核苷酸。

計算所有miRNA的核苷酸組成，核苷酸中腺嘌呤（A）為23.5%，尿嘧啶（U）為26.4%，胞嘧啶（C）為18.6%，鳥嘌呤（G）為31.5%（圖1）。

2.3 鹽堿處理和對照青稞的miRNA表達模式

為了鑒定在鹽堿脅迫下差異表達的miRNA，基于歸一化表達水平的統計比較，對差異表達模式進行了分析。我們通過比較鹽堿處理和對照樣品分析了0226在5個應激階段（2、8、24、48和72 h）中的miRNA表達模式，發現上調和下調的顯著差異表達的miRNA在8和48 h階段增加，在24 h階段減少（圖2-A）。對于0119發現了相同的表達模式（圖2-B）。

比較了0119和0226品種中差異表達的miRNA（圖3），大多數miRNA在0119或0226品種中被特異性上調或下調。

2.4 青稞miRNA的靶基因功能預測

miRNA的功能通常是靶基因的轉錄后調控，由于與特定mRNA的序列互補性，導致蛋白質翻譯抑制或裂解[25]，總共363個差異表達的miRNA具有潛在的靶標候選物，映射到了Swiss Prot數據庫的靶標基因。miRNA預測的靶標數量差異很大，范圍為1～2 107，總共為363個miRNA確定了7 608個唯一的目標候選物。大多數miRNA（224個）具有多個靶基因（2～100個），83個miRNA僅具有1個靶基因，56個miRNA每個具有100多個靶基因。miRNA及其靶基因的一些例子顯示了相反的表達模式。該靶標預測分析表明，已鑒定的靶標基因調控了廣泛的生物學過程。所有這些目標基因均按照基因本體（GO）數據庫進行功能分類（表1）。在細胞成分（cellular component）類別中，前10個GO類別顯示大 多 數 目 標 基 因 與 膜 相 關（GO：0016020，GO：0016021，GO：0031224，GO：0044425和GO：0012505等）。在生物過程（biological process）類別中，最富集的條目與運輸相關（GO：0044765，GO：0006811，GO：0006810，GO：0006812，GO：0055085和GO：0030001）。在分子功能（molecular function）類別中，轉運蛋白活性（GO：0005215）、離子結合（GO：0005507和GO：0030001）、過氧化物酶活性（GO：0004601）、氧化還原酶活性（GO：0016684）和抗氧化活性（GO：0016209）是顯著富集。分析發現，與“膜”相關的miRNA-mRNA調控網絡中，Thb-miRn724a-3p、Thb-miRn724-1-5p和Thb-miRn1158-3p位于網絡的中心。KEGG通路富集分析表明，最富集的6個通路展示在表2中，其中“植物-病原體相互作用”和“植物激素信號轉導”分別富集到248和149個靶向基因。此外，靶基因GO富集分析顯示，目標基因在2 h階段顯著富集了氧化還原酶活性（GO：0016682，GO：0016491和GO：0016679）和陽離子結合（GO：0043169）。24 h時富集的GO條目是ATP酶活性（GO：0042623，GO：0042626和GO：0043492）和水解酶活性（GO：0016820）。在48 h的靶基因GO分析顯示出金屬簇結合（GO：0051540）、核苷酸結合（GO：0000166）和核苷三磷酸酶活性（GO：0017111）富集，而在72 h條件下的靶基因富集氧化還原酶活性（GO：0016682，GO：0016679和GO：0016491）和銅離子結合（GO：0005507）。

表1 所有鹽堿反應性miRNA靶基因的GO功能分類

表2 所有鹽堿反應性miRNA靶基因的KEGG通路富集的類別

3 討論與結論

已有研究報道了對鹽堿條件的兩階段生長反應[26-27]，在玉米和小麥品種中已經證實了這種兩階段生長反應[28-30]。第一階段主要的作用是水分脅迫或滲透作用，由于鹽的積累，第二階段的主要作用是離子性的。在本研究中，顯著差異表達的miRNA數量在24 h下降。這種現象表明第一階段可能發生在青稞用200 mmol/L Na+培養液處理的前24 h。GO富集分析表明，目標基因在2 h階段主要為氧化還原酶活性相關的蛋白，隨后在24、48及72 h，ATP酶活性、水解酶活性、核苷酸結合及銅離子結合相關的蛋白呈現顯著的表達。這些結果與先前研究報道的生理現象[31-33]一致，這是對植物鹽堿脅迫的反應。

本研究鑒定出的大部分miRNA是物種特異性的。在擬南芥[10]、玉米[11]和水稻[34]中，已鑒定出許多與植物對鹽脅迫反應相關的miRNA。我們鑒定了13 757個miRNA，其中只有225個在大麥中保守，2 227個對鹽堿脅迫有響應。這些miRNA可以與靶mRNA的非翻譯區（UTR）或編碼序列（CDS）互補結合，從而降低mRNA的表達水平[1，35]。盡管表達的miRNA可以負調控其靶標mRNA，但由于調控關系的復雜性，其靶標可能表現出多樣化的表達模式。根據miRNA-mRNA調控網絡的靶標基因的GO注釋，我們推斷“膜”的形成與鹽堿脅迫應答有著密切的關系。在以前的研究中，為預測鹽堿脅迫相關的靶標，發現許多靶標是SPL（類表面活性物質B蛋白）、MYB（myb域蛋白）、ARF（生長素應答因子）、AP2（花瓣分生組織發育蛋白APETALA2）、NAC（NAC結構域蛋白）和NF-Y（核轉錄因子Y）。NAC在鹽度、寒冷、ABA和干旱脅迫下受到廣泛調節[36-38]。在這里，鹽和干旱脅迫廣泛調節了一個NF-YA（核轉錄因子Y亞基A）靶標[39,34]；確定了一組介導ARF調節的miRNA，ARF可能在植物對各種非生物脅迫（例如鹽，冷和干旱脅迫）的反應中起重要作用[40-45]；還鑒定了其他基因，包括調節植物生長和發育的SPL、MYB和AP2；認為許多編碼重要酶和蛋白質的靶基因在鹽脅迫反應中也起著重要作用。我們還發現，miRNA交叉調控了幾個靶基因，且單個miRNA可以調控多個靶基因。靈活的調控模式表明，這些相互作用可以調控特定發育階段和生活環境中的特定靶標。該研究可能有助于其他作物及植物鹽堿逆境應答研究方法的完善，以提高農作物適應各種環境條件的能力。