基于RNA-Seq挖掘桃低溫應答相關基因的研究

李小蘭 郝蘭蘭 張 帆 王 鴻,*

(1 甘肅省農(nóng)業(yè)科學院林果花卉研究所，甘肅 蘭州 730070；2 甘肅農(nóng)業(yè)大學園藝學院，甘肅 蘭州 730070)

桃(PrunuspersicaL.)原產(chǎn)于我國[1]，屬薔薇科(Rosaceae)李屬(Prunus)桃亞屬(Amygdalus)，是我國的主栽果樹之一[2]。溫度是影響桃樹生長發(fā)育以及地理分布的一個重要環(huán)境因子。我國北方地區(qū)冬季氣候嚴寒,如果發(fā)生凍害，會給果樹生產(chǎn)造成嚴重的損失。此外，“倒春寒”的發(fā)生，也會使花芽受到損傷，嚴重影響果樹產(chǎn)量。近年來，轉錄組測序RNA-Seq為植物響應低溫脅迫的研究提供了一個新的視角，可以高效、快捷地反映生物在某一特定壞境和時間下基因的表達情況[3]，該技術已成為解析植物應答低溫脅迫的分子機制以及挖掘功能基因的重要手段[4-5]。

植物不同的器官組織在低溫誘導下表現(xiàn)出不同的表型和抗寒反應[6]。張麗[7]通過對低溫儲藏過程中兩種桃果實進行轉錄組分析，發(fā)現(xiàn)ID為Prupe.1G220700的基因可能在桃的低溫儲藏過程中發(fā)揮重要作用。前人對已經(jīng)開花的桃樹枝條進行不同時間的低溫處理，隨后分離雌蕊進行高通量測序，發(fā)現(xiàn)與碳水化合物水解和低溫響應相關的基因表達量顯著上調[8]。Niu等[9]和Renaut等[10]分別對一年生枝條和樹干皮層組織進行不同程度的低溫誘導，并對獲得的差異基因進行生物信息學分析，發(fā)現(xiàn)與碳水化合物代謝途徑以及信號轉導途徑相關的基因在桃樹抗寒過程中發(fā)揮重要作用。前人已經(jīng)從不同層面揭示桃樹的花[8]、果實[11-12]、枝條[9]和樹干皮層組織[10]在低溫誘導下的基因表達水平變化，但針對低溫脅迫下桃葉片應激反應的研究較少。葉片作為重要的營養(yǎng)器官，在植物生長發(fā)育過程中發(fā)揮重要作用。因此，研究低溫脅迫下桃葉片的抗寒反應，明確響應低溫脅迫的分子調控機制，可為抗寒桃品種的培育提供理論基礎。

甘肅河西走廊地區(qū)平均最低氣溫約為-27℃[13]，引進的桃品種無法適應這種環(huán)境。然而，屬于李光桃抗寒品種群的丁家壩油桃品種能很好地適應這種環(huán)境，進化出優(yōu)良的耐寒特性。本研究選取低溫耐受型丁家壩[14]和敏感型加納巖[14]桃葉片作為試驗材料，比較其在低溫脅迫下的生理生化水平變化，并利用生物信息學方法對丁家壩和加納巖桃葉片在正常和低溫條件下的轉錄組進行比較分析，鑒定其響應低溫的相關基因及代謝途徑，并對相關代謝途徑進行深入研究，以期闡明桃樹葉片抗寒的分子機制，為研究桃抗寒的分子機制提供新線索，也為桃抗寒育種提供新的候選基因。

1 材料與方法

1.1 試驗材料

試驗材料取自甘肅省農(nóng)業(yè)科學院桃園中長勢優(yōu)良的冷耐受型丁家壩(D)和冷敏感型加納巖(K)七年生桃樹葉片。

1.2 試驗處理

2021年7月下旬，于甘肅省農(nóng)業(yè)科學院桃園各采集12枝粗細相近、長短相同、朝北直立生長的丁家壩和加納巖一年生帶葉枝條。將每個品種的枝條隨機分為4組，用保鮮袋密封，放在4℃培養(yǎng)箱中進行低溫處理。在處理0(D0/K0)、2(D2/K2)、6(D6/K6)和12 h(D12/K12)時，每個品種隨機選取一組枝條采集葉片，每組3個枝條分別作為3個重復。選取從上往下數(shù)第6到第9片葉作為試驗材料，立即用液氮冷凍(相對電導率用鮮葉測定)，在-80℃下保存?zhèn)溆谩?/p>

