藜麥葉片響應低氮脅迫基因的可變剪接事件分析

楊利艷，楊小蘭，朱滿喜，張玉榮，楊雅舒，王創云

(1.山西師范大學生命科學學院，中國 臨汾 041000；2.山西農業大學農學院，中國 太谷 030800)

藜麥(ChenopodiumquinoaWilld)為一年生雙子葉植物，莧科藜亞科藜屬，是一種起源于南美洲安第斯區的古老物種[1]。藜麥種子富含無谷蛋白以及人體無法合成的8種必需氨基酸，且鐵、鈣、磷含量豐富，具有良好的營養健康效益[2，3]。藜麥不僅有極高的營養價值，還具有很強的抗逆性，因此受到世界的關注并被廣泛引種，是未來最具潛力的農作物之一[4]。

氮元素(N)是細胞結構不可或缺的一部分，擔負著植物的生長發育及機體調控的重要責任[5]。氮素約占植物干重的0.3%～5%，缺氮或低氮是導致農作物生長不良、減產的重要原因[6]。氮脅迫條件下，植物會通過根的伸長來增強氮的吸收[7]，在甘藍型油菜[8]、小麥[9]、水稻[10]中皆有報道。然而目前對藜麥低氮響應的研究鮮有報道。

可變剪接(alternative splicing，AS)是一種轉錄后調節機制。該機制可使蛋白的最終產物表現出不同的功能和結構特性，或者在相同的細胞中可能由于表達水平的不同而導致表型發生變化，因此對增強機體內轉錄組的可塑性、基因的表達以及蛋白質結構的多樣性起重要作用[11,12]。目前，利用RNA-seq 技術對可變剪接的分析有不少報道，可變剪接基因(alternatively spliced genes，ASG)在一些生物中占較高比例，如人(95%)、線蟲(71%)[13,14]，在禽類[15]、甲殼類動物[16]中也發現了新的剪接變體。作為模式植物，擬南芥AS事件的發生比例也較高，為61%[17]，約有42%的內含子參與擬南芥的生長發育及脅迫響應[18,19]。KHOKHAR等[20]發現在不同生態類型的擬南芥中，許多反式剪接數量性狀位點(sqtl)在晝夜節律鐘、開花和應激反應基因中富集，表明差異可變剪接在調節這些重要過程中發揮了作用。KANNO等[21]在擬南芥中通過對一個prp4ka突變體的定量磷酸化蛋白組學分析，檢測到幾個富含絲氨酸/精氨酸的蛋白的磷酸化變化，這些蛋白調節組成型和選擇性剪接以及其他剪接相關的因素。因此，研究可變剪接有利于對基因的深入了解。本研究對不同氮處理的藜麥葉片進行轉錄組測序，統計并分析可變剪接基因，旨在挖掘響應低氮脅迫的ASG，從而探究AS對植物低氮響應基因的調控作用。

1 材料與方法

1.1 材料

藜麥ZK7種子來自于山西省農業科學院作物科學研究所。將種子浸泡在10% NaClO溶液中消毒15 min，用蒸餾水沖洗3～5次，置于放有濕潤濾紙的培養皿中，在25 ℃下培養發芽。待芽長約0.5 cm時，將其移入干凈的砂子中置于光照培養箱育苗，光暗周期16 h/8 h，溫度18～22 ℃。待藜麥幼苗生長至2葉，采用霍格蘭全營養液培養，每3天澆一次，長至6葉期進行以下處理。

1.2 實驗設置

選取長勢相同的幼苗進行兩組處理，即對照(CK)：全營養液培養(N濃度為7.1 mmol·L-1)；低氮培養(LN)：氮元素濃度為全營養液的20% (N濃度1.42 mmol·L-1)。每個處理設置3個重復，處理5天后取相同葉位材料，液氮速凍后保存。提取總RNA，轉錄組數據委托生物公司測定。

1.3 總RNA的提取與質量檢測

采用試劑盒RNAprep Pure Plant Kit(TIANGEN)的方法，分別提取藜麥葉片的總RNA。用Nanodrop ND-1000檢測RNA的濃度，Caliper LabChip GX精確檢測RNA的純度和完整性(RIN)。

1.4 RNA-seq文庫構建和測序

