楊利艷， 張玉榮， 楊雅舒， 王美霞， 陳保國(guó)， 趙麗， 張麗光， 王創(chuàng)云

(1.山西師范大學(xué)生命科學(xué)學(xué)院,山西 臨汾 041004;2.山西省農(nóng)業(yè)科學(xué)院作物科學(xué)研究所,山西 太原 030031)

同一mRNA前體在產(chǎn)生成熟mRNA的過(guò)程中會(huì)在不同的排列方式下剪除內(nèi)含子，這種機(jī)制最早由GILBERT[1]在1978年首先提出，稱之為可變剪接(alternative splicing, AS) ?？勺兗艚颖徽J(rèn)為是導(dǎo)致蛋白質(zhì)功能多樣性的重要原因之一,它使一個(gè)基因產(chǎn)生結(jié)構(gòu)和功能不同的mRNA和蛋白質(zhì)[2-3]。選擇性剪接事件在真核生物基因的轉(zhuǎn)錄過(guò)程中廣泛發(fā)生，盡管最初有研究認(rèn)為高等生物僅大約5%的基因發(fā)生可變剪接[4]，但隨后在人和其他幾種重要的模式生物研究中發(fā)現(xiàn)，可變剪接的發(fā)生概率遠(yuǎn)大于5%。在人類(Homosapiens)、線蟲(chóng)(Caenorhabditiselegans)、果蠅(Drosophilamelanogaster)、擬南芥(Arabidopsisthaliana)中發(fā)生AS事件的比例分別高達(dá)95%，71%，60.7%，61%[5-8]。基因組數(shù)據(jù)顯示，由于選擇性剪接的發(fā)生，擬南芥中約有42%的內(nèi)含子基因參與了擬南芥的發(fā)育、生物及非生物脅迫反應(yīng)[9-10]?？勺兗艚油ㄟ^(guò)產(chǎn)生不同轉(zhuǎn)錄本，導(dǎo)致蛋白質(zhì)的結(jié)構(gòu)存在差異，并可能在蛋白質(zhì)的功能和特異性方面進(jìn)行修飾。可變剪接轉(zhuǎn)錄本編碼的蛋白質(zhì)中，約有超過(guò)一半與信號(hào)轉(zhuǎn)導(dǎo)有關(guān)，且大部分可變剪接區(qū)域包含了蛋白質(zhì)和其結(jié)合模體的相互作用位點(diǎn)，這些模體包括配體、DNA/RNA和其他的蛋白質(zhì)[11]。通過(guò)可變剪接通過(guò)形成不同轉(zhuǎn)錄本通過(guò)調(diào)節(jié)信號(hào)途徑中某些關(guān)鍵基因的表達(dá)，從而影響動(dòng)植物生長(zhǎng)發(fā)育、生理代謝以及抗逆能力。煙草中N抗性基因是植物抗性基因LRR的成員，通過(guò)選擇性剪接可以形成2個(gè)轉(zhuǎn)錄本NS，NL。研究發(fā)現(xiàn)，轉(zhuǎn)cDNA-NL和cDNA-NS煙草均不能對(duì)煙草花葉病毒(TMV)表現(xiàn)出完全抗性，只有cDNA-NS-內(nèi)含子和含有N基因組序列的轉(zhuǎn)基因植株，同時(shí)編碼NS和NL轉(zhuǎn)錄本才表現(xiàn)出對(duì)TMV的完全抗性，即煙草N基因發(fā)生可變剪接所形成的2個(gè)轉(zhuǎn)錄本對(duì)于TMV的完全抗性是必需的[12]。擬南芥中蛋白磷酸激酶HAB1在ABA信號(hào)轉(zhuǎn)導(dǎo)途徑中起著關(guān)鍵作用，但若RBM25剪接因子丟失，其HAB1基因的表達(dá)量會(huì)下降，從而阻止ABA信號(hào)傳導(dǎo)[13]。此外，可變剪接參與了擬南芥開(kāi)花時(shí)間的控制，flowering time control locus a(FCA)前mRNA的選擇性剪接參與了花期轉(zhuǎn)換的調(diào)控[14-15]。由此可見(jiàn)，可變剪接是產(chǎn)生蛋白質(zhì)結(jié)構(gòu)和功能多樣性的一個(gè)重要機(jī)制，鑒定基因可變剪接事件有利于全面研究基因的功能。近幾年，已有多種方法用以鑒定AS事件[16]，包括基因芯片及RNA-seq 技術(shù)[17-20]。目前，利用RNA-Seq技術(shù)已分析了多個(gè)物種基因組的AS事件，包括人[5]、小鼠(Musmusculus)[21]、牛(Bovini)[22]、擬南芥[8]、中華獼猴桃(Actinidiachinensis)[23]等。本研究分別對(duì)黏蟲(chóng)取食前后玉米的葉片進(jìn)行轉(zhuǎn)錄組測(cè)序，將兩者發(fā)生可變剪接事件的基因進(jìn)行比較分析，旨在挖掘響應(yīng)抗蟲(chóng)取食、同時(shí)發(fā)生AS事件的基因，從而探究可變剪接對(duì)植物抗蟲(chóng)響應(yīng)基因的調(diào)控作用。

