堿脅迫下榆葉梅差異基因表達(dá)分析

劉 佳 田 云 鄧家林

(1四川省農(nóng)業(yè)科學(xué)院園藝研究所,四川 成都 610066；2農(nóng)業(yè)農(nóng)村部西南地區(qū)園藝作物生物學(xué)與種質(zhì)創(chuàng)制重點(diǎn)實(shí)驗(yàn)室,四川 成都 610066)

我國(guó)川中丘陵地區(qū)為中國(guó)李(Prunus salicinaLindl.)的主要種植區(qū),此區(qū)域的土壤廣泛由石灰性紫色土組成,pH 值大多介于7.69～8.47 之間[1]。堿性土壤上生長(zhǎng)的果樹容易出現(xiàn)缺鐵性黃化癥狀[2],使其生長(zhǎng)受抑制,影響開花或坐果,甚至導(dǎo)致樹體死亡[3]。近年來(lái),學(xué)者們對(duì)植物如何在鹽堿土地上正常生長(zhǎng)進(jìn)行了大量研究[4-5],發(fā)現(xiàn)選擇適宜的耐鹽堿植物是高效利用鹽堿土壤的主要方法[5]。果樹大多通過嫁接進(jìn)行繁殖,所以果樹的耐堿性實(shí)際取決于砧木的耐堿性。此外,砧木的品種也直接影響著果樹壽命的長(zhǎng)短、結(jié)果實(shí)的早晚以及產(chǎn)量的高低,因此篩選適宜的耐堿砧木是增強(qiáng)果樹耐堿性的有效途徑[6]。

榆葉梅(Prunus trilobaLindl.),又稱小桃紅或“鸞枝”,是薔薇科(Rosaceae)李屬(PrunusL.)落葉觀花灌木或小喬木,在我國(guó)的東北、華北和西北地區(qū)廣泛栽植[7-8]。榆葉梅具有很多品種,在我國(guó)已有300 多年的栽培歷史[7,9]。它是一種具有2 n＝8 X＝64 核型的八倍體植物,具有適應(yīng)性極強(qiáng)、耐貧瘠和低水平管理、耐寒、耐旱以及抗病性較強(qiáng)等特點(diǎn)[7,10-11]。此外,榆葉梅還是一種天然的耐鹽堿性植物,在鹽堿土壤(pH 值8.8,鹽含量0.3%)中生長(zhǎng)良好,與中國(guó)李嫁接后表現(xiàn)出親和性好、果實(shí)早熟、產(chǎn)量高等特點(diǎn),是川中丘陵地區(qū)李樹的良好砧木[10]。

土壤鹽化、堿化或鹽堿化是植物常遭遇的幾種非生物逆境脅迫[12-13]。研究表明,堿脅迫對(duì)植物造成的傷害遠(yuǎn)大于鹽脅迫[13-14]。目前有關(guān)植物對(duì)鹽堿脅迫響應(yīng)的報(bào)道主要集中于鹽脅迫響應(yīng)機(jī)制[15],而關(guān)于植物對(duì)堿脅迫響應(yīng)的報(bào)道較少,急需探明植物響應(yīng)堿脅迫的各種潛在機(jī)制(如細(xì)胞生理生化與分子抗逆機(jī)制等)。近年來(lái),多種大規(guī)模平行測(cè)序平臺(tái)的飛速發(fā)展已成為探明植物響應(yīng)堿脅迫分子機(jī)制的重要手段。其中以RNA-Seq 技術(shù)為基礎(chǔ)的轉(zhuǎn)錄組測(cè)序,因具有耗時(shí)少、高通量、高靈敏度、操作簡(jiǎn)單、低成本等優(yōu)點(diǎn),在眾多領(lǐng)域中已被大量推廣運(yùn)用。且該技術(shù)能夠快速有效地對(duì)植物某一特定過程或階段處理后的基因表達(dá)情況進(jìn)行診斷,發(fā)掘參與脅迫響應(yīng)的功能基因并預(yù)測(cè)基因的功能,使其成為目前常用于研究植物抗逆作用機(jī)制的強(qiáng)大手段[15-16]。RNA-Seq 技術(shù)已被應(yīng)用于研究多種植物的抗堿響應(yīng)機(jī)制[17-20]。

