段詩瑤，田 佶，張 杰

(北京農(nóng)學(xué)院植物科學(xué)技術(shù)學(xué)院/農(nóng)業(yè)應(yīng)用新技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102206)

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

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

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

1 材料與方法

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

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

1.2 建庫及分析流程

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

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

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

1.3 可變剪接事件鑒定

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

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

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

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

2 結(jié)果與分析

2.1 測序質(zhì)量分析

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

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

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

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

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

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

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

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

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

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

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

3 討 論

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

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

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

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