1 材料與方法

1.1 材料

玉米(Zeamays)自交系B 73由中國(guó)農(nóng)業(yè)科學(xué)院作物科學(xué)研究所提供。分別用2%次氯酸鈉溶液消毒種子20 min，用蒸餾水清洗干凈后，置于濕潤(rùn)濾紙上，于22～23 ℃培養(yǎng)箱，萌發(fā)后，置于25 ℃、相對(duì)濕度65%～80%、光暗周期16 h/8 h的光照培養(yǎng)箱，長(zhǎng)出2 葉后，移至花盆，常規(guī)管理至5～6葉。

黏蟲(chóng)(Mythimnaseparata)來(lái)自中國(guó)農(nóng)業(yè)科學(xué)院植保所培養(yǎng)的2～3齡幼蟲(chóng)。

1.2 方法

待幼苗長(zhǎng)至六葉期，選長(zhǎng)勢(shì)一致、健康的幼苗，每株幼苗即為1個(gè)重復(fù)，每組設(shè)置3個(gè)重復(fù)。每株材料上接種4頭黏蟲(chóng)，使其均勻分布在植株上，12 h后取材，以不接蟲(chóng)材料為對(duì)照(M0)。各處理組(M12)植株置于用木架(長(zhǎng)50 cm、寬50 cm、高80 cm)中進(jìn)行取食試驗(yàn)，外罩一層薄紗以防黏蟲(chóng)逃逸。12 h后取相同葉位材料并用液氮和干冰保存，提取總RNA，由生物公司進(jìn)行轉(zhuǎn)錄組測(cè)序。

1.3 RNA的提取

采用Trizol法分別提取玉米試驗(yàn)組M12、對(duì)照組M0中葉片總RNA。將所提取的RNA以1% 的瓊脂糖凝膠進(jìn)行檢測(cè)，并利用超微量紫外分光光度計(jì)對(duì)RNA進(jìn)行定量分析。

1.4 RNA-seq文庫(kù)構(gòu)建和測(cè)序

RNA-seq文庫(kù)構(gòu)建和測(cè)序均由上海生工生物工程有限公司完成，分別取3個(gè)重復(fù)樣品的總RNA等量混合后構(gòu)建2個(gè)文庫(kù)，即M12(試驗(yàn)組)、M0(對(duì)照組)，進(jìn)行Illumina HiseqTM測(cè)序。測(cè)序得到的原始數(shù)據(jù)，里面含有帶接頭的、低質(zhì)量的序列。為了保證信息分析質(zhì)量，使用Trimmomatic進(jìn)行對(duì)原始數(shù)據(jù)進(jìn)行過(guò)濾，去除帶N堿基的序列(N表示無(wú)法確定堿基信息)、 去除reads中的接頭序列、去除低質(zhì)堿基(Q<20)、使用滑窗法去除reads尾部質(zhì)量值在20以下的堿基(窗口大小為5 bp)，得到Clean數(shù)據(jù)。本研究?jī)H使用得到的RNA-seq測(cè)序數(shù)據(jù)鑒定響應(yīng)黏蟲(chóng)脅迫的發(fā)生可變剪接事件的差異表達(dá)基因，分析M12與M0發(fā)生可變剪接事件的共性和特異性，研究可變剪接對(duì)抗蟲(chóng)的調(diào)控作用。

1.5 可變剪接事件的鑒定