本研究采用Illumina HiSeqTM2500 高通量測(cè)序技術(shù),對(duì)堿脅迫0、1、3、6、12 h 下的榆葉梅葉片進(jìn)行de novo轉(zhuǎn)錄組測(cè)序,并對(duì)得到的大量轉(zhuǎn)錄本進(jìn)行從頭組裝,獲得大量unigenes 差異表達(dá)(DEGs),對(duì)其進(jìn)行基因GO 和KEGG 富集分析,同時(shí)對(duì)7 個(gè)堿性脅迫反應(yīng)相關(guān)基因進(jìn)行RT-qPCR 驗(yàn)證,以期進(jìn)一步了解榆葉梅堿脅迫分子的響應(yīng)機(jī)制,同時(shí)為研究以榆葉梅為砧木的中國(guó)李在堿脅迫下的適應(yīng)機(jī)制奠定分子研究基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料及堿脅迫處理

2014年11月將榆葉梅種子(來(lái)源于新疆伊犁植物園)于沙中層積保存,2015年2月中旬,將經(jīng)沙藏萌動(dòng)的榆葉梅種子播種于底部直徑25 cm 的帶孔瓦盆。每個(gè)瓦盆中裝入等量的培養(yǎng)土,培養(yǎng)土由泥炭土與草炭土按同等比例混合而成。待榆葉梅幼苗的高度長(zhǎng)至約15 cm 后,選擇大小、長(zhǎng)勢(shì)基本一致的150 株盆栽幼苗作為試驗(yàn)材料,對(duì)其進(jìn)行堿脅迫處理。

堿脅迫處理的植株用堿液澆灌,堿液由9∶1摩爾比的NaHCO3和Na2CO3溶液混合構(gòu)成,其濃度和pH值分別為50 mmol·L-1和9.11±0.104 ,堿液處理之前土壤pH 值為7.4,處理后土壤pH 值為9.0。處理組材料于晴天下午16：00 用該堿液澆灌土壤,直至缽底流出堿液開始計(jì)時(shí),對(duì)照同期使用清水澆灌(澆灌后土壤pH 值為7.1),用透明塑料薄膜對(duì)這些材料進(jìn)行避雨遮擋[21]。

1.2 RNA 提取、cDNA 文庫(kù)構(gòu)建及Illumina 測(cè)序

在堿脅迫處理后的0、1、3、6、12 h 分別對(duì)榆葉梅中、上部功能葉片(第2～第3 節(jié))進(jìn)行采樣。每次采集的樣品存放于液氮中,以備后續(xù)RNA 的提取。其中每個(gè)時(shí)間點(diǎn)進(jìn)行2 個(gè)生物學(xué)重復(fù),每個(gè)重復(fù)由10 株材料的葉片混合而成。本試驗(yàn)共10 個(gè)樣本材料?？俁NA 提取采用北京艾德萊(aidlab)生物科技有限公司的通用植物RNA 提取試劑盒。分別使用Nanodrop 2000 UV-Vis 分光光度計(jì)(Thermo Fisher Scientific 公司,美國(guó))、Qubit 2.0 熒光計(jì)(Life Technologies,美國(guó))和Agilent 2100 生物分析儀(安捷倫科技有限公司,美國(guó))檢測(cè)RNA 樣品的濃度、純度和完整性。對(duì)符合試驗(yàn)?zāi)康牡腞NA 樣品用于文庫(kù)制備。

提取的10 個(gè)供試榆葉梅RNA 樣本送至蘇州金唯智(GENEWIZ)公司進(jìn)行cDNA 文庫(kù)的構(gòu)建和轉(zhuǎn)錄組測(cè)序[在Illumina HiSeq 2500 平臺(tái)上采用雙末端(PE)測(cè)序法進(jìn)行高通量測(cè)序]。

1.3 De novo 轉(zhuǎn)錄組組裝

