江西鉛山紅芽芋和青稈芋的轉(zhuǎn)錄組比較分析

尹明華，曹 晴，陳 紅，鄧思宇，鄧燕梅

(上饒師范學(xué)院 生命科學(xué)學(xué)院，江西 上饒 334001)

江西鉛山紅芽芋種植歷史較為悠久，1 300多年前東晉時期江西鉛山縣就已開始種植紅芽芋。至明代，萬歷《鉛書》記載“齋菜有莧、茄、萵苣、芋、蓮、姜、蘿卜，皆能散穢去積”。更有《鉛山縣志》記錄河口郊區(qū)農(nóng)民于清咸豐五年(1855年)從事紅芽芋生產(chǎn)，遠(yuǎn)銷長江流域等地[1]。江西鉛山紅芽芋因芽紅肉白而名，藥食兼優(yōu)，為江西省名優(yōu)農(nóng)產(chǎn)品，2013年核準(zhǔn)為國家地理標(biāo)志農(nóng)產(chǎn)品[2]。食用，富含淀粉、蛋白質(zhì)和脂肪，鮮美細(xì)膩，粉而不黏，味甘可口，不愧為芋類佳品[3]；藥用，性平味甘，益脾補(bǔ)胃，調(diào)氣解乏，含多種微量元素，可增強(qiáng)免疫、防癌抑瘤，含氟較高，可潔齒防齟護(hù)齒[4]。研究表明，紅芽芋芋稈長而肥大，約占芋頭重量30%，對人體有用的營養(yǎng)物質(zhì)較為豐富，且含有優(yōu)良的食用纖維，但其刺激性澀味較為強(qiáng)烈，容易造成刺舌、刺喉感，芋稈中一水合草酸鈣針晶的特殊晶形及含量與芋稈澀味相關(guān)，因此，紅芽芋芋頭采收后芋稈大多作為廢料處理，一般不用作食用[5]。但江西鉛山縣有種本土種植的青稈芋，地下芋頭不大無收獲價值，而地面部分繁殖快，割去也可快速長出，且芋稈肥大粗壯，食用無澀味，一般用作餐桌的美味佳肴。目前關(guān)于江西鉛山紅芽芋和江西鉛山青稈芋的轉(zhuǎn)錄組分析，尚無相關(guān)研究報道。

隨著第2代測序技術(shù)的發(fā)展，轉(zhuǎn)錄組測序(RNA-Seq)已逐漸成為植物分子生物學(xué)研究的重要手段，廣泛應(yīng)用于功能基因挖掘、分子標(biāo)記開發(fā)、代謝通路和調(diào)控機(jī)制研究等研究[6-7]。轉(zhuǎn)錄組測序能夠獲取植物特定組織或器官在特定時期的基因表達(dá)信息，有利于挖掘品種之間的差異基因[8]。轉(zhuǎn)錄組分析成本低、數(shù)據(jù)量大、效率高、準(zhǔn)確性高，可對組織或者細(xì)胞中所有RNA進(jìn)行測序，并通過讀段(reads)的拼接和豐度統(tǒng)計獲得相應(yīng)的轉(zhuǎn)錄本序列信息及其表達(dá)水平[9]，這些結(jié)果將為園藝植物分子標(biāo)記開發(fā)提供重要序列信息，也為園藝植物種質(zhì)資源遺傳多樣性分析、功能基因挖掘以及分子輔助育種等工作提供了重要基礎(chǔ)[10]。江西鉛山紅芽芋芋稈刺激性澀味較為強(qiáng)烈，而青稈芋芋稈無澀味可作餐桌的美味佳肴，為了解其原因，本研究將以江西鉛山紅芽芋為對象，以江西鉛山青稈芋為對照，利用高通量測序技術(shù)進(jìn)行轉(zhuǎn)錄組基因測序，以期獲得葉片中差異基因的種類、表達(dá)量及差異表達(dá)信息，探究江西鉛山紅芽芋芋稈具刺激性澀味而青稈芋芋稈食用無澀味的相關(guān)基因差異表達(dá)情況，為江西鉛山紅芽芋和江西鉛山青稈芋的品種鑒定和分子育種奠定基礎(chǔ)。

1 材料與方法