使用HISAT2將質(zhì)控后的測(cè)序序列與參考基因組進(jìn)行比對(duì)，再運(yùn)用 Cufflinks 軟件[24-25]預(yù)測(cè)基因模型，使用ASprofile分析[26]，并根據(jù)每個(gè)樣本預(yù)測(cè)出來(lái)的基因模型，對(duì)可變剪切事件進(jìn)行分類。ASprofile分析存在12種可變剪接事件類型(圖1)，即單外顯子跳躍(skipped exon，SKIP)、單外顯子跳躍(模糊邊界)(approximate SKIP，XSKIP)、多外顯子跳躍(multi-exon SKIP，MSKIP)、多外顯子跳躍(模糊邊界)(approximate MSKIP，XMSKIP)、單內(nèi)含子滯留(intron retention，IR)、單內(nèi)含子滯留(模糊邊界)(approximate IR，XIR)、多內(nèi)含子滯留(multi-IR，MIR)、多內(nèi)含子滯留(模糊邊界)(approximate MIR， XMIR)、可變5′ 或3′ 端剪切(alternative exon ends，AE)、可變 5′ 或3′ 端剪切(模糊邊界)(approximate AE，XAE)、第1個(gè)外顯子可變剪切(alternative 5′first exon，TSS)、最后1個(gè)外顯子可變剪切(alternative 3′ last exon，TTS)。

A．單外顯子跳躍和多外顯子跳躍(灰色和紅色均代表外顯子)；B. 單內(nèi)含子滯留和多內(nèi)含子滯留(灰色代表外顯子、紅色代表內(nèi)含子)；C. 可變5′ 或3′ 端剪切(灰色和紅色均代表外顯子)；D. 第1個(gè)外顯子可變剪切(灰色和紅色均代表外顯子)；E. 最后1個(gè)外顯子可變剪切(灰色和紅色均代表外顯子)。

A. Skipped exon (SKIP) and multi-exon SKIP(MSKIP)(Both gray and red represent exons); B. Intron retention(IR) and multi-IR (MIR)(Gray for exons and red for introns); C. Alternative exon ends (AE)(Both gray and red represent exons); D. Alternative 5′ first exon(TSS)(Both gray and red represent exons); E. Alternative 3′last exon (TTS)(Both gray and red represent exons).

圖1 可變剪接事件的分類

Fig.1 Classification of alternative splicing events

1.6 抗蟲(chóng)相關(guān)基因分析及GO功能富集分析

為篩選出響應(yīng)抗蟲(chóng)脅迫的可變剪接基因，對(duì)M12中特有可變剪接基因用DEseq進(jìn)行差異分析。為了得到表達(dá)顯著差異的基因，將篩選條件為qvalue<0.05(q為校正后的P值<0.05)且差異倍數(shù)|FoldChange|> 2(表達(dá)量差異倍數(shù)>2)。根據(jù)差異基因的表達(dá)量FPKM(fragments per kilobase of exon per million reads mapped)值，采用Heatmap軟件繪制它們?cè)诟鳂颖局械谋磉_(dá)量熱圖。

為確定其基因功能，將顯著差異表達(dá)的可變剪接基因ID比對(duì)到DAVID數(shù)據(jù)庫(kù)，獲得根據(jù)GO分類的分子功能(molecular function)、生物學(xué)過(guò)程(biological process)和細(xì)胞組成(cellular component)的GO編號(hào)，再使用clusterProfiler進(jìn)行功能富集分析(qvalue<0.05，認(rèn)為該功能存在顯著富集)。

2 結(jié)果與分析

2.1 RNA的提取和檢測(cè)

將所提取的RNA以1% 的瓊脂糖凝膠進(jìn)行檢測(cè)(圖2)，并使用超微量紫外分光光度計(jì)對(duì)RNA進(jìn)行定量測(cè)定，其中A260/A280均為1.9～2.1。

2.2 可變剪接事件的鑒定

將M0與M12 2種材料的轉(zhuǎn)錄組數(shù)據(jù)分別與玉米基因組進(jìn)行比對(duì)后，M0中15 701個(gè)基因發(fā)生了可剪接事件，共發(fā)生79 363個(gè)， M12中11 791個(gè)基因發(fā)生了AS， 共發(fā)生了39 385個(gè)。

在不同的可變剪接事件類型中，M0以TSS，TTS，IR，AE 4種類型為主，數(shù)量分別為21 148，18 805，10 252，8 283，占所有可變剪接事件的26.65%，23.69%，12.92%，10.44%，其他的8種類型占比為1.2% ～5.7%。在經(jīng)黏蟲(chóng)取食后，M12主要發(fā)生的可變剪接事件類型沒(méi)有變化，仍為T(mén)SS，TTS，IR，AE等4種，但可變剪接事件發(fā)生數(shù)量有所減少，分別為14 421，13 625，2 735，2 706，其比例為36.62%，34.59%，6.94%，6.87%，余下8種AS事件類型比例為0.2% ～ 4.3%(表1，圖3)。