對(duì)獲得的FASTq 格式的大量原始測(cè)序數(shù)據(jù)(raw data)采用數(shù)據(jù)質(zhì)量統(tǒng)計(jì)軟件Trimmomatic(v0.30)[22]進(jìn)行處理,去除其中的接頭序列以及低質(zhì)量reads(3′或5′末端質(zhì)量值低于20 或者含N 的堿基；短于75 bp的序列),以此獲得用于后續(xù)信息分析的clean data。使用fastqc(http：//www.bioinformatics.babraham.ac.uk/projects/fastqc/)軟件對(duì)供試的10 個(gè)cDNA 文庫(kù)的clean data 進(jìn)行質(zhì)量評(píng)估,評(píng)估內(nèi)容：Q20、Q30、GC 和N含量。隨后,采用Trinity 軟件(版本r2013-02-25)[23]實(shí)現(xiàn)clean data 的重新組裝,再使用序列聚類軟件TGICL[24]做進(jìn)一步序列拼接和冗余處理,以產(chǎn)生更長(zhǎng)、更完整的非冗余序列,稱為unigenes。最后,使用BLAST(http：//blast.ncbi.nlm.nih.gov/Blast.cgi)(e 值＜1e-5)將unigenes 比對(duì)到NCBI 非冗余蛋白數(shù)據(jù)庫(kù)Nr[National Center for Biotechnology Information(NCBI)non-redundant protein sequences]中,以確定unigene 的序列方向并得到最佳的功能注釋。

1.4 基因表達(dá)定量與差異表達(dá)分析

每個(gè)樣品獲得的unigenes 均采用RSEM 軟件(Ⅴ1.2.4)[25]中 的 FPKM(fragments per kilobase per million reads)方法對(duì)其基因表達(dá)量水平進(jìn)行估算。采用基于負(fù)二項(xiàng)分布模型的Bioconductor 軟件包中的DESeq(Ⅴ1.14.0)[26],對(duì)不同堿脅迫處理下樣品在基因表達(dá)水平分析中得到的expected_count 數(shù)據(jù)進(jìn)行處理,以此獲得在堿脅迫下榆葉梅的DEGs。P值用于基因差異表達(dá)檢驗(yàn)。使用FDR(False discovery rate)方法確定P值的閾值。按照差異顯著性標(biāo)準(zhǔn)[FDR 閾值≤0.05,且｜log2(fold change)｜＞1]篩選顯著性差異表達(dá)的基因。隨后,采用基于Walleius non-central hyper-geometric distribution 的R package GOseq 對(duì)差異表達(dá)顯著的基因進(jìn)行GO 功能顯著性富集分析,P值閾值設(shè)置為0.01[27]。此外,對(duì)差異表達(dá)顯著的基因采用KOBAS 軟件[28]進(jìn)行KEGG Pathway 富集分析,以此確定顯著性DEGs 參與的最主要信號(hào)轉(zhuǎn)導(dǎo)途徑和生化代謝途徑。

1.5 RT-qPCR 驗(yàn)證

為了驗(yàn)證RNA-Seq 檢測(cè)到的DEGs,對(duì)其中7 個(gè)與堿脅迫相關(guān)的差異表達(dá)顯著基因進(jìn)行了RT-qPCR分析。用1.2 中的RNA 提取方法對(duì)堿脅迫后0、3、6、12 h 的榆葉梅葉片進(jìn)行RNA 提取。使用北京Aidlab公司反轉(zhuǎn)錄試劑盒(TUREscript 1st Stand cDNA SYNTHESIS Kit)進(jìn)行cDNA 的合成,采用40 μL 反應(yīng)體系,包括500 ng 總RNA、8 μL 5×RT Reaction Mix、2 μL Oligo(dT)引物、2 μL TUREscript H- RTase/RI Mix,RNase Free dH2O 補(bǔ)足至40 μL。反轉(zhuǎn)錄反應(yīng)程序：42℃、60 min,65℃、10 min,待反應(yīng)結(jié)束后,得到cDNA,于-20℃保存。特異性引物情況詳見表1。采用qTOWER2.2 實(shí)時(shí)PCR 系統(tǒng)(DBI? Bioscience,德國(guó))進(jìn)行RT-qPCR,其中所有反應(yīng)均具有3 個(gè)生物學(xué)重復(fù)和技術(shù)重復(fù)。RT-qPCR 反應(yīng)體系最終體積為10 μL：5 μL 1×2×SYBR ? Green Supermix(DBI ?Bioscience,德國(guó))、0.5 μL 200 nmol·L-1正向引物、0.5 μL 200 nmol·L-1反向引物、1 μL cDNA 模板、3 μL ddH2O。RT-qPCR 反應(yīng)程序：95℃預(yù)變性3 min；95℃變性10 s；60℃擴(kuò)增30 s,循環(huán)39 次；溶解曲線分析,即隨后將溫度從60℃逐漸升至95℃,每次加1℃,每加1℃停留4 s。以ACTIN 基因作為內(nèi)參基因,并通過2-△△Ct法進(jìn)行數(shù)據(jù)分析。

表1 RT-qPCR 所用基因引物信息Table 1 Primer information for RT-qPCR

2 結(jié)果與分析

2.1 測(cè)序結(jié)果、組裝與注釋

采用Illumina HiSeq 2500 高通量測(cè)序技術(shù)對(duì)榆葉梅堿脅迫處理不同時(shí)間下的葉片進(jìn)行轉(zhuǎn)錄組測(cè)序,共獲得了595 894 216 條原始序列數(shù)據(jù)(raw reads)。將這些raw data 提交到NCBI 的SRA 數(shù)據(jù)庫(kù),其登記號(hào)為SRP093875。經(jīng)過嚴(yán)格的質(zhì)量評(píng)估和數(shù)據(jù)過濾后,共得到529 847 752 條高質(zhì)量clean reads(約62.10 Gb)。每個(gè)樣本clean reads 的Q20、Q30 和GC 含量分別都大于99.16%、97.69%和45.16%(表2)。結(jié)果表明獲得的clean reads 具有足夠高的質(zhì)量用于后續(xù)研究。使用Trinity 和TGICL 軟件對(duì)clean reads 進(jìn)行從頭組裝,共獲得了124 786 個(gè)unigenes,其平均長(zhǎng)度為668.65 bp,N50 長(zhǎng)度為1 111 bp。在組裝的unigenes中,長(zhǎng)度為200～500 bp 的unigenes 最多,有80 060個(gè),占unigenes 總數(shù)的64.16%；其次為500～1 000 bp的unigenes,為22 426 個(gè),占總數(shù)的17.97%(圖1)。

本研究將獲得的全部unigenes 與Nr 數(shù)據(jù)庫(kù)進(jìn)行了比對(duì)注釋。結(jié)果顯示,有52 691條 unigenes(42.23%)被成功注釋功能。其中對(duì)注釋到的近源物種進(jìn)行統(tǒng)計(jì),發(fā)現(xiàn)占總體(52 691 條unigenes)45.41%、35.30%、2.56%、1.94%和1.61%的unigenes分別匹配到梅(Prunus mume)、碧桃(Prunus persica)、蘋果(Malus domestica)、川桑(Morus notabilis)和鴨梨(Pyrus×bretschneideri)中,表明榆葉梅與這幾類植物具有較近的親緣關(guān)系(表3)。此外,這些匹配到的物種除川桑外,均屬于薔薇科,也表明本研究獲得的序列被正確地組裝和注釋。

表2 序列分析統(tǒng)計(jì)結(jié)果Table 2 Summary of sequences analysis

表3 被注釋最多unigenes 的5 個(gè)近源物種Table 3 The five relative species annotated with the most unigenes

2.2 基因表達(dá)定量分析與DEGs 的識(shí)別

采用FPKM 法估計(jì)了每個(gè)樣品的基因表達(dá)水平。在這10 個(gè)樣本中,FPKM 值位于0～1 之間的unigenes占總體的20.27%～32.49%,而FPKM 值超過60 的unigenes 僅占總體的1.24%～1.57%(表4)。

堿脅迫下的榆葉梅葉片中有8 948 條unigenes(占總體的7.17%)被鑒定為顯著性DEGs,其中0、186、3 211 和8 347 個(gè)顯著性DEGs 分別于1 h vs.0 h、3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 比較中識(shí)別(圖2),12 h vs.0 h 比較中篩選出最多的顯著性DEGs。在8 948個(gè)顯著性DEGs 中,3 h vs.0 h、6 h vs.0 h、12 h vs.0 h中分別有181、1 999、4 895 個(gè)DEGs 上調(diào)；3 h vs.0 h、6 h vs.0 h、12 h vs.0 h 中分別有5、1 212 和3 452 個(gè)DEGs 下調(diào)(圖2)。

圖1 組裝的unigenes 的長(zhǎng)度分布圖Fig.1 Pie chart for length distribution of the assembled unigenes

2.3 DEGs 功能分類

2.3.1 GO 富集分析 為了進(jìn)一步描述連續(xù)堿脅迫時(shí)間下的基因表達(dá)變化,對(duì)3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 中得到的顯著性DEGs 進(jìn)行了GO 功能顯著性富集分析(P值≤0.01),以期找到與整個(gè)基因組背景相比,在DEGs 中顯著富集的GO 條目。分析顯示,GO富集包含3 個(gè)ontology,分別為細(xì)胞組分、分子功能以及生物過程。在3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h中共分別識(shí)別了186、3 211 和8 347 個(gè)顯著性DEGs,其中分別有0、24 和72 個(gè)顯著性DEGs 被注釋到GO細(xì)胞組分中；分別有9、146 和307 個(gè)DEGs 被注釋到GO 分子功能中；分別有11、86 和184 個(gè)顯著性DEGs在GO 生物過程中注釋功能。在細(xì)胞組分類型中(圖3),3 h vs.0 h 中沒有顯著性DEGs 注釋功能,而在6 h vs.0 h 和12 h vs.0 h 中,細(xì)胞[16(66.67%)；53(73.61%)]、細(xì)胞成分[16(66.67%)；53(73.61%)]、細(xì)胞內(nèi)部分[13(54.17%)；47(65.28%)]和細(xì)胞內(nèi)[13(54.17%)；47(65.28%)]條目中所注釋到的DEGs 數(shù)量較多。在細(xì)胞分子功能類型(圖4)中,核酸結(jié)合轉(zhuǎn)錄因子活性[4,(44.44%)]和序列特異性DNA結(jié)合轉(zhuǎn)錄因子活性[4,(44.44%)]是3 h vs.0 h 中注釋最多DEGs 的兩個(gè)GO 條目；對(duì)于6 h vs.0 h 和12 h vs.0 h,注釋DEGs 較多的條目有結(jié)合活性[70(47.95%)；152(65.28%)]、催化活性[51(34.93%)；106(34.53%)]、雜環(huán)化合物結(jié)合[38(26.03%)、79(27.73%)]以及有機(jī)環(huán)狀化合物結(jié)合[38(26.03%)；79(27.73%)]。在生物學(xué)過程功能類型(圖5)中,3 h vs.0 h 中注釋最多DEGs 的GO 條目為單生物代謝過程[5,(45.45%)]、氧脂代謝過程[2,(18.18%)]以及氧脂生物合成過程(oxylipin biosynthetic process)[2,(18.18%)]；而在6 h vs.0 h 和12 h vs.0 h 中,代謝過程[49(56.98%)；101(54.89%)]、細(xì)胞進(jìn)程[28(32.56%)；78(42.39%)]、有機(jī)物質(zhì)代謝過程[29(33.72%)；70(38.04%)]和初級(jí)代謝過程[27(31.40%)；63(34.24%)]所富集的DEGs 最多?？梢?隨著堿脅迫處理時(shí)間的延長(zhǎng),更多的DEGs 被注釋到功能,其中細(xì)胞組分類型中的“細(xì)胞”和“細(xì)胞成分”,分子功能類型中的“結(jié)合活性”和“催化活性”,以及生物學(xué)過程中的“代謝過程”和“細(xì)胞進(jìn)程”注釋到的DEGs 最多,表明榆葉梅的內(nèi)部組織和器官在堿脅迫下經(jīng)歷了多種合成代謝過程,以此來(lái)抵御和適應(yīng)堿脅迫環(huán)境。

表4 不同表達(dá)水平區(qū)間的基因數(shù)量統(tǒng)計(jì)表Table 4 Gene number of each sample in different interval of expression level

圖2 堿脅迫處理下不同時(shí)間點(diǎn)的DEGs 數(shù)量與0 h 的對(duì)比圖Fig.2 Number of DEGs at different time points after under alkaline stress compared with 0 hour

此外,在6 h vs.0 h 中,28 個(gè)基因(c57635_g1_i1、c64583_g1_i1、c64900_g1_i1、c71351_g1_i1、c72539_g1_i1、c73498_g1_i4、c73839_g1_i2、c73869_g1_i2、c74641_g1_i1、c75112_g1_i1、c76093_g2_i1、c76646_g1_i1、c78345_g1_i1、c78766_g1_i1、c79591_g1_i1、c79591_g1_i2、c79591_g1_i3、c81104_g2_i1、c82384_g1_i1、c82500_g1_i1、c83082_g2_i1、c83291_g1_i3、c83618_g1_i4、c83799_g2_i1、c86513_g2_i1、c87340_g1_i1、c87510_g1_i2 和c87918_g1_i2)比對(duì)到對(duì)刺激的應(yīng)答(response to stimulus)GO 條目中,表明這些基因可能在堿脅迫響應(yīng)機(jī)制中起著重要的作用。

2.3.2 KEGG pathway 富集分析 由圖6可知,3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的基因KEGG 分析結(jié)果存在著較大的差異。3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 中分別共有22、403 和1 006 個(gè)DEGs 被注釋到42、41 和47 個(gè)代謝途徑中。對(duì)于3 h vs.0 h,主要涉及的代謝途徑分支有碳代謝途徑[3,(13.64%)]、原核生物碳固定途徑[3,(13.64%)]和萜類骨架生物合成途徑[3,(13.64%)]。對(duì)于6 h vs.0 h,前三大富集途徑依次為碳代謝途徑[37,(9.18%)]、淀粉和蔗糖代謝途徑[23,(5.71%)]和植物-病原體相互作用途徑[21,(5.21%)]。在12 h vs.0 h 分配的47 個(gè)代謝途徑分支中,碳代謝途徑[77,(7.65%)]也是最富集的途徑,其次是氨基酸的生物合成途徑[65,(6.46%)]、淀粉和蔗糖代謝途徑[53,(5.27%)]及植物激素信號(hào)轉(zhuǎn)導(dǎo)途徑[46,(4.57%)]。

圖3 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 細(xì)胞組分功能類型對(duì)比圖Fig.3 Comparative distribution of GO enrichment of terms related to cellular components for 3 h vs.0 h， 6 h vs.0 h and 12 h vs.0 h

3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的基因KEGG 分析結(jié)果表明,榆葉梅在這些時(shí)間下均持續(xù)涉及了以下幾個(gè)重要的代謝途徑：碳代謝途徑、光合生物碳固定途徑、苯丙素生物合成、絲氨酸和蘇氨酸代謝、甲烷代謝途徑、丙酮酸代謝途徑、乙醛酸和二羧酸代謝途徑、萜類骨架生物合成途徑、α-亞麻酸代謝、色氨酸代謝途徑、檸檬烯和蒎烯降解途徑。在堿脅迫3、6、12 h 條件下,均涉及“碳代謝途徑”的相關(guān)基因?yàn)閏86254_g5_i1、c86254_g2_i1、c58505_g1_i1。對(duì)這些碳代謝相關(guān)基因進(jìn)行KEGG 注釋表明,c86254_g5_i1 和c86254_ g2 _ i1 均為磷酸烯醇式丙酮酸羧化酶(phosphoenolpyruvate carboxylase,PEPC,EC 4.1.1.31)基因,而c58505_g1_i1 為乙酰CoA 酰基轉(zhuǎn)移酶(acetyl-CoA C-acetyltransferase,AACT)基因。