1.1 材料

江西鉛山紅芽芋試管苗(HYY組)和江西鉛山青稈芋試管苗(QGY組)。

1.2 方法

1.2.1 RNA提取

參照說明書用Trizol試劑提取江西鉛山紅芽芋(HYY)和江西鉛山青稈芋(QGY)試管苗(全株)總RNA。所提RNA的質(zhì)量由 Nanodrop2000和瓊脂糖凝膠電泳檢測，Agilent2100測定RIN值。單次建庫若RNA總量5 μg，濃度≥200 ng·μL-1，D260/D280介于1.8～2.2即為可用。

1.2.2 Illumina Hiseq測序

RNA測序送由上海美吉生物科技股份有限公司完成。利用帶有Oligo(dT)的磁珠將3′端帶有Poly A尾結(jié)構(gòu)的mRNA從總RNA分離出來。然后利用Truseq TM RNA sample prep kit試劑盒構(gòu)建cDNA文庫。利用Illumina hiseq 4000 進(jìn)行測序和分析。

1.2.3 數(shù)據(jù)分析

使用Trinity軟件(https://github.com/trinityrnaseq/trinityrnaseq/wiki)對clean reads進(jìn)行denovo拼接，形成一個不能在兩端擴(kuò)展的Unigene序列。利用Trinity軟件組裝得到轉(zhuǎn)錄本與參考蛋白質(zhì)數(shù)據(jù)庫NR進(jìn)行BLAST比對。使用軟件Trans Decode對組裝結(jié)果進(jìn)行 CDS 預(yù)測。再將序列注釋到Swiss-Protein、KEGG和GO數(shù)據(jù)庫中獲得基因功能信息。利用RPKM(每百萬個reads映射到外顯子區(qū)的每千個堿基上的reads數(shù))對Unigene的read count 值進(jìn)行標(biāo)準(zhǔn)化。根據(jù)錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR≤0.001)和表達(dá)量倍數(shù)變化(|foldchange|>1.5)的閾值篩選出兩個樣本間差異表達(dá)的基因(DEGs)。對顯著性差異表達(dá)基因進(jìn)行KEGG pathway顯著性富集分析。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)產(chǎn)出統(tǒng)計

經(jīng)過測序質(zhì)量控制，HYY組Raw reads為41 952 830，Raw bases為6 334 877 330，Clean reads為41 298 676，Clean bases為6 147 303 730，Q20和Q30分別為98.79%和96.13%，GC含量為53.09%；QGY組Raw reads為46 252 016，Raw bases為6984054416，Clean reads為45 551 860，Clean bases為6 745 871 513，誤差率為0.0228%，Q20和Q30分別為998.92%和96.45%，GC含量為51.17%。

2.2 轉(zhuǎn)錄組從頭組裝結(jié)果評估

HYY組和QGY組組裝得到的Unigene和Transcript的序列條數(shù)分別為51 941和91 532；組裝得到的所有Unigenes和Transcripts的堿基個數(shù)分別為44 181 196和99 976 031；組裝得到的最長Unigene和Transcript長度均為8 734 bp；組裝得到的最短Unigene和Transcript長度均為201 bp；組裝得到的所有Unigene和Transcript平均長度分別為850.60和1 092.25 bp。

2.3 測序數(shù)據(jù)與組裝結(jié)果比對

HYY組Clean reads為41 298 676，Mapped reads為31 959 832，Mapped ratio為77.39%。QGY組Clean reads為45 551 860，Mapped reads為30 904 580，Mapped ratio為67.84%。

2.4 表達(dá)量統(tǒng)計

圖1展示各樣本表達(dá)Unigene/Transcript的表達(dá)量分布整體情況和表達(dá)量密度分布情況。HYY組和QGY組表達(dá)量密度在0～0.7。

圖1 HYY組和QGY組表達(dá)量分布盒形圖(A)和密度圖(B)Fig.1 Box diagram (A) and density diagram (B) of expression quantity distribution in HYY group and QGY group

2.5 表達(dá)量樣本間Venn分析

從圖2可知，HYY組和QGY組表達(dá)的共有基因數(shù)為14 038，HYY組單獨表達(dá)的基因數(shù)為17 843，QGY組單獨表達(dá)的基因數(shù)為18 634。

圖2 HYY組和QGY組表達(dá)量樣本間Venn圖Fig.2 Venn map between HYY group and QGY group

2.6 表達(dá)量樣本間相關(guān)性分析

由圖3可知，HYY組和QGY組表達(dá)量的相關(guān)系數(shù)為0.035，表明基因在HYY組和QGY組2個樣本間的表達(dá)量相似度越低，即樣本間相關(guān)性不夠。

圖3 HYY組和QGY組表達(dá)量樣本間相關(guān)性熱圖Fig.3 Thermogram of correlation between HYY group and QGY group

2.7 差異表達(dá)分析

HYY組和QGY組共產(chǎn)生差異表達(dá)基因(differentially expressed gene，DEG) 32 555個，與QGY組比較，HYY組上調(diào)基因15 887個，下調(diào)基因16 668個(圖4)。從表1可知，HYY組和QGY組前10個差異表達(dá)基因為假設(shè)蛋白GLYMA_13G018000、類絲氨酸/蘇氨酸蛋白激酶D6PK、植物/T3A5-60蛋白質(zhì)、非特征蛋白LOC110018203亞型X3、過硫雙加氧酶ETHE1同系物(線粒體樣)、肉桂酸-4-羥化酶、生長素應(yīng)答蛋白IAA17、磷脂酰肌醇4-激酶γ7、非特征蛋白LOC105044622亞型X2。

圖4 HYY組和QGY組表達(dá)量差異火山圖(A)和表達(dá)量差異散點圖(B)Fig.4 Volcano map (A) and scatter map (B) of expression difference between HYY group and QGY group

表1 HYY組和QGY組表達(dá)量差異統(tǒng)計表(前10)Table 1 Statistical table of expression difference between hyy group and QGY group (top10)

2.8 差異表達(dá)基因GO富集分析

HYY組和QGY組差異表達(dá)基因GO富集分析結(jié)果如圖5和表2所示。從圖5和表2可知，通過富集，多糖代謝過程、基因沉默、果膠分解代謝過程、碳水化合物代謝過程、生長素激活的信號通路、細(xì)胞外區(qū)、膜的固有成分、膜的組成成分、膜部件、凋亡細(xì)胞、作用于糖基鍵水解酶活性、水解O-糖基化合物的水解酶活性、蛋白質(zhì)二聚活性、DNA結(jié)合、氧化還原酶活性對二酚及其相關(guān)物質(zhì)的作用等GO term富集顯著。

表2 部分差異表達(dá)基因GO富集結(jié)果Table 2 GO enrichment results of some differentially expressed genes

圖5 HYY組和QGY組差異表達(dá)基因GO富集分析柱形圖Fig.5 Column chart of GO enrichment analysis of differentially expressed genes in HYY group and QGY group

2.9 差異表達(dá)基因KEGG注釋

從圖6和表3可知，HYY組和QGY組差異表達(dá)基因富集的KEGG通路主要有植物激素信號轉(zhuǎn)導(dǎo)、植物病原相互作用、半胱氨酸與蛋氨酸代謝、苯丙酸生物合成、皮質(zhì)醇合成與分泌、過氧化物酶體、植物MAPK(絲裂原活化蛋白激酶)信號途徑、不匹配修復(fù)、淀粉和蔗糖代謝、硫代謝、戊糖和葡萄糖醛酸相互轉(zhuǎn)化、C型凝集素受體信號通路、甘油酯代謝、不飽和脂肪酸的生物合成、α-亞麻酸代謝、二萜生物合成、類黃酮生物合成、脂肪酸延長、光合作用、吲哚生物堿生物合成等。其中植物激素信號轉(zhuǎn)導(dǎo)、植物病原相互作用、植物MAPK信號途徑、淀粉和蔗糖代謝、苯丙酸生物合成、半胱氨酸與蛋氨酸代謝、過氧化物酶體注釋的基因數(shù)均超過80,這可能是導(dǎo)致參試品種差異的主要原因。

圖6 差異表達(dá)基因KEGG富集分析柱形圖(A)和氣泡圖(B)Fig.6 Column chart (A) and bubble chart (B) of enrichment analysis of differential expression gene KEGG

表3 差異表達(dá)基因KEGG富集部分分析Table 3 Analysis of enrichment of differential expression gene KEGG

2.10 與芋稈澀味相關(guān)差異表達(dá)基因的篩選

芋稈中所含的草酸鈣針晶為芋稈主要刺激性澀味成分之一。由表4可知，與芋稈澀味相關(guān)差異表達(dá)基因如線粒體樣鈣單轉(zhuǎn)運(yùn)蛋白2(calcium uniporter protein 2，mitochondrial-like)基因上調(diào)，而L-抗壞血酸氧化酶(L-ascorbate oxidase)基因下調(diào)，導(dǎo)致形成草酸鈣針晶可能性下降。

表4 與芋稈澀味相關(guān)差異表達(dá)基因Table 4 Differentially expressed genes related to astringency of taro

3 討論

轉(zhuǎn)錄組一般指所有表達(dá)基因的身份以及其轉(zhuǎn)錄水平，廣義上指細(xì)胞內(nèi)在某一生理條件下的所有轉(zhuǎn)錄產(chǎn)物集合(mRNA、核糖體RNA、轉(zhuǎn)運(yùn)RNA及非編碼RNA)；狹義上特指所有mRNA的集合[11]。轉(zhuǎn)錄組測序有助于新基因的發(fā)現(xiàn)、基因功能注釋，基因差異表達(dá)和分子標(biāo)記篩選[12]。通過第二代測序——Illumina測序技術(shù)，能夠全面快速地獲得植物組織在某一狀態(tài)下幾乎所有轉(zhuǎn)錄本序列信息，已經(jīng)廣泛用于模式和非模式植物的研究[13]。

KEGG數(shù)據(jù)庫是全面研究基因產(chǎn)物在細(xì)胞中的代謝途徑和研究基因產(chǎn)物的功能的數(shù)據(jù)庫，基于KEGG的代謝途徑的分析有助于更進(jìn)一步了解基因的生物學(xué)功能[14]。紫花苜蓿差異基因的GO功能分析和KEGG途徑分析表明，核糖體、剪接體、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工和碳代謝途徑可能導(dǎo)致2個品種紫花苜蓿的蛋白質(zhì)含量差異，而植物激素信號轉(zhuǎn)導(dǎo)代謝途徑可能是導(dǎo)致2個紫花苜蓿品種中性洗滌纖維和酸性洗滌纖維含量差異的重要原因之一[15]。黃果、滿天星和紫果3個品種百香果的KEGG富集分析表明，黃果與滿天星的細(xì)胞胞吞作用有顯著差異，黃果與紫果的玉米素生物合成和糖酵解途徑差異基因富集顯著，滿天星與紫果的纈氨酸、亮氨酸和異亮氨酸生物合成中存在較大差異[16]。次生代謝物生物合成和萜類和聚酮化合物代謝途徑可能與4個茶樹品種(藪北、湘妃翠、牛皮茶和鐵香茗)的苯甲醇合成有關(guān)，編碼茶樹苯甲醇合成途徑可能有5條關(guān)鍵酶β-葡萄糖苷酶基因[17]。在本實驗中，通過對江西鉛山紅芽芋和江西鉛山青稈芋的轉(zhuǎn)錄組分析，結(jié)果表明，HYY組和QGY組差異表達(dá)基因富集的KEGG通路主要有植物激素信號轉(zhuǎn)導(dǎo)、植物病原相互作用、半胱氨酸與蛋氨酸代謝、苯丙酸生物合成、皮質(zhì)醇合成與分泌、過氧化物酶體、植物MAPK信號途徑、不匹配修復(fù)、淀粉和蔗糖代謝、硫代謝、戊糖和葡萄糖醛酸相互轉(zhuǎn)化、C型凝集素受體信號通路、甘油酯代謝、不飽和脂肪酸的生物合成、α-亞麻酸代謝、二萜生物合成、類黃酮生物合成、脂肪酸延長、光合作用、吲哚生物堿生物合成等。其中植物激素信號轉(zhuǎn)導(dǎo)(175個基因)、植物病原相互作用(151個基因)、植物MAPK(絲裂原活化蛋白激酶)信號途徑(118個基因)、淀粉和蔗糖代謝(118個基因)、苯丙酸生物合成(95個基因)、半胱氨酸與蛋氨酸代謝(84個基因)、過氧化物酶體(81個基因)注釋的基因數(shù)較多，這可能是導(dǎo)致參試品種差異的主要原因。

植物激素實際上是一類信號小分子，植物的生長發(fā)育需要通過植物激素信號轉(zhuǎn)導(dǎo)系統(tǒng)進(jìn)行協(xié)調(diào)。一些因素包括內(nèi)因或外因可能會誘導(dǎo)一系列植物激素基因的表達(dá)，繼而作用于激素受體或組件，最終顯現(xiàn)出不同性狀[18]。江西鉛山紅芽芋和江西鉛山青稈芋可能由于其植物激素(生長素、細(xì)胞分裂素、赤霉素、脫落酸、乙烯、油菜素內(nèi)酯、茉莉酸和水楊酸等)在這些激素信號轉(zhuǎn)導(dǎo)系統(tǒng)中有的受體或關(guān)鍵組分因互作(interact)或交互應(yīng)答(crosstalk)會產(chǎn)生協(xié)同或拮抗的作用而使信號途徑網(wǎng)絡(luò)化，從而對江西鉛山紅芽芋和江西鉛山青稈芋的生長發(fā)育、防御、環(huán)境適應(yīng)等產(chǎn)生作用[19]。植物病原真菌是造成植物病害的主要病原菌之一，也是影響園藝植物產(chǎn)量和品質(zhì)最主要因素之一。研究表明，植物細(xì)胞信號調(diào)控的普遍機(jī)制之一是蛋白磷酸化修飾，在植物—病原微生物互作過程中，關(guān)鍵調(diào)控蛋白的磷酸化狀態(tài)直接影響到免疫信號的激活[20]。在江西鉛山紅芽芋和江西鉛山青稈芋的植物—病原微生物互作過程中，可能依賴于蛋白激酶和蛋白磷酸酶來調(diào)節(jié)多種蛋白的磷酸化狀態(tài)最終調(diào)控江西鉛山紅芽芋和江西鉛山青稈芋的免疫途徑，并且病原微生物效應(yīng)蛋白可能還可改變植物代謝以滿足侵染病原微生物的營養(yǎng)需要，這可能也是影響江西鉛山紅芽芋和江西鉛山青稈芋品質(zhì)差異的一個原因。MAPK信號轉(zhuǎn)導(dǎo)通路由高度保守的絲/蘇氨酸蛋白激酶MAPKKK、MAPKK和MAPK組成，廣泛存在于真核生物體內(nèi)[21]。近年來研究表明，在植物中，該信號轉(zhuǎn)導(dǎo)通路通過其3個組分將胞外信號經(jīng)細(xì)胞膜上的受體傳遞至下游應(yīng)答因子，在個體生長發(fā)育及防御反應(yīng)中發(fā)揮重要作用[22]。江西鉛山紅芽芋和江西鉛山青稈芋的生長發(fā)育需要各細(xì)胞、組織和器官間的精準(zhǔn)協(xié)調(diào)，蛋白質(zhì)磷酸化可能參與了所有細(xì)胞信號轉(zhuǎn)導(dǎo)機(jī)制，參與了細(xì)胞間以及細(xì)胞與外界信號間的交流，MAPK信號轉(zhuǎn)導(dǎo)通路在江西鉛山紅芽芋和江西鉛山青稈芋生長發(fā)育及器官發(fā)生中可能起到了核心作用。碳水化合物是植物光合作用的最終產(chǎn)物，碳水化合物含量高低體現(xiàn)了光合作用效率高低[23]。蔗糖和淀粉是光合作用同化產(chǎn)物的主要分配形式，蔗糖是光合產(chǎn)物的主要運(yùn)輸和貯存形式，為植物的生長發(fā)育提供能量，而淀粉是一種暫貯形式，供植物夜晚呼吸消耗[24]。淀粉和蔗糖的積累和轉(zhuǎn)換不僅與物質(zhì)基礎(chǔ)的形成有關(guān)，而且也與產(chǎn)量和品質(zhì)密切相關(guān)[25]。江西鉛山紅芽芋和江西鉛山青稈芋可能通過增加體內(nèi)蔗糖—淀粉代謝關(guān)鍵酶基因的表達(dá)，進(jìn)而形成兩者的品種特性和獨特品質(zhì)。苯丙酸類代謝是抗菌途徑之一，通過苯丙酸類代謝能合成多種次生酚類物質(zhì)，多屬于植保素類物質(zhì)(如酚類化合物、香豆酸酯類、類黃酮、木質(zhì)素等)能抑制病原菌的生長，對植物抗病具有重要作用[26]。江西鉛山紅芽芋和江西鉛山青稈芋的差異可能也在于苯丙氨酸解氨酶活性的大小以及木質(zhì)素和綠原酸含量的高低。研究表明，植物蛋白普遍存在氨基酸不均衡的現(xiàn)象，有些植物缺乏蛋氨酸和半胱氨酸等含硫氨基酸[27]。蛋氨酸和半胱氨酸參與合成含硫維生素(如維生素 B1、生物素和硫辛酸等)、含硫輔酶(如輔酶A)、形成谷胱甘肽以及許多含硫防御性有機(jī)物(如金屬硫蛋白、植物鰲合肽、蒜素和芥子油苷等)和含硫信號分子(如結(jié)瘤因子和膨脹素)等在新陳代謝上重要的有機(jī)化合物[28]。江西鉛山紅芽芋和江西鉛山青稈芋可能通過蛋氨酸和半胱氨酸代謝廣泛參與種芋萌發(fā)、幼苗發(fā)育、 脅迫響應(yīng)和組織分化衰老等過程。過氧化物酶體是高度動態(tài)、代謝活躍的細(xì)胞器，主要參與脂肪酸等脂質(zhì)的代謝及產(chǎn)生和清除不同的活性氧(ROS)[29]。因此，過氧化物酶體的動態(tài)平衡由ROS、過氧化物酶體蛋白酶及自噬過程調(diào)節(jié)，對于維持細(xì)胞的氧化還原平衡至關(guān)重要[30]。近年來，越來越多的證據(jù)表明，過氧化物酶體是ROS的重要來源，包含一套復(fù)雜的抗氧化防御機(jī)制來調(diào)控ROS的積累，ROS作為植物激素信號的第二信使，可以調(diào)節(jié)植物發(fā)育過程和應(yīng)答反應(yīng)[31]。過氧化物酶體、線粒體和葉綠體與ROS依賴的信號網(wǎng)絡(luò)之間的相互作用是充分發(fā)揮江西鉛山紅芽芋和江西鉛山青稈芋過氧化物酶體功能的關(guān)鍵。隨著生物信息學(xué)發(fā)展與更新，有望開發(fā)出對江西鉛山紅芽芋和江西鉛山青稈芋過氧化物酶體蛋白功能的預(yù)測軟件，及通過信息技術(shù)與分子技術(shù)的結(jié)合可以多層面對ROS介導(dǎo)的途徑如何與江西鉛山紅芽芋和江西鉛山青稈芋生長發(fā)育、自噬過程及脅迫應(yīng)答信號網(wǎng)絡(luò)的整合有更深認(rèn)識。

芋稈中所含的草酸鈣針晶為芋稈主要刺激性澀味成分之一[5]。與芋稈澀味相關(guān)差異表達(dá)基因主要有線粒體樣鈣單轉(zhuǎn)運(yùn)蛋白2(calcium uniporter protein 2，mitochondrial-like)基因和L-抗壞血酸氧化酶(L-ascorbate oxidase)基因[32]。江西鉛山青稈芋通過上調(diào)線粒體樣鈣單轉(zhuǎn)運(yùn)蛋白2(calcium uniporter protein 2，mitochondrial-like)基因活性，造成胞質(zhì)鈣離子濃度下降，并再通過下調(diào)L-抗壞血酸氧化酶(L-ascorbate oxidase)基因活性，這樣以抗壞血酸作為前體物質(zhì)通過氧化降解后生成草酸的可能也下降，最終導(dǎo)致江西鉛山青稈芋形成草酸鈣針晶可能性下降。這也是江西鉛山紅芽芋芋稈具刺激性澀味而青稈芋芋稈食用無澀味的根本原因。