M：DNA mark；1, 2, 3: 樣品RNA 。

表1 M0與M12基因組中鑒定的可變剪接事件

圖3 M0與M12中各類型可變剪接事件所占比例

對(duì)M0與M12中發(fā)生可變剪接的基因進(jìn)行統(tǒng)計(jì)和比較發(fā)現(xiàn)，M0與M12之間有10 670個(gè)可變剪接事件是重合的。然而，在M0中存在5 031個(gè)新基因發(fā)生可變剪接，這些基因在M12中沒(méi)有發(fā)生AS事件。在黏蟲(chóng)取食后的M12基因組中存在相同現(xiàn)象，有1 121個(gè)新基因發(fā)生可變剪接，這些基因沒(méi)有在M0中發(fā)生可變剪接事件(圖4)。

2.3 M12特有AS事件基因差異表達(dá)分析

對(duì)M12 中1 121個(gè)可變剪接基因進(jìn)行差異表達(dá)分析，得到177個(gè)顯著差異表達(dá)基因。差異基因的聚類分析見(jiàn)圖5。其中，上調(diào)基因有142個(gè)，下調(diào)基因有35個(gè)。在這177個(gè)差異表達(dá)基因發(fā)生的可變剪接事件類型中，主要為第1個(gè)外顯子可變剪切(TSS)和最后1個(gè)外顯子可變剪切(TTS)，表明玉米響應(yīng)抗蟲(chóng)基因的表達(dá)與這2種剪切方式之間存在密切的聯(lián)系。

值得關(guān)注的是，M12中的177個(gè)發(fā)生AS事件的差異表達(dá)基因有多個(gè)基因?yàn)榭瓜x(chóng)相關(guān)基因。這些基因分別是Zm00001d048736,Zm00001d004555。Zm00001d022453，主要參與生物堿的合成；Zm00001d039394，主要參與萜類物質(zhì)的合成；Zm00001d032079，主要參與植物-病原互作信號(hào)通路。還有一些與RNA剪接、RNA轉(zhuǎn)運(yùn)、RNA降解相關(guān)的基因。

2.4 可變剪接基因GO注釋及功能富集分析

為了解可變剪接對(duì)玉米抗蟲(chóng)能力的調(diào)控，對(duì)經(jīng)黏蟲(chóng)取食的M12基因組中這差異表達(dá)的177個(gè)發(fā)生AS事件的基因進(jìn)行了GO注釋和功能分析，統(tǒng)計(jì)P≤1×10-10的GO項(xiàng)，發(fā)現(xiàn)有132個(gè)基因存在GO注釋(圖6)。其中，在生物學(xué)過(guò)程、分子功能和細(xì)胞構(gòu)成中顯著(qvalue<0.05)富集的GO項(xiàng)分別有84，93，25個(gè)。在生物學(xué)過(guò)程中，發(fā)生可變剪接的差異基因主要集中于氧化還原過(guò)程(GO:0055114)和代謝過(guò)程(GO:0008152)；在分子功能中，較多可變剪接基因被富集于金屬離子的結(jié)合(GO:0046872)、氧化還原酶的活性(GO:0016491)相關(guān)GO項(xiàng)中；在細(xì)胞構(gòu)成中，發(fā)生可變剪接的基因主要集中于膜的組成部分(GO:0016021)。

圖4 M0和M12共有和特有發(fā)生可變剪接基因數(shù)量

數(shù)字代表注釋到該功能下的可變剪接基因數(shù)目.A：生物學(xué)過(guò)程；B：細(xì)胞成分；C：分子功能。The numbers represent the AS gene number annotated to the corresponding function.A: Biological process; B: Cellular components; C: Molecular functions.

3 結(jié)論與討論

可變剪接在植物的生長(zhǎng)發(fā)育、病蟲(chóng)害防御中扮演著重要的角色。KANG等[27]發(fā)現(xiàn)，擬南芥miR400基因在位于內(nèi)含子306 bp所發(fā)生的1個(gè)可變剪接事件是由熱應(yīng)激引起的。玉米ZmPP2C26基因轉(zhuǎn)錄后加工過(guò)程中會(huì)發(fā)生內(nèi)含子保留型可變剪接，產(chǎn)生2個(gè)轉(zhuǎn)錄本，第1個(gè)轉(zhuǎn)錄本內(nèi)含子盡管沒(méi)有改變其編碼蛋白的PP2C功能，但有可能使其作用由細(xì)胞核和細(xì)胞膜延伸至葉綠體[28]。