對(duì)28 個(gè)比對(duì)到對(duì)刺激應(yīng)答GO 條目中的基因進(jìn)行KEGG 注釋表明,c71351_g1_i1 是一個(gè)RNA 聚合酶相關(guān)蛋白,而c73869_g1_i2 是一個(gè)HSP20 家族蛋白。對(duì)這些基因進(jìn)行COG 功能注釋,發(fā)現(xiàn)有5 個(gè)基因比對(duì)到“翻譯后修飾,蛋白質(zhì)折疊和分子伴侶類基因”(c73869_g1_i2,c74641_g1_i1,c78766_g1_i1),“復(fù)制、重組和修復(fù)類基因(replication,recombination and repair)”(c82384_g1_i1),和“核苷酸轉(zhuǎn)運(yùn)和代謝類基因(nucleotide transport and metabolism)”(c83082_g2_i1)中。

2.4 RT-qPCR 驗(yàn)證

采用基因特異性引物的RT-qPCR 分析對(duì)RNASeq 檢測(cè)到的與堿脅迫相關(guān)的DEGs 進(jìn)行了驗(yàn)證。RTqPCR 驗(yàn)證選取了7 個(gè)候選基因(c78498_g1_i1、c80260_g1_i1、c81662_g1_i1、c84090_g1_i1、c85502_g2_i1、c85846_g1_i3 和c88652_g1_i2),并將ACTIN 基因作為數(shù)據(jù)歸一化的內(nèi)參基因。RNA-Seq 數(shù)據(jù)分析顯示,這7 個(gè)基因均屬于上調(diào)基因,且基因的大多數(shù)表達(dá)模式與通過RT-qPCR 驗(yàn)證得到的表達(dá)模式基本一致,且RT-qPCR 驗(yàn)證中這些基因在3 種堿脅迫時(shí)間處理下均呈極顯著性差異表達(dá)(P＜0.01)(圖7)。

