觀賞海棠‘火焰’果皮中可變剪接基因分析

段詩瑤，田 佶，張 杰

(北京農學院植物科學技術學院/農業(yè)應用新技術北京市重點實驗室，北京 102206)

蘋果(Malusdomestica)是中國栽培面積最廣的果樹[1]。蘋果果實的品質包括果實的大小、形狀、色澤、香氣、糖分和有機酸的比例等，其直接影響蘋果的經濟價值。在產量大幅度提升的同時，蘋果品質的下降嚴重阻礙中國蘋果產業(yè)的發(fā)展[2-3]。提高蘋果果實品質是維持其產業(yè)發(fā)展的重要發(fā)展方向，而果皮色澤是果實品質的重要判斷標準。蘋果果實發(fā)育從開花到成熟一般需要150 d左右，是一個漫長的過程，有大量基因參與其成熟的調控。通過基因組與表觀組甲基化測序分析，發(fā)現(xiàn)SPL13、ACS8等基因在果實發(fā)育早期通過控制細胞數(shù)量的變化而可能對成熟期果實大小產生影響[4]。MdMYB9和MdMYB11則通過與bHLH3、bHLH33相互作用來促進蘋果果皮花色素苷的合成[5-6]。利用單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)標記，可以了解到蘋果基因組單核苷酸的變異情況。蘋果MdCoLBD1/2基因序列通過SNP分析，顯示出LBD基因家族各基因的特異性，進而導致各基因功能的特異性[7]。SNP4299和SNP4432與蘋果抗炭疽菌葉枯病相關的基因位點緊密連鎖[8]。

可變剪接在植物生長發(fā)育過程中起著非常重要的調控作用，其通過不同的剪接方式使同一個mRNA前體產生兩個或更多成熟mRNA，進而產生多種不同異構體來增加蛋白的多樣性[9-10]。擬南芥中開花阻遏因子FLOWERINGLOCUSM通過可變剪接產生一種被稱為FLM-β的剪接體，其受溫度的調節(jié)進而影響擬南芥的開花時間[11]。番茄MADS-box基因Lemads1在轉錄時存在4種可變剪接體用于調節(jié)萼片及果實的形成[12]。可變剪接不僅可以通過產生蛋白異構體來調控植物生長發(fā)育，還可以通過改變轉錄本的可讀框使終止密碼子提前產生，進入降解途徑。目前針對蘋果果皮發(fā)育進程中的可變剪切研究較少[13]。

觀賞海棠作為砧木在蘋果屬植物的生產及育種中發(fā)揮重要作用。其葉片和果皮的多色性為研究果皮發(fā)育過程提供良好的研究材料。該研究通過轉錄組測序(RNA Sequencing，RNA-seq)技術對觀賞海棠‘火焰’5個發(fā)育階段果皮進行轉錄組測序分析，比較‘火焰’在不同發(fā)育階段與果皮成熟相關的共有及特有可變剪接基因，為探索蘋果果皮發(fā)育的遺傳調控機理提供新的思路，進而為蘋果屬植物的育種奠定理論基礎。

1 材料與方法

1.1 測序數(shù)據(jù)

測序數(shù)據(jù)為課題組前期測定的蘋果屬觀賞海棠‘火焰’5個果實發(fā)育時期測序數(shù)據(jù)(NCBI登錄號：PRJNA546083)。觀賞海棠‘火焰’果實選自北京農學院觀賞海棠種質資源圃(40.l°N，116.6°E)長勢大小一致的東南方向著生的5個不同發(fā)育時期的綠果品種‘火焰’(Maluscv. ‘Flame’)。在花后35、60、95、120和150 d進行果實的采摘，分別命名為S1、S2、S3、S4和S5。

1.2 建庫及分析流程

分別利用Qubit?2.0熒光儀(Life Technologies，CA)和Nano Photometer?2000分光光度計(IMPLEN,CA)對觀賞海棠果實樣品中提取的總RNA進行濃度和純度的測定。質檢合格后利用富含Oligo dT的磁珠富集mRNA，經過mRNA片段化后，將mRNA反轉成cDNA，連接adaptor上機測序(Illumina Hiseq x Ten Miseq測序儀)。下機后Raw data通過Trimmomatic[14]軟件過濾得到Clean data。

直接提取基因的蛋白序列，用eggNOG-Mapper[15]進行功能注釋后，構建數(shù)據(jù)庫。