本試驗(yàn)中，與M0相比，黏蟲(chóng)取食后的M12基因組中有1 121個(gè)新基因發(fā)生可變剪接。對(duì)差異表達(dá)的可變剪接基因進(jìn)行GO分析表明，分子功能結(jié)果主要與氧化還原酶活性、轉(zhuǎn)移酶活性、水解酶活性、DNA結(jié)合、ATP結(jié)合相關(guān)，而這些功能均與植株的抗性相關(guān)。如研究發(fā)現(xiàn)，過(guò)表達(dá)2-半胱氨酸氧化還原酶基因(2-Cysperoxiredoxins,Zj2-CP)的擬南芥與野生型擬南芥相比，過(guò)表達(dá)Zj2-CP基因植株增加了對(duì)干旱和鹽脅迫的敏感性[29]，依賴于2OG-Fe(II)的氧化還原酶可以特異性調(diào)控種子的萌發(fā)過(guò)程[30]。對(duì)擬南芥精氨酸酶(AtArgAH1和AtArgAH2)基因與植物耐鹽性之間的關(guān)系分析發(fā)現(xiàn)，2個(gè)精氨酸酶基因(AtArgAH1和AtArgAH2)表達(dá)與鹽脅迫存在一定的應(yīng)答關(guān)系[31]。生物過(guò)程結(jié)果主要與氧化還原過(guò)程、代謝過(guò)程、以DNA為模板的轉(zhuǎn)錄及調(diào)控相關(guān)。這些功能與次生代謝物的產(chǎn)生、植株抗性基因表達(dá)調(diào)控相關(guān)，從而增加植物的抵抗力。植物能夠產(chǎn)生類萜、生物堿、多炔、類黃酮、精油、獨(dú)特的氨基酸和糖類等多種性質(zhì)各異的次生代謝物, 這些代謝物在作物保護(hù)上有著廣泛用途[32]。

本研究發(fā)現(xiàn)的Zm00001d048736基因發(fā)生2個(gè)AE可變剪接事件，TTS和TSS事件各1個(gè)，主要參與生物堿的合成。大多數(shù)生物堿都具有毒性作用，是許多植物殺蟲(chóng)劑中最常見(jiàn)的主要成分[33]。生物堿對(duì)昆蟲(chóng)的影響是多方面的，可影響昆蟲(chóng)的神經(jīng)系統(tǒng)，刺激或抑制昆蟲(chóng)取食，影響昆蟲(chóng)產(chǎn)卵等[34]。Zm00001d039394基因發(fā)生TTS和TSS事件各1個(gè)，主要參與萜類物質(zhì)的合成。利馬豆葉暴露于蜘蛛螨取食植物后所誘導(dǎo)的任何萜類揮發(fā)物質(zhì)，可誘導(dǎo)利馬豆葉PR蛋白以及LOX，PAL和FPS[1]等基因的表達(dá)[35-36]，將植物暴露于外源茉莉酸中，這些基因的表達(dá)水平也有所提高，表明這些基因會(huì)對(duì)植物的損傷脅迫作出響應(yīng)[37-38]，即萜類物質(zhì)可以誘導(dǎo)植物防御基因的表達(dá)。還有一些與RNA剪接、RNA轉(zhuǎn)運(yùn)、RNA降解相關(guān)的基因也發(fā)生可變剪接，這些基因通過(guò)產(chǎn)生不同轉(zhuǎn)錄本對(duì)其功能產(chǎn)生影響，通過(guò)影響基因的轉(zhuǎn)錄翻譯水平而發(fā)揮功能。以上表明，可變剪接對(duì)于植物抗蟲(chóng)的調(diào)控不僅發(fā)生于植物的誘導(dǎo)抗性，還對(duì)參與轉(zhuǎn)錄的基因有所調(diào)控。

本試驗(yàn)利用轉(zhuǎn)錄組樣品的高通量測(cè)序數(shù)據(jù)鑒定了玉米蟲(chóng)食脅迫后所特異發(fā)生可變剪接的差異表達(dá)基因，通過(guò)尋找可變剪接事件的發(fā)生規(guī)律和可變剪接參與植物基因表達(dá)調(diào)控的途徑，為探究可變剪接對(duì)植物抗蟲(chóng)能力的調(diào)控，挖掘植物對(duì)蟲(chóng)害的適應(yīng)機(jī)制提供了基礎(chǔ)研究數(shù)據(jù)。