圖4 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 分子功能類型對(duì)比圖Fig.4 Comparative distribution of GO enrichment of terms related to molecular function for 3 h vs.0 h， 6 h vs.0 h and 12 h vs.0 h

3 討論

本研究運(yùn)用Illumina HiSeqTM2500 技術(shù)對(duì)無(wú)參考基因組的榆葉梅堿脅迫處理不同時(shí)間下的葉片進(jìn)行了轉(zhuǎn)錄組測(cè)序,共獲得了將近5 300 萬(wàn)條高質(zhì)量clean reads,并經(jīng)過de novo組裝后獲得了124 786 條榆葉梅unigenes,獲得的大量榆葉梅轉(zhuǎn)錄本序列將為榆葉梅物種提供重要的基因信息,可用于堿脅迫下重要功能基因的挖掘。本研究對(duì)堿脅迫下榆葉梅的DEGs 進(jìn)行研究,發(fā)現(xiàn)榆葉梅葉片在12 h vs.0 h 比較中篩選出最多的顯著性DEGs,這與星星草(Puccinellia tenuiflor)[29]和野生大豆(Glycine soja)[30]根中所獲得的結(jié)果不一致,其星星草(Puccinellia tenuiflor)和野生大豆(Glycine soja)在6 h 篩選出的DEGs 量比在12 h 和24 h 更高,這可能是因?yàn)槿~細(xì)胞響應(yīng)外界刺激的速度慢于根細(xì)胞。本研究在堿脅迫3 h 下未檢測(cè)到顯著性DEGs,這與野生大豆(G.soja)[17]根中所得結(jié)果類似。此外,本研究3 h vs.0 h、6 h vs.0 h、12 h vs.0 h 中呈上調(diào)的DEGs 數(shù)量比下調(diào)的多,這與在多種非生物脅迫下的擬南芥轉(zhuǎn)錄組結(jié)果基本一致[30]。