1.3 葉片生理生化指標測定

采用磺基水楊酸—酸性茚三酮法測定游離脯氨酸(Proline, Pro)含量[15]；采用蒽酮比色法測定可溶性糖含量[16]；考馬斯亮藍G-250染色法測定可溶性蛋白含量[17]；采用硫代巴比妥酸(thiobarbituric acid, TBA)反應測定丙二醛(malondiadehyde, MDA)含量[18]，并使用Yun等[19]描述的方法計算電解質滲透指數(shù)(electrolyte leakage index, ELI)，使用POD試劑盒(北京索萊寶科技有限公司)測定過氧化物酶含量(peroxidase, POD)。數(shù)據(jù)采用Excel 2010和SPSS 20.0進行整理分析。

1.4 RNA提取、文庫構建及測序

用植物總RNA提取試劑盒(中國天根生物技術有限公司，北京)提取總RNA。轉錄組測序由南京派森諾基因科技有限公司完成。采用瓊脂糖凝膠電泳法檢測RNA的完整性及是否存在DNA污染。使用納米分光光度計(IMPLEN，德國)檢測RNA純度。建庫中使用的建庫試劑盒為NEBNext? UltraTM RNA Library Prep Kit(Illumina，美國)。庫檢合格后，對不同文庫進行Illumina HiSeq X Ten PE150高通量測序。

1.5 數(shù)據(jù)處理及拼裝

對獲得的原始讀數(shù)(raw reads)進行過濾，去除帶接頭(adapter)的reads、含N的reads以及低質量的reads。同時，對干凈讀數(shù)(clean reads)進行Q20和Q30含量計算。使用HISAT2 v2.1.0[20]將配對末端clean reads與參考基因組比對。使用featureCounts(1.5.0-p3)計算映射到每個基因的讀數(shù)，然后根據(jù)基因的長度計算每個基因的FPKM(每千個堿基轉錄每百萬映射讀取的fragments，F(xiàn)ragments Per Kilobase of exon model per Million mapped fragments)。

1.6 差異基因富集分析

采用DESeq[21]對基因表達進行差異分析，|log2FoldChange|>1、P<0.05作為差異基因篩選條件。使用topGO[22]進行GO富集分析，找出差異基因顯著富集的GO term。使用clusterProfiler[23]軟件進行京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)通路富集分析，探究P<0.05的顯著富集通路。

1.7 實時熒光定量PCR(quantitative real-time PCR，qRT-PCR)表達量驗證

選取15個差異基因通過qRT-PCR對測序結果進行驗證，使用與測序同批次的RNA樣本，按照Prime Script TM RT reagent Kit反轉錄試劑盒(TaKaRa, 大連)說明書合成cDNA。使用Primer 3.0在線設計引物，由上海生物工程股份有限公司合成，Actin作為內參基因(表1)。使用Light Cycler?96 Real-Time PCR System PCR儀(Roche，瑞士)進行擴增，擴增體系為20 μL：上下游引物各1 μL，cDNA 2 μL，ddH2O 6 μL，TB Green-II 10 μL。反應程序：95℃變性5 min；95℃變性10 s，60℃退火10 s，72℃延伸10 s，40個循環(huán)。本試驗進行3次生物學重復，目的基因表達量的統(tǒng)計方法采用相對定量法(2-△△CT)。

表1 本研究所需的引物信息

2 結果與分析

2.1 桃葉片在冷處理下的生理生化變化

冷處理下，POD可以阻止植物細胞氧化損傷，減輕脅迫危害。由圖1-A可知，D和K桃葉片在冷處理12 h的POD活性與對照相比分別上升了58.94%和33.66%。

MDA含量和ELI與植物細胞膜受損傷程度及植物抗逆性強弱相關。D和K葉片的MDA含量(圖1-B)和ELI(圖1-C)均隨著冷處理時間的延長呈上升趨勢，且D的增長幅度小于K，表明D的細胞膜不易被破壞，抗寒性強。

注：不同小寫字母表示同一品種在不同時間低溫處理下差異顯著(P<0.05)；*和**分別代表不同品種在相同時間低溫處理下差異顯著(P<0.05)和極顯著(P<0.01)。

脯氨酸(Pro)、可溶性糖和可溶性蛋白含量的積累可以增強植物對逆境脅迫的耐受性，是植物在逆境下采取的自我保護措施。總體來看，D和K的脯氨酸、可溶性糖以及可溶性蛋白含量在冷處理12 h時與對照相比顯著上升，且D的增長幅度大于K。