RNA-seq文庫構建和測序由北京百邁客生物科技公司完成，將對照組(CK)和低氮脅迫組(LN)兩個樣品檢測合格后構建文庫，之后進行Illumina HiseqTM測序。測序去除含有接頭和低質量的Reads，獲得高質量Clean數據。本研究通過分析已獲得的Clean數據鑒定在低氮脅迫條件下發生AS事件的差異表達基因，從中找尋CK和LN中ASG的共性和差異性，進而理解AS對低氮脅迫的調控作用。

1.5 可變剪接事件的鑒定

采用StringTie[22]對Hisat2的比對結果進行拼接，通過ASprofile[23]軟件分析可變剪接類型，細分為12類，并獲取CK和LN中相應的可變剪接類型及表達量[24]。

1.6 低氮脅迫相關基因GO分類及KEGG功能富集分析

為篩選響應低氮脅迫的ASG，對LN中特有ASG進行差異分析[25]。并以Fold Change(表達量差異倍數)≥2，FDR(錯誤發現率)<0.01且q<0.05(q為校正后的P值)為篩選標準，由于不同處理組基因存在差異，根據CK和LN組中基因的表達量值繪制各樣本中的基因表達量熱圖，采用Heatmap軟件作圖。GO分類根據GO數據庫(Gene Ontology Consortium)獲得分子功能( molecular function )、生物學過程(biological process )和細胞組成(cellular component)的GO編號，功能富集采用cluster Profiler分析(q<0.05)。KEGG通路分析以Pathway為單位，應用超幾何檢驗，統計并分析其與整個基因組背景相比顯著富集的Pathway。

2 結果與分析

2.1 RNA的提取與檢測

將提取的RNA(RIN值≥7.0)采用Caliper Lab Chip GX檢測，28S/18S均在2.0～2.1，A260/A280均在2.1～2.2(圖1)。

圖1 CK與LN總RNA的GX結果Fig. 1 The GX results of total RNA in CK and LN

2.2 可變剪接事件的鑒定

對CK和LN兩個樣品的轉錄組數據分析后發現CK中26 094個基因發生了AS，共發生68 764個，LN中25 562個基因發生AS，共發生66 611個。

在各種AS類型中，CK組以TSS和TTS兩種類型為主，分別占所有可變剪接事件的40.95%和40.87%；此外AE，IR及SKIP三種類型分別占所有剪接事件的7.06%，3.97%及3.98%；剩下7種類型共占3.17%。低氮處理組中，LN組的可變剪接事件仍以TSS和TTS兩種類型為主，但較CK組有所增加，分別占41.45%和41.11%；AE，IR及SKIP三種類型分別占6.18%，4.58%及3.68%，其中IR剪接類型較CK組增加15.4%；剩下7種類型共占3.05%(表1)。

表1 CK和LN基因組中鑒定的可變剪接事件類型及數量

圖2 CK和LN中可變剪接重合及特有的基因數量 Fig. 2 Number of overlapped and specific alternative splicing genes in CK and LN

對CK與LN中發生的ASG進行整理對比后發現，有24 520個AS發生重合。然而，有1 573個新基因僅在CK中發生AS。同時，LN組中有1 041個新基因發生AS，在CK中卻沒有發生(圖2)。

2.3 LN特有可變剪接基因差異表達分析

對LN中1 041個ASG進行分析，差異表達顯著的有147個基因，對其聚類分析見圖3。這147個基因中AS類型以TSS和TTS為主，其中，上調基因有72個，下調基因有74個，還有一個基因(NW_018743732.1)在TSS事件中下調，在TTS事件中上調。表明藜麥在低氮條件下的基因表達與這兩種可變剪接事件密切相關。同時，這147個可變剪接基因的功能主要有無機離子運輸和傳遞，次生代謝物的合成(NW_018742204.1及NW_018742461.1)、DNA解旋酶及修復蛋白(NW_018742412.1)、翻譯后修飾(NW_018742673.1)、信號轉導(NW_018743014.1)、細胞骨架的形成(NW_018742814.1)及細胞周期調控(NW_018742966.1)。

圖3 LN差異可變剪接基因表達聚類熱圖 Fig. 3 Cluster thermography of differentially expressed genes

2.4 LN中差異可變剪接基因GO注釋及功能分析