DEGs 的功能分類將闡明其在細(xì)胞內(nèi)的分子功能和生物代謝途徑,對(duì)挖掘榆葉梅抗堿基因具有重要的意義。本研究中,GO 功能顯著性富集分析表明,3 h vs.0 h 中只有少量顯著性DEGs 注釋到GO 條目中,這些GO 條目?jī)H包括分子功能類型中的“核酸結(jié)合轉(zhuǎn)錄因子活性”和“序列特異性DNA 結(jié)合轉(zhuǎn)錄因子活性”,以及生物學(xué)過程類型中的“單-生物代謝過程”、“氧脂代謝過程”和“氧脂生物合成過程”。這些GO 條目反映了榆葉梅對(duì)堿脅迫的早期響應(yīng),以及這些條目可能誘導(dǎo)產(chǎn)生后期的級(jí)聯(lián)反應(yīng)。因此,這些基因在早期抗堿機(jī)制中可能起著重要的作用。本試驗(yàn)中,KEGG pathway 富集分析結(jié)果表明,榆葉梅葉片中代謝途徑隨堿脅迫處理時(shí)間的延長(zhǎng)而有所變化。但在堿脅迫處理3、6 和12 h 條件下,榆葉梅均持續(xù)涉及了幾個(gè)重要的代謝途徑。其中“碳代謝途徑”涉及最多的DEGs,表明榆葉梅為了應(yīng)對(duì)堿脅迫需要大量的能量供應(yīng),以及堿脅迫可能主要通過改變碳代謝相關(guān)基因的表達(dá)來(lái)影響基因轉(zhuǎn)錄情況。本試驗(yàn)發(fā)現(xiàn)PEPC和AACT兩個(gè)基因在逆境脅迫中具有重要作用。PEPC 主要在植物細(xì)胞中參與植物的光合碳同化等重要代謝途徑,以及調(diào)節(jié)細(xì)胞pH 值,保持離子平衡,調(diào)節(jié)氣孔保衛(wèi)細(xì)胞的滲透壓及運(yùn)動(dòng)等[31]。AACT 為萜類化合物生物合成甲羥戊酸(mevalonic acid,MVA)途徑的起始酶,催化使2個(gè)分子的乙酰CoA 縮合為乙酰乙酰CoA,主要參與植物的呼吸作用等[32]。此外,一些代謝途徑僅在某一處理時(shí)期涉及有DEGs。在堿脅迫處理3 h 時(shí),發(fā)現(xiàn)一個(gè)特異基因(c58505_g1_i1)涉及了“鈣信號(hào)途徑(calcium signaling pathway)”,另一個(gè)特異基因(c89089_g2_i1)涉及了“ABC 轉(zhuǎn)運(yùn)蛋白途徑(ABC transporters pathway)”。鈣信號(hào)途徑是生物體內(nèi)生理生化反應(yīng)的重要控制途徑,研究證明其能間接影響植物的生理活動(dòng),如氣孔關(guān)閉、花粉管發(fā)育,根系伸長(zhǎng)等[33]。鈣信號(hào)途徑在堿脅迫處理3 h 時(shí)的特異響應(yīng),有可能與幫助榆葉梅應(yīng)對(duì)堿脅迫,形成胞內(nèi)鈣離子濃度振蕩變化,進(jìn)而引發(fā)氣孔關(guān)閉,抵御離子毒害有關(guān)。ABC(ATP-Binding Cassette)轉(zhuǎn)運(yùn)蛋白是目前已知最大、功能最廣泛的蛋白家族,其可能涉及植物次生代謝物的跨膜轉(zhuǎn)運(yùn),從而保護(hù)植物免受環(huán)境脅迫的傷害[34-36]。這些參與早期特異性途徑的基因在榆葉梅早期堿性響應(yīng)機(jī)制中可能起著重要的作用,并可能在后續(xù)植物響應(yīng)堿脅迫中誘導(dǎo)產(chǎn)生級(jí)聯(lián)反應(yīng)。本研究中,在堿脅迫處理6、12 h 條件下,開始涉及許多重要代謝途徑,包括淀粉和蔗糖代謝途徑、植物激素信號(hào)轉(zhuǎn)導(dǎo)途徑、糖酵解/糖異生途徑、光合作用途徑、類胡蘿卜素生物合成途徑,以及光合作用-天線蛋白途徑等,表明這些參與代謝途徑的基因隨著堿脅迫時(shí)間的增加而呈多樣性表達(dá)模式,并且榆葉梅在堿脅迫適應(yīng)過程中存在快速、多樣的信號(hào)感知、轉(zhuǎn)導(dǎo)以及相應(yīng)復(fù)雜的分子機(jī)制。