2.2 轉錄組測序數(shù)據(jù)統(tǒng)計與質量評估

分別對D和K兩種桃葉片經(jīng)4℃處理0、2、6和12 h并進行取樣，每個樣本進行3次生物學重復，共建立24個cDNA文庫(表2)。經(jīng)過高通量測序后，樣本產(chǎn)生的原始序列(raw reads)數(shù)為4 300萬～5 700萬，過濾后的序列(clean reads)為4 000萬～5 300萬。比對到參考基因組(GCF_000 346 465.2)的序列總數(shù)占clean reads的87.59%～97.00%。對每個樣本的clean reads進行統(tǒng)計后發(fā)現(xiàn)Q20平均達到97.61%，Q30值平均達到 93.27%。以上結果說明測序所得的cDNA文庫質量較高，可以進一步進行生物信息學分析。

表2 測序數(shù)據(jù)的統(tǒng)計與質量評估

2.3 低溫脅迫下差異表達基因的分析

為了篩選與低溫響應相關的基因，分別對兩種桃葉片不同時間段低溫處理產(chǎn)生的差異基因進行比較分析(圖2)。D在低溫脅迫2、6和12 h時與對照相比分別有659、1 632和2 458個基因上調，534、1 064和1 951個基因下調(圖2-A)，3個處理時間段共同表達的DEGs有569個(圖2-B)。與對照相比，K在低溫脅迫2、6和12 h時分別有407、1 742和2 016個基因上調，506、946和1 815個基因下調(圖2-A)，共同表達的差異基因有505個(圖2-B)。

圖2 兩種桃葉片在不同時間低溫處理下的差異基因柱形圖(A)和韋恩圖(B)

2.4 差異基因GO富集分析

在D對比組中，569個共同表達的差異基因富集在1 097個GO項(基因功能國際標準分類體系，Gene Ontology)，主要包括分子功能(288個)、生物過程(714個)和細胞組分(95個)。根據(jù)錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)值選擇富集最顯著的20個GO項作圖(圖3-A)。顯著富集中的前5個GO項為：屬于分子功能的轉錄調節(jié)活性(GO:0140110)、DNA結合轉錄因子活性(GO:0003700)和分子功能(GO:0003674)以及屬于細胞組分的膜的固有成分(GO:0031224)和屬于生物過程的生物學過程(GO:0008150)。

注: 1: 轉錄調節(jié)活性; 2: 轉錄，DNA模板; 3: 序列特異性DNA結合; 4: 轉錄調控，DNA模板; 5: RNA代謝過程的調控; 6: RNA生物合成過程的調控; 7: 初級代謝過程的調節(jié); 8: 含堿基的化合物代謝過程的調控; 9: 核酸模板轉錄的調控; 10: 分子功能; 11: 膜; 12: 膜的固有成分; 13: 膜的組成部分; 14: 雜環(huán)化合物結合; 15: DNA結合轉錄因子活性; 16: DNA結合; 17: 細胞組成; 18: 細胞解剖實體; 19: 生物過程; 20: 結合; 21: 轉移己糖基的轉移酶活性; 22: 東莨菪堿β-葡萄糖苷酶活性; 23: 細胞過程的調節(jié); 24: 細胞生物合 成過程的調節(jié); 25: 生物合成過程的調控; 26: 葡萄糖苷酶活性; 27: β-葡萄糖苷酶活性。

在K對比組中，505個共同表達的差異基因富集到1 021個GO項，其中分子功能(268個)、細胞組分(79個)以及生物過程(674個)。挑選FDR值最小即富集最顯著的20個GO項進行展示(圖3-B)。前5個顯著富集的GO項為：MF的DNA結合轉錄因子活性(GO:0003700)、轉錄調節(jié)活性(GO:0140110)、序列特異性DNA結合蛋白(GO:0043565)和β-葡萄糖苷酶活性(GO:0008422)以及BP的生物學過程(GO:0008150)。

2.5 差異基因KEGG富集分析

為了進一步分析DEGs的生物學功能，采用KEGG數(shù)據(jù)庫對DEGs進行通路富集分析。在D對比組中，569個差異基因涉及57個KEGG途徑，包括細胞過程、遺傳信息處理、環(huán)境信息處理、新陳代謝、有機體系統(tǒng)等五個方面，大部分途徑屬于新陳代謝類別。顯著富集的代謝通路有植物與病原菌互作(plant-pathogen interaction)、MAPK信號通路-植物(MAPK signaling pathway-plant)和倍半萜和三萜生物合成(Sesquiterpenoid and triterpenoid biosynthesis)(圖4-A)，分別涉及11、8和3個相關基因。

