劉 佳 田 云 鄧家林

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

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

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

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

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

1 材料與方法

1.1 試驗材料及堿脅迫處理

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

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

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

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

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

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

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

1.4 基因表達定量與差異表達分析

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

1.5 RT-qPCR 驗證

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

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

2 結(jié)果與分析

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

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

本研究將獲得的全部unigenes 與Nr 數(shù)據(jù)庫進行了比對注釋。結(jié)果顯示,有52 691條 unigenes(42.23%)被成功注釋功能。其中對注釋到的近源物種進行統(tǒng)計,發(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)計結(jié)果Table 2 Summary of sequences analysis

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

2.2 基因表達定量分析與DEGs 的識別

采用FPKM 法估計了每個樣品的基因表達水平。在這10 個樣本中,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 個顯著性DEGs 分別于1 h vs.0 h、3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 比較中識別(圖2),12 h vs.0 h 比較中篩選出最多的顯著性DEGs。在8 948個顯著性DEGs 中,3 h vs.0 h、6 h vs.0 h、12 h vs.0 h中分別有181、1 999、4 895 個DEGs 上調(diào)；3 h vs.0 h、6 h vs.0 h、12 h vs.0 h 中分別有5、1 212 和3 452 個DEGs 下調(diào)(圖2)。

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

2.3 DEGs 功能分類

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

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

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

此外,在6 h vs.0 h 中,28 個基因(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)比對到對刺激的應答(response to stimulus)GO 條目中,表明這些基因可能在堿脅迫響應機制中起著重要的作用。

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 個DEGs 被注釋到42、41 和47 個代謝途徑中。對于3 h vs.0 h,主要涉及的代謝途徑分支有碳代謝途徑[3,(13.64%)]、原核生物碳固定途徑[3,(13.64%)]和萜類骨架生物合成途徑[3,(13.64%)]。對于6 h vs.0 h,前三大富集途徑依次為碳代謝途徑[37,(9.18%)]、淀粉和蔗糖代謝途徑[23,(5.71%)]和植物-病原體相互作用途徑[21,(5.21%)]。在12 h vs.0 h 分配的47 個代謝途徑分支中,碳代謝途徑[77,(7.65%)]也是最富集的途徑,其次是氨基酸的生物合成途徑[65,(6.46%)]、淀粉和蔗糖代謝途徑[53,(5.27%)]及植物激素信號轉(zhuǎn)導途徑[46,(4.57%)]。

圖3 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 細胞組分功能類型對比圖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é)果表明,榆葉梅在這些時間下均持續(xù)涉及了以下幾個重要的代謝途徑：碳代謝途徑、光合生物碳固定途徑、苯丙素生物合成、絲氨酸和蘇氨酸代謝、甲烷代謝途徑、丙酮酸代謝途徑、乙醛酸和二羧酸代謝途徑、萜類骨架生物合成途徑、α-亞麻酸代謝、色氨酸代謝途徑、檸檬烯和蒎烯降解途徑。在堿脅迫3、6、12 h 條件下,均涉及“碳代謝途徑”的相關(guān)基因為c86254_g5_i1、c86254_g2_i1、c58505_g1_i1。對這些碳代謝相關(guān)基因進行KEGG 注釋表明,c86254_g5_i1 和c86254_ g2 _ i1 均為磷酸烯醇式丙酮酸羧化酶(phosphoenolpyruvate carboxylase,PEPC,EC 4.1.1.31)基因,而c58505_g1_i1 為乙酰CoA ?；D(zhuǎn)移酶(acetyl-CoA C-acetyltransferase,AACT)基因。

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

2.4 RT-qPCR 驗證

采用基因特異性引物的RT-qPCR 分析對RNASeq 檢測到的與堿脅迫相關(guān)的DEGs 進行了驗證。RTqPCR 驗證選取了7 個候選基因(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 個基因均屬于上調(diào)基因,且基因的大多數(shù)表達模式與通過RT-qPCR 驗證得到的表達模式基本一致,且RT-qPCR 驗證中這些基因在3 種堿脅迫時間處理下均呈極顯著性差異表達(P＜0.01)(圖7)。

圖4 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 分子功能類型對比圖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 討論

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

圖5 3 h vs.0 h、6 h vs.0 h 和12 h vs.0 h 的GO 生物學過程類型對比圖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

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

圖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

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

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

4 結(jié)論

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