利用HISAT2[16]軟件構建參考基因組的索引，將質量控制后的高質量數(shù)據(jù)序列與參考蘋果基因組序列(http://www.rosaceae.org或http://www.ncbi.nlm.nih.gov)進行比對，通過StringTie將比對上的reads進行組裝和定量。

1.3 可變剪接事件鑒定

基于各樣品reads與參考基因組序列的Hisat2比對結果，使用GATK[17]軟件識別測序樣品與參考基因組間的單堿基錯配，識別潛在的SNP位點，并分析這些SNP位點是否影響基因的表達水平或者蛋白產物的種類。變異注釋(SNP、InDel)和預測變異影響是利用SnpEff[18]軟件。根據(jù)變異位點在參考基因組上的位置以及參考基因組上的基因位置信息，獲得變異位點在基因組發(fā)生的區(qū)域(基因間區(qū)、基因區(qū)或CDS區(qū)等)，以及變異產生的影響(同義突變或非同義突變等)。

根據(jù)SNP位點堿基替換方式的不同，可以將SNP位點分為轉換(Transition)和顛換(Transversion)兩種類型[19]。根據(jù)SNP位點的等位基因(Allele)數(shù)目，可以將SNP位點分為純合型SNP位點(只有一個等位基因)和雜合型SNP位點(兩個或多個等位基因)[20]。采用StringTie[21]對Hisat2的比對結果進行拼接，通過ASprofile[22]軟件獲取每個樣品存在的可變剪接類型及相應表達量并將可變剪接類型分為可變5′端或3′端剪接(Alternative exon ends-5′，3′，or both，AE)，內含子滯留(Intron retention，IR)和外顯子跳躍(Skipped exon，SKIP)，第一個外顯子可變剪接(Alternative 5′ first exon-transcription start site，TSS)和最后一個外顯子可變剪接(Alternative 3′ last exon-transcription terminal site，TTS)4種類型。

1.4 發(fā)生可變剪接基因的注釋分析

KEGG富集采用KOBAS[23]軟件進行分析。將5個時期(每個時期3次生物學重復)的共同發(fā)生可變剪接基因注釋到KEGG數(shù)據(jù)庫，分析其代謝通路。

2 結果與分析

2.1 測序質量分析

對‘火焰’5個發(fā)育階段(花后的35、60、95、120、150 d，分別對應S1、S2、S3、S4和S5時期)的果皮進行測序分析。Clean reads的總長度10 430 201～13 583 562，有82.96%到88.64%的reads被映射到參考基因組上，見表1。GC含量均維持在47.10%～48.15%范圍內。每個時期的3個生物學重復高度一致，且所有樣本測序所得的純凈數(shù)據(jù)Q30值均在90%以上(表1)。RNA-seq數(shù)據(jù)質量較好，滿足后續(xù)分析需要。

表1 測序數(shù)據(jù)質量統(tǒng)計Tab.1 Statistical analysis of RNA-seq data

2.2 觀賞海棠‘火焰’果皮可變剪接事件類型鑒定

分析可變剪接事件(圖1)發(fā)現(xiàn)在‘火焰’果皮發(fā)育的5個不同時期均存在可變5′端或3′端剪接，內含子滯留和外顯子跳躍，第一個外顯子可變剪接和最后一個外顯子可變剪接這4種可變剪接類型。

其中，第一個外顯子可變剪接和最后一個外顯子可變剪接占總可變剪接事件的絕大多數(shù)，且二者比例非常接近。在S1期，這兩種類型共占比88.44%；在之后的4個時期中，共占比均高達92%以上。可變5′端或3′端剪接，內含子滯留和外顯子跳躍發(fā)生較少。

在整個發(fā)育過程中，發(fā)生可變剪接的基因數(shù)量(圖1)趨勢較平穩(wěn)，隨著發(fā)育階段的不同略有減少。S1期發(fā)生可變剪接的基因數(shù)量最多為46 952個，S5期最少為43 669個，減少的數(shù)量占比7.19%。

2.3 觀賞海棠‘火焰’果皮可變剪接基因數(shù)量分析

在‘火焰’果皮發(fā)育過程中,每個時期發(fā)生可變剪接的基因總數(shù)不同，各類型在不同時期的基因數(shù)

也不同。通過構建維恩圖(圖2)，共有15 506個基因是5個發(fā)育時期共同擁有的，占總可變剪接基因數(shù)(23 062)的67.24%。每兩個相鄰時期(S1和S2、S2和S3、S3和S4、S4和S5)共有的可變剪接基因數(shù)分別為18 298、18 217、18 326及16 946個。每個發(fā)育階段特有的基因數(shù)依次是432、261、373、270和775個。

2.4 觀賞海棠‘火焰’果皮發(fā)育中共有可變剪接基因的KEGG功能富集

將‘火焰’果皮整個發(fā)育過程中共有的15 506個可變剪接基因進行KEGG功能富集，選擇顯著富集基因最多的20條代謝通路(圖3)。顯著富集基因最多的是嘌呤代謝，其次是氧化磷酸化、糖酵解/糖異生和嘧啶代謝，同時發(fā)現(xiàn)檸檬酸循環(huán)、磷酸戊糖途徑、光合作用等通路富集在果皮發(fā)育過程中。在果皮發(fā)育過程中，可變剪接過程可能主要與NADH、ATP的形成有關，參與呼吸作用、光合作用等多種生理過程。

3 討 論

可變剪接是廣泛存在于生物體內基因表達過程中的一種重要的調控機制，且在植物的各個組織器官及發(fā)育過程中有重要作用[24]。同一個pre-mRNA經過可變剪接可以形成不同mRNA剪接異構體，從而增加蛋白的多樣性[25]。馮雅嵐等[26]借助RNA-seq證明可變剪接主要在轉錄后水平對植物發(fā)育和逆境脅迫響應進行調控，進而在發(fā)育和非生物脅迫響應中發(fā)揮重要作用。孫鴻等[27]發(fā)現(xiàn)小麥中的4個ATG18基因均具有2種可變剪接方式，在抗白粉菌侵染和抗高鹽、干旱、低溫、黑暗和缺氮等逆境中起到重要作用。張盼娃等[28]對干旱脅迫下對玉米蛋白磷酸酶2C基因ZmPP2C26兩個可變剪接體的功能進行分析，證明兩個剪接體均可增加植物對干旱脅迫的敏感性。

可變剪接在植物生長發(fā)育過程中起著十分重要的作用，其在增加植物抗逆性方面被廣泛驗證[5-6,29-33]。可變剪接基因參與蘋果果皮發(fā)育調控鮮有研究。大多對果皮發(fā)育的研究是集中在轉錄因子的調控作用上，MYB-bHLH-WD40復合物、NAC轉錄因子家族、WRKY轉錄因子家族以及部分鋅酯蛋白及乙烯合成相關基因能夠廣泛參與蘋果果皮的發(fā)育調控[5-6,34-37]。蘆筍兩性花發(fā)育過程中的可變剪切的分析證明可變剪切不僅能夠在增加植物抗性中起重要作用，而且其對植物的整個發(fā)育過程非常關鍵[38]。

該研究中利用RNA-seq測序技術分析觀賞海棠‘火焰’果皮在5個不同發(fā)育時期中發(fā)生可變剪接的基因，在不同的發(fā)育時期，主要發(fā)生可變剪切類型是第一個外顯子可變剪接和最后一個外顯子可變剪接，且二者比例非常接近。發(fā)生可變剪接的基因數(shù)量隨著果皮的發(fā)育有所減少，但幅度不大，說明可變剪接在觀賞海棠果皮的整個發(fā)育過程中普遍存在，且可能主要在果實發(fā)育的前期起十分重要的作用。

為了進一步明確這些可變剪接基因的作用，通過KEGG對‘火焰’果皮發(fā)育過程中5個階段共有的可變剪接基因進行功能富集，發(fā)現(xiàn)顯著富集的代謝通路主要是嘌呤代謝，氧化磷酸化，糖酵解/糖異生和嘧啶代謝。氧化磷酸化與多種生物化學代謝過程有著密切關系，尤其是糖酵解、檸檬酸循環(huán)、β氧化等。另外，果糖和甘露糖代謝、半乳糖代謝、磷酸戊糖途徑等代謝通路也被發(fā)現(xiàn)。糖類主要負責提供生物體所需的能量，也可作為信號分子影響糖運輸方向、寄生抗病性等多樣的細胞反應[37]。這些證明可變剪接確實在植物免疫及增加抗性中存在關聯(lián)[36-39]，并且可能借助信號轉導在其他生理過程中發(fā)揮重要作用。