K中的505個DEGs被分配到47個KEGG途徑中，包括細胞過程、遺傳信息處理、環(huán)境信息處理、新陳代謝、有機體系統(tǒng)等五個方面。最顯著富集的途徑是氰基氨基酸代謝(cyanoamino acid metabolism)和內質網(wǎng)蛋白加工(protein processing in endoplasmic reticulum)(圖4-B)。以上結果表明，低溫脅迫對兩種抗寒性不同的基因型的影響不同。比較D和K富集到的KEGG通路，重點研究D中顯著富集的代謝通路，如植物與病原菌互作(plant-pathogen interaction)和MAPK信號通路-植物(MAPK signaling pathway-plant)。

注: 1: 胞吞作用; 2: MAPK信號通路-植物; 3: 植物激素信號轉導; 4: 磷脂酰肌醇信號系統(tǒng); 5: 內質網(wǎng)中的蛋白質加工; 6: 倍半萜和三萜生物合成; 7: 萜類化合物生物合成支柱; 8: 甜菜色素生物合成; 9: α-亞麻酸代謝; 10: 油菜素類固醇生物合成; 11: 酮體的合成與降解; 12: 二萜生物合成; 13: 類胡蘿卜素生物合成; 14: 苯丙氨酸、酪氨酸和色氨酸生物合成; 15: 卟啉與葉綠素代謝; 16: 核黃素代謝; 17: 半乳糖代謝; 18: 硫胺素代謝; 19: 亞油酸; 20: 晝夜節(jié)律-植物; 21: 光合作用-天線蛋白; 22: 生物素代謝; 23: 類黃酮生物合成; 24: 丁酸代謝; 25: 葉酸生物合成; 26: 莨菪堿、哌啶和吡啶生物堿的生物合成; 27: 淀粉和蔗糖代謝; 28: 角質、木栓質和蠟的生物合成; 29: 植物與病原菌互作; 30: 苯丙烷的生物合成; 31: 蛋白質輸出; 32: 氰基氨基酸代謝; 33: 芥子油苷生物合成; 34: N-聚糖生物合成; 35: 精氨酸和脯氨酸代謝; 36: 賴氨酸生物合成; 37: 不飽和脂肪酸生物合成; 38: 煙酸鹽和煙酰胺代謝; 39: 色氨酸代謝; 40: 苯丙烷代 謝; 41: 脂肪酸生物合成; 42: 脂肪酸降解。

2.6 丁家壩抗寒相關基因

對丁家壩顯著富集的代謝通路進行分析表明，在植物與病原菌互作途徑，由真菌幾丁質、Ca2+、細菌鞭毛等不同PAMPs激發(fā)的反應中，CaM、CDPK、CML和WRKY等相關基因轉錄本顯著上調，促進ROS和一氧化氮(NO)的產(chǎn)生，參與丁家壩抗寒反應。此外，MAPK信號通路-植物代謝通路中FLS2受體識別來自細菌鞭毛蛋白的flg22肽，激發(fā)植物的MAPK級聯(lián)反應，引起細胞內轉錄及代謝水平變化，進一步響應低溫脅迫。對植物與病原菌互作和MAPK信號通路-植物途徑涉及的基因進行熱圖分析(圖5)，結果表明CPK26、CALM7、CML45、CML31、CML24、CML27、CML38、WRKY22、WRKY24和WRKY33等相關基因在丁家壩抗寒反應中發(fā)揮重要作用。

圖5 丁家壩抗寒相關基因熱圖

2.7 qRT-PCR驗證

為驗證轉錄組測序數(shù)據(jù)的可靠性，從轉錄組鑒定到的差異基因中選取了15個可能參與低溫響應的基因進行實時熒光定量分析。如圖6所示，qRT-PCR分析的基因表達模式與RNA-Seq方法分析的基因表達模式相似，證實了RNA-Seq數(shù)據(jù)的可靠性。

注：qRT-PCR：直方圖；FPKM：折線圖。

3 討論