對LN基因組中147個可變剪接基因進行GO注釋和功能分析，發現39個GO注釋的基因(圖4) 。其中，顯著富集(q<0. 05)的GO項在生物學過程中有81個，在細胞構成中有26個，分子功能中有60個。在生物學過程中，發生可變剪接的差異基因主要集中于代謝過程(GO:0008152)和刺激反應(GO:0050896)；在細胞構成中，發生可變剪接的差異基因主要集中于核(GO:0005634)、細胞質(GO:0005737)、細胞壁(GO:0005618)以及質外體(GO:0048046)上；在分子功能中，發生可變剪接的差異基因主要富集于轉移酶活性(GO:0016740)和水解酶活性(GO:0016787)相關GO項。

2.5 差異可變剪接基因KEGG通路分析

對LN中147個差異可變剪接基因進行KEGG通路分析(表2)，發現顯著富集到DNA復制、異源末端鏈接、抗壞血酸和醛酸代謝、錯配修復系統等通路上。DNA復制、異源末端鏈接參與新細胞的合成；抗壞血酸和醛酸代謝參與植物逆境響應、細胞生長；錯配修復系統參與DNA雙鏈上錯配堿基對的矯正，泛素介導的蛋白水解調控細胞周期，囊泡運輸的相互作用參與分子的轉運，谷胱甘肽代謝參與抗氧化和解毒作用。

注:數字代表在該功能下注釋的基因數；A：生物學過程；B：細胞構成；C：分子功能。圖4 LN中特有可變剪接差異基因的GO注釋Fig. 4 GO annotations for specific differential expressed alternative splicing genes in LN

3 討論

選擇性剪接在廣大真核生物產生潛在的基因組信息方面發揮著重要作用，它是一個重要的轉錄后調節機制，可以增加蛋白質多樣性，影響mRNA的穩定性[26]。然而，關于AS如何參與藜麥對低氮脅迫的響應，人們知之甚少。在本實驗中，筆者在LN組中共識別到66 611個AS事件，包括TSS，TTS，AE，IR，SKIP等12種。其中TSS及TTS事件最為豐富，其次是AE，IR和SKIP。與CK組相比，鑒別出低氮處理組(LN)基因組中有1 041個新基因發生AS，這些基因的生物學過程富集在代謝過程和刺激反應等過程中。LEVIATAN等[27]發現擬南芥在低溫脅迫下，其表達水平保持相對不變，但其剪接方式在異構體的相對豐度上有明顯變化。大多數低溫調控的選擇性剪接將提前終止密碼子(PTC)引入到轉錄本中，通過mRNA衰變(NMD)過程產生潛在的降解目標，這可能導致新蛋白質功能改變或產生負面影響。FANG等[28]發現水稻中11個OsAAT選擇性剪接基因，其中1個基因剪接變異的表達受氮(N)的調控，5個基因在各種形態氮濃度較高時表達上調，9個基因在硝酸鹽濃度較高時表達上調。剪接變異體的mRNA水平在水稻葉和穗發育過程中也受到調控，結果表明水稻OsAAT家族的可變剪接變異體在不同氮素條件下的表達率均有變化。本實驗發現，藜麥響應低氮的差異可變剪接基因富集于與轉移酶活性和水解酶活性相關的通路。GIARETTA等[29]研究表明由GGT1和GGT2編碼的Apoplastic-谷氨酰胺轉移酶活性對植物營養和生殖發育具有重要意義，因此筆者推測轉移酶對藜麥緩解低氮脅迫起一定作用。對LN中差異可變剪接基因進行KEGG分析，發現顯著富集到DNA復制、異源末端鏈接、抗壞血酸和醛酸代謝、錯配修復系統等通路上，這些信號通路與藜麥氮吸收密切相關，表明在藜麥生長過程中可能通過差異可變剪接基因調節氮吸收來響應低氮脅迫。

4 結論

本試驗利用轉錄組樣品的高通量測序數據鑒定了藜麥響應低氮的特異可變剪接的差異表達基因，藜麥葉片主要通過TSS和TTS兩種剪接方式形成差異基因，其中NW_018742556.1,NW_018742616.1,NW_018744348.1,NW_018742616.1,NW_018742733.1等基因通過DNA復制和異源末端鏈接等通路來響應低氮。