圖5 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 生物學(xué)過程類型對(duì)比圖Fig.5 Comparative distribution of GO enrichment of terms related to biological processes for 3 h vs.0 h， 6 h vs.0 h and 12 h vs.0 h

本試驗(yàn)結(jié)果表明,RT-qPCR 所得基因的表達(dá)模式與RNA-Seq 數(shù)據(jù)結(jié)果類似,但它們間仍存在一些差異,這可能是RNA-Seq 和RT-qPCR 兩種方法使用的算法和原理不同而致[37-38],這與Yates 等[39]和Huang等[40]的研究結(jié)果類似。此外,這7 個(gè)基因在3 種堿脅迫時(shí)間處理下均呈極顯著性差異表達(dá)(P＜0.01),表明它們可能在響應(yīng)堿脅迫過程中發(fā)揮作用?？傊?本研究從RT-qPCR 和RNA-Seq 方法中獲得了類似的基因表達(dá)模式,證實(shí)了測(cè)序結(jié)果的可靠性和準(zhǔn)確性。對(duì)于重要的差異表達(dá)基因,后續(xù)還將進(jìn)一步用northern blot對(duì)其進(jìn)行驗(yàn)證。

圖6 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h DEGs 的KEGG 富集分析Fig.6 KEGG enrichment analysis of the DEGs of 3 h vs.0 h， 6 h vs.0 h and 12 h vs.0 h