低溫是植物生長發(fā)育過程中常見的環(huán)境脅迫因子，植物通過改變其形態(tài)和生理生化特性來適應低溫脅迫[24]。前人研究表明，植物對低溫脅迫的響應與許多基因的表達有關[25]。轉錄組測序技術的快速發(fā)展為解析植物響應低溫脅迫的分子機制和挖掘抗寒關鍵基因提供了新的途徑[26-27]。本研究對冷處理下的耐冷型丁家壩和敏感型加納巖桃葉片的生理生化水平變化進行研究，并利用轉錄組技術對其冷處理下的轉錄組進行對比和分析，以鑒定出與丁家壩抗寒性相關的調控機制及關鍵基因。結果表明，丁家壩和加納巖桃葉片的電解質外滲率和MDA含量整體呈上升趨勢，在冷處理6 h時出現(xiàn)激增。與冷處理2 h相比，6 h丁家壩和加納巖的ELI分別增加了45%和53%，MDA含量分別增加21%和24%。這與Megha等[28]的研究中關于橄欖型油菜ELI隨著冷處理時間的延長逐漸增加的結果相同；加納巖在冷處理各時期的MDA含量均高于丁家壩，這與低溫脅迫下甜瓜不同品種間的耐低溫性與MDA含量成顯著負相關的研究結果一致[29]。MDA含量的積累是ROS過度積累導致的，抗氧化物酶在植物活性氧解毒過程中發(fā)揮重要作用[30]。本研究中，冷處理12時，MDA增長率在與6 h相比有所下降，推測可能與POD活性在冷處理12 h時與6 h相比顯著上升有關。前人研究表明，植物體內可溶性糖、可溶性蛋白和脯氨酸等重要滲透調節(jié)物質的含量與植物自身的抗寒性成顯著正相關[31-32]。本研究中，丁家壩和加納巖的可溶性蛋白、可溶性糖以及脯氨酸的含量隨著冷處理時間的延長逐漸升高，且丁家壩的增長幅度高于加納巖，與前人研究結果一致。

轉錄組分析結果表明，丁家壩在冷處理下的差異基因顯著富集在植物與病原菌互作、MAPK信號通路-植物和倍半萜和三萜生物合成等代謝通路。信號轉導途徑是冷感機制和遺傳反應之間的聯(lián)系[33]。前人已經(jīng)研究并闡明了可能與植物冷響應相關的信號通路[34]。Niu等[9]對丁家壩休眠枝在不同程度低溫脅迫下的轉錄組研究結果表明，差異基因最顯著富集的途徑是植物與病原菌互作，并認為Ca2+信號轉導途徑在丁家壩抗寒過程中發(fā)揮重要作用。本研究中，由真菌幾丁質、Ca2+、細菌鞭毛等不同的PAMPs激發(fā)的反應中，CDPK和CaM/CML的轉錄本顯著上調。CDPK基因上調可促進ROS的積累，是植物早期防衛(wèi)反應的標志之一[35]；CaM/CML的上調加速NO的產(chǎn)生，促發(fā)過敏性壞死反應和抗病反應[36]。此外，F(xiàn)LS2受體識別來自細菌鞭毛蛋白的flg22肽，激發(fā)丁家壩的早期防衛(wèi)反應后進一步激發(fā)MAPK級聯(lián)反應，將細胞表面的信號放大并傳遞到細胞內，引起細胞內的轉錄及代謝水平變化，從而響應低溫脅迫。這與田玉珍等[37]對天汪一號蘋果花芽響應不同時間低溫脅迫的研究結果相同，該研究發(fā)現(xiàn)差異基因最顯著富集的途徑為植物與病原菌互作，并表明Ca2+在誘導CaM/CML上調表達參與低溫響應的同時也激發(fā)MAPK級聯(lián)通路中的WRKY22、WRKY33和WRKY25等基因參與低溫響應。本研究中，MAPK級聯(lián)反應的功能酶基因MPK3/6以及WRKY22、WRKY24和WRKY33等相關基因在冷處理下的丁家壩葉片中具有較高表達，參與丁家壩對低溫的響應過程。WRKY22、WRKY24和WRKY33分別在菜心和篤斯越橘的研究中被證實響應低溫脅迫[38-39]。類似的低溫響應機制在Niu等[9]研究休眠枝響應不同時間低溫脅迫的試驗中也有印證。

4 結論

本研究對低溫耐受型丁家壩和敏感型加納巖在冷處理下的生理生化水平變化和轉錄組進行比較分析，發(fā)現(xiàn)植物與病原菌互作和MAPK信號通路以及CDPK、CaM、CML、MPK3/6、WRKY22、WRKY24和WRKY33等相關基因在丁家壩抗寒過程中發(fā)揮重要作用。研究結果有助于了解丁家壩響應冷脅迫的分子機理，為桃耐冷品種的培育提供理論依據(jù)。