葉片作為植物進(jìn)行光合作用和生命活動(dòng)的主要器官,是對(duì)逆境刺激反應(yīng)最為敏感的部位。本研究主要通過葉片的分子變化來(lái)探索榆葉梅應(yīng)答堿脅迫的機(jī)制。目前,在薔薇科果樹上有研究報(bào)道[41],水分脅迫下,新根內(nèi)皮層細(xì)胞出現(xiàn)細(xì)胞壁增厚現(xiàn)象,對(duì)于榆葉梅根系在堿脅迫下的響應(yīng)機(jī)制將進(jìn)行后續(xù)研究。本研究發(fā)現(xiàn)了一些堿脅迫下差異表達(dá)的基因,可能對(duì)應(yīng)調(diào)控不同的物質(zhì)或途徑,為將來(lái)采用分子生物學(xué)的手段將單個(gè)抗性基因?qū)胫仓牦w,獲得抗堿轉(zhuǎn)基因材料、培育耐堿果樹新品種創(chuàng)造前提條件。此外,本研究是針對(duì)榆葉梅砧木本身進(jìn)行堿脅迫響應(yīng)機(jī)制的分析,而嫁接苗的抗性機(jī)制與砧木又不完全相同,后續(xù)會(huì)以榆葉梅為砧木嫁接不同李品種,并與毛桃砧木嫁接苗對(duì)比研究,進(jìn)一步探討砧木對(duì)李抗堿性的影響及其內(nèi)在機(jī)理。

圖7 榆葉梅7 個(gè)差異表達(dá)耐堿基因的RNA-Seq 和RT-qPCR 表達(dá)模式圖Fig.7 Expression patterns of seven genes selected from alkalinity-responsive DEGs of Prunus triloba by RNA-Seq and RT-qPCR

4 結(jié)論

本研究對(duì)榆葉梅在單一堿脅迫下的分子響應(yīng)機(jī)制進(jìn)行了探索。使用Illumina 測(cè)序平臺(tái)對(duì)榆葉梅在短期堿脅迫下的葉片進(jìn)行了de novo轉(zhuǎn)錄組分析,獲得了大量的DEGs,可用于榆葉梅今后的分子育種中?；谶@些DEGs,進(jìn)行GO 富集分析發(fā)現(xiàn),其中28 個(gè)基因可能在早期堿脅迫響應(yīng)機(jī)制中起著重要的作用；而對(duì)這些DEGs 進(jìn)行KEGG 富集分析發(fā)現(xiàn),不同堿處理時(shí)間下的KEGG 途徑具有顯著差異。此外,采用RT-qPCR 法驗(yàn)證了7 個(gè)堿性相關(guān)基因的表達(dá)模式,證實(shí)了RNASeq 結(jié)果的可靠性。本研究獲得了關(guān)于榆葉梅對(duì)堿脅迫的短期分子響應(yīng)機(jī)制的相關(guān)數(shù)據(jù),可為了解榆葉梅及其他相關(guān)果樹的短期堿脅迫分子適應(yīng)機(jī)制奠定基礎(chǔ)。