基于RNA-seq的番茄青枯病響應基因鑒定與表達分析

史建磊，熊自立，蘇世聞，王克磊，宰文珊

(1.溫州科技職業學院，浙江 溫州 325006；2.福建農林大學 農學院，福建 福州 350002)

青枯病是一種毀滅性的土傳維管束病害，又稱細菌性枯萎病，病原為青枯雷爾氏菌(Ralstoniasolanacearum，Rs)，具有明顯的生理分化和寄主多樣性。作為世界上最重要的植物病原細菌之一，青枯菌可從植物傷口或次生根的根冠侵入，一定條件下大量繁殖并快速傳播，破壞并阻塞維管束，影響水分運輸，從而導致植株萎蔫死亡。該病主要發生在熱帶亞熱帶氣候濕潤、土壤呈酸性的地區，在我國主要發生在長江流域及以南的廣大地區，但受全球氣候變暖等的影響，有向北蔓延的趨勢，是造成多種作物減產、甚至絕收的主要原因，特別是對茄科作物番茄的危害尤為嚴重。番茄(Solanumlycopersicum)是全球重要蔬菜，但遺傳基礎狹窄，經常受到各種病害的影響。由于連作重茬、土壤富集病菌和高溫高濕氣候，我國番茄產區青枯病日趨嚴重，已成為番茄生產的主要障礙。

隨著二代測序技術(Next generation sequencing，NGS)的快速發展和廣泛應用，轉錄組測序(RNA-seq)已成為研究基因表達變化和表達調控網絡的有效手段，為抗病候選基因的鑒定和分子標記開發提供了基礎。目前，RNA-seq技術已在番茄多種生物或非生物脅迫研究中得到廣泛應用[1-8]。關于青枯病響應基因的轉錄組鑒定，Ishihara等[9]利用基因芯片技術對莖刺法接種Rs的番茄莖進行轉錄組分析，發現致病、激素信號和木質素生物合成相關的140個基因在抗病番茄LS-89中上調表達，β-1，3-葡聚糖酶基因尤為明顯；French等[10]對傷根法接種Rs的番茄根進行了轉錄組分析，發現抗病番茄Hawaii7996根部防御基因更早更強被激活，同時生長素信號途徑被抑制；Zuluaga等[11]對傷根法接種Rs的野生馬鈴薯根進行轉錄組分析，在抗感材料中分別發現221，644個差異表達基因，前者誘導產生大量新的轉錄本，后者激素相關基因更突出；Chen等[12]對傷根法接種Rs的茄子莖進行轉錄組分析，發現抗病材料分別有1 137，9 048個基因上調或下調表達，而感病材料對應數目是6 087，5 832個；謝冰等[13]葉轉錄組分析發現，Y87作為砧木能夠提高接穗HD PAL活性和木質素、多酚及病程相關蛋白等多種抗性基因表達量，從而提高煙草青枯病抗性。

面對各種病原微生物傷害，植物在漫長的進化過程中形成了完善的自我保護機制，抗病基因(R基因)在識別和抵御病原菌入侵方面具有重要作用。盡管前人已在番茄青枯病響應方面進行了相關研究，但不同材料、不同接種方法、不同取樣部位等獲得的研究結果不盡一致，有必要拓展和豐富這一課題。

本研究利用高通量測序技術對傷根法接種Rs后的抗感番茄自交系地上部進行轉錄組分析，結合基因功能注釋，挖掘青枯病響應相關基因并進行RT-qPCR驗證，以期為候選基因抗源的有效利用提供理論指導。

1 材料和方法

1.1 試驗材料

番茄AH13112111為溫州農業科學院番茄課題組2013年春從荷蘭瑞克斯旺公司引進品種愛好與自育AH1雜交選育的高代自交系，黃色小果型，苗期人工接種鑒定為高抗青枯病；G149351121為以色列海澤拉公司同步引進的品種3951與自育G14雜交選育的高代自交系，紅色普通型，對青枯病表現為高感。青枯菌株(生理小種1生化變種Ⅲ)從田間發病番茄植株上自行分離。

1.2 試驗方法

1.2.1 建庫與測序 番茄幼苗長到4～6片真葉時用傷根法接種青枯菌，接后48 h植株萎蔫前分別取幼嫩組織，錫箔紙包裹，立即置于液氮中，再放入-80 ℃超低溫冰箱中保存備用[14]。設置3次生物學重復。參照植物TRIzol總RNA提取試劑盒(Sangon Biotech)說明書提取RNA。用帶有Oligo(dT)的磁珠富集mRNA；加入Fragmentation Buffer將mRNA進行隨機打斷；以mRNA為模板，用六堿基隨機引物合成第1條cDNA鏈，然后加入緩沖液、dNTPs、RNase H和DNA Polymerase Ⅰ合成第2條cDNA鏈，利用AMPure XP beads純化cDNA；純化的雙鏈cDNA再進行末端修復、加A尾并連接測序接頭，然后用AMPure XP beads進行片段大小選擇；最后通過PCR富集得到cDNA文庫。文庫質控合格后進行Illumina轉錄組測序(Biomarker Technologies)。

1.2.2 數據質量控制 去除原始數據中帶有接頭序列和低質量的reads，得到Clean reads。計算Clean reads的Q30值和GC含量。從茄科基因組數據庫(Solanaceae Genomics Network，https：//solgenomics.net/organism/Solanum_lycopersicum/genome)下載番茄參考基因組(ITAG4.0)的序列文件和注釋文件。使用HISAT2[15]進行測序序列的比對。比對分析完成后利用StringTie[16]對符合要求的reads進行組裝和定量。基于所選參考基因組序列，使用StringTie對mapped reads進行拼接，并與原有的基因組注釋信息進行比較，尋找原來未被注釋的轉錄區，發掘該物種的新轉錄本和新基因。

1.2.3 差異基因篩選 以差異倍數(Fold Change，FC)≥2且錯誤發現率(False Discovery Rate，FDR)<0.01作為標準，利用DESeq2[17]進行樣品組間的差異表達分析，篩選差異表達基因。通過Blast比對工具，將測序基因與蛋白相鄰類聚簇數據庫(Cluster of orthologous groups of proteins，COG)[18]、基因本體論數據庫(Gene ontology consortium，GO)[19]、東京基因與基因組百科全書(Kyoto encyclopedia of genes and genomes，KEGG)[20]、蛋白質真核同源數據庫(Eukaryotic orthologous groups，KOG)[21]、非冗余蛋白數據庫(Non-redundant protein datebase，NR)[22]、Pfam[23]和蛋白質序列數據庫(Swiss prot protein datebase，Swiss-Prot)[24]進行序列比對，獲得基因的注釋信息。

1.2.4 基因定量驗證 以SlRPL2為內參基因，利用Primer Premier 5.0軟件設計定量引物(表1)。RT-qPCR反應體系包括2 μL cDNA，0.4 μL PCR primer，10 μL SYBR，7.2 μL ddH2O。擴增程序為95 ℃ 3 min；95 ℃ 5 s，60 ℃ 30 s，45個循環，整體反應在StepOne Plus型熒光定量PCR儀(ABI，Foster，CA，美國)中進行，采用2-ΔΔCt法計算基因相對表達量。用SPSS 17.0對數據進行單因素方差分析，用最小顯著極差法(LSD)進行多重比較和差異顯著性分析。

表1 部分RT-qPCR引物Tab.1 Part of the primers used for RT-qPCR

2 結果與分析

2.1 測序數據質量評估

經檢測，樣本RNA的OD260/280均為2.2，OD260/230≥1.8，28S/18S ≥2.0，RIN≥8，表明所提RNA滿足后續分析需求。抗感番茄接種前后12個樣本原始測序數據去除接頭序列和低質量序列，共獲得75.02 Gb Clean reads，各樣本凈堿基數均不小于5.73 Gb，GC含量在43.27%～44.17%，Q30堿基百分比均超過了94%(表2)。將各樣本的Clean reads與指定參考基因組進行序列比對，比對效率為95.00%～96.31%，且比對到唯一位置的reads均超過了90%(表3)。基于基因表達量，重復樣本間相關系數在0.947～0.993(表4)，生物學重復良好。

表2 番茄各樣本RNA測序數據統計Tab.2 RNA sequencing data statistics of tomato libraries

表3 番茄各樣本Clean reads與參考基因組比對Tab.3 Comparing Clean reads onto the reference genome in each tomato library

表4 生物學重復樣本的相關性Tab.4 Correlation between biological replicates

2.2 差異表達基因篩選

本研究測序共獲得35 371個基因，其中新基因1 296個。經log2FC和FDR過濾，在抗感番茄中分別獲得970，695個DEGs，二者重復基因353個，DEG總數1 312個(新基因32個)，占基因總數的3.71%。其中，上調基因分別為457，450個，重復214個，總計693個；下調基因分別為513，245個，重復137個，總計621個(表5)。2種材料接種后均上調或下調的基因分別為214，137個；抗病材料上調而感病材料下調的基因數目是0，反之是2(Solyc04g074440.1和Solyc12g096780.2)；抗病材料上調且感病材料不變的基因是198個，反之是216個；抗病材料下調且感病材料不變的基因是316個，反之是104個，后6種表達模式基因總數為836個(新基因18個)(表6)。

表5 差異表達基因數目Tab.5 Number of differentially expressed genes

2.3 差異表達基因注釋

本研究注釋基因共31 903個，注釋率為90.20%。其中，抗感番茄DEG在GO和KEGG等八大數據庫的注釋數目分別為951，687個，注釋率均超過了98%；總注釋DEG是1 290個，注釋率是98.32%。這些注釋基因包括4個NBS類抗病基因、6個植物-病原互作基因、11個植物激素信號轉導基因、22個防御反應基因、32個蛋白激酶、65個轉錄因子及其他多個重要功能因子。

表6 差異表達基因的模式分布Tab.6 Pattern distribution of differentially expressed genes

1 312個DEGs GO注釋基因為973個，注釋率74.16%。將這些基因劃分成47個功能組，根據功能組的不同，又分為生物過程(Biological process，BP，19個亞類)、細胞組分(Cellular component，CC，15個亞類)和分子功能(Molecular function，MF，13個亞類)三大類(圖1)。在生物過程中，主要涉及代謝過程、細胞過程、單生物體過程、生物調控和刺激響應等亞類；對于細胞組分，主要涉及細胞、細胞部分、膜、細胞器、膜部分、細胞器部分和大分子復合物等亞類；結合催化活性在分子功能中非常顯著，其次是轉運活性和核酸結合轉錄因子活性，再次是電子載體活動、抗氧化活性、結構分子活性和信號傳導活性。

1.代謝過程；2.細胞過程；3.單生物體過程；4.生物調控；5.刺激響應；6.定位；7.細胞組成與發生；8.信號；9.發育過程 ；10.多細胞生物體過程；11.生殖過程；12.生物群體過程；13.生長；14.繁殖；15.免疫系統過程；16.細胞死亡；17.節律過程；18.移動；19.生物粘附；20.細胞；21.細胞部分；22.膜；23.細胞器；24.膜部分；25.細胞器部分；26.大分子復合物；27.膜內腔；28.胞間區；29.細胞連接；30.共質體；31.胞間區部分；32.顆粒；33.顆粒部分；34.類核；35.結合；36.催化活性；37.轉運活性；38.核酸結合轉錄因子活性；39.結構分子活性；40.電子載體活性；41.分子傳導活性；42.信號傳導活性；43.抗氧化活性；44.轉錄因子活性-蛋白結合；45.養分存儲活性；46.蛋白質標簽；47.金屬伴侶活性。

本研究209個(15.93%)DEGs被注釋到88個代謝通路中，分為細胞過程、環境信息處理、遺傳信息處理、代謝和生物系統五大類代謝途徑(圖2-A)。其中，植物激素信號轉導和苯丙烷生物合成相關基因分別為16個，淀粉和蔗糖代謝相關基因13個，植病互作相關基因11個，同時注釋還包括核酸、糖、脂肪、蛋白、萜類和黃酮類等的生物合成與代謝相關基因。依據q-value和富集因子大小，取前20個最顯著富集的代謝途徑(圖2-B)，主要包括DNA復制(9個DEGs)、二萜合成(7個DEGs)、組氨酸代謝(6個DEGs)、苯丙烷生物合成(16個DEGs)、類固醇生物合成(5個DEGs)、類黃酮生物合成(6個DEGs)、萜類骨架合成(5個DEGs)、谷胱甘肽代謝(9個DEGs)。另外，植病互作(rank 38)、淀粉和蔗糖代謝(rank 46)、植物激素信號轉導(rank 47)相關基因分別為11，13和16個。

圖2 差異表達基因KEGG注釋(A)與富集(B)Fig.2 KEGG annotation(A)and enrichment(B)of differentially expressed genes

植病互作通路在抗感番茄接種后顯著富集，抗病番茄KEGG分析該通路中有7個基因差異表達，感病番茄該通路中有6個基因差異表達，2個材料中共有差異表達基因11個(圖3)。該通路中抗病番茄接種后CNGCs(Solyc02g086980.3和Solyc05g050380.4)、EIX1/2(Solyc07g008620.1)、FLS2(Solyc02g031790.3、Solyc02g068820.3和Solyc02g070890.3)和WRKY25/33(Solyc09g014990.4)上調表達；感病番茄接種后CaM/CML(Solyc04g058170.1)、FLS2(Solyc02g068820.3)、PIK1(Solyc06g075550.3)和RPM1(Solyc01g073985.1)上調表達，CNGCs(Solyc05g050380.4和Solyc09g007840.3)既有上調表達又有下調表達。

圖3 植病互作通路差異表達基因KEGG富集Fig.3 KEGG enrichment of differentially expressed genes in the plant-pathogen interaction pathway

2.4 差異表達基因啟動子分析

基于基因功能注釋，選取74個DEGs進行啟動子分析。結果發現這些基因啟動子主要包括轉錄起始(TATA-box和CAAT-box)、防御與脅迫(W box、TC-rich repeats和AT-rich sequence)、干旱(MBS)和低溫(LTR)、MeJA(CGTCA-motif和TGACG-motif)、ABA(ABRE)、SA(TCA-element)、GA(TATC-box、Pbox和GARE-motif)和AUX(TGA-element、TGA-box、AuxRR-core和AuxRE)、光響應及厭氧誘導等元件(圖4)。其中，29個DEGs含42個防御與脅迫響應相關元件。

2.5 差異表達基因的誘導表達

選取上述50個DEGs進行RT-qPCR驗證。與轉錄組測序相比，28個基因青枯菌接種前后的表達變化趨勢一致，結合差異顯著性分析，11個基因表達調控模式一致，其中6個基因完全一致(圖5)。番茄接種青枯菌后，Solyc02g086980.3和Solyc04g011670.3在不同材料中相對表達量均極顯著上調(P<0.01)，抗病材料中分別提高2.62，2.36倍，感病材料分別提高1.06，1.34倍，為正調控基因；其他4個基因(Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1)表達量在感病材料中極顯著上調(P<0.01)，分別提高1.16，1.19，1.63，0.87倍，但在抗病材料中變化不顯著(P>0.05)，為負調控基因。推測Solyc02g086980.3和Solyc04g011670.3可能參與抗病反應，Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1可能參與感病反應。

A.已知基因；B.新基因。A.Known genes；B.New genes.

3 結論與討論

番茄青枯病是由青枯菌引起的一種土傳維管束病害，獲得抗病相關基因并理解其作用機制非常重要，而RNA-seq技術為這一問題的解決提供了有效途徑。本研究測序共獲得35 371個基因，其中1 296個新基因，1 312個DEGs，836個材料特異DEGs，并對其進行了功能注釋和RT-qPCR驗證，這為目的基因篩選、功能驗證和相關機制解析提供了基礎。

3.1 差異表達基因的篩選與注釋

本研究測序共獲得75.02 Gb高質量數據，各樣本均達到5.73 Gb，Q30堿基占比均大于94%，過濾數據與番茄參考基因組比對率均不小于95%，說明樣品不存在污染且參考基因組選擇合適，所有樣本測序質量良好，數據可靠，可用于后續分析。獲得的35 371個基因中注釋基因數目是31 903個，略多于Heinz 1706的測序基因數30 855個[25]，可能源自材料的不同和技術手段的改進。基于所選參考基因組序列，共發掘到1 296個新基因，說明RNA-seq技術在新基因挖掘上的可行性。生物學重復可有效消除基因表達在不同個體間的可變性，本研究重復樣本皮爾遜相關系數大于0.9，表明生物學重復良好。抗感番茄接種前后DEG數目多于Ishihara等[9]的鑒定結果，但少于French等[10]的鑒定結果，主要歸因于接種方法、取樣部位和基因表達計算的不同。有研究表明，地上部早期對青枯病的防御反應遠沒有地下部強烈，但萎蔫病癥出現后卻有大量基因差異表達[26]。抗病番茄DEG要多于感病番茄，尤其是下調基因，說明抗病性主要與病菌侵染后某些基因下調表達有關。不同材料間DEG數目少于接種前后的1 312個，說明差異表達主要來自接種。感病番茄接種前后DEG少于接種前抗感材料間的DEG，且不同材料接種后DEG變少，說明感病番茄對青枯菌侵染的響應不如抗病番茄劇烈，這和前人研究結果一致[9-10]。抗感材料特異表達的836個DEGs可作為病害響應基因和候選目的基因的篩選基礎。

R和S分別表示抗病和感病番茄。不同大小寫字母分別表示0.01和0.05水平下的差異顯著性。R and S represent resistant and susceptible tomato lines，respectively.Different upper and lowercase letters indicate statistically significant differences at the 0.01 and 0.05 level，respectively.

本研究通過抗感材料轉錄組分析表明，病菌侵染會啟動植物相關基因的表達及多條代謝通路。對這些DEG進行GO富集分析，抗病反應主要發生在膜和細胞內部，同時也發生在細胞連接和細胞外部區域等部位，通過調節催化活性、結合、轉運活性、轉錄因子活性、信號轉導活性以及抗氧化活性等分子功能，引起以細胞過程、代謝過程、單生物體過程和生物調控為主，定位、刺激響應、細胞成分組成、信號、免疫系統過程等一系列生物過程的協同變化，抵抗病原菌的入侵。通過KEGG分析表明，DEG在植物激素信號轉導、淀粉和蔗糖代謝、DNA復制、嘧啶代謝、嘌呤代謝、植病互作、谷胱甘肽代謝、苯丙烷生物合成、氨基酸合成代謝、玉米素生物合成、泛素化蛋白水解、類黃酮生物合成、萜類物質合成、異喹啉生物堿生物合成、芪類化合物合成、ABC轉運蛋白、內吞作用、角質蠟質生物合成、油菜素甾醇生物合成、VB1/6和VC代謝和生理節律等多個代謝通路富集，說明這些代謝基因可能參與了抗病調控。本研究選取番茄幼苗地上部進行分析，其對青枯病的防御反應既與地下部部分重疊，又具有自身獨特的代謝途徑[10]。

基于基因功能注釋，50個DEGs用來RT-qPCR驗證(這些基因普遍具有轉錄起始和脅迫響應相關元件)，結果發現，28個基因在番茄接種青枯菌前后表達變化趨勢與轉錄組數據一致。推測Solyc02g086980.3和Solyc04g011670.3可能參與抗病反應，Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1可能參與感病反應。2種策略的相關性不強，基因表達結果不完全一致，主要歸因于計算原理和表示方法的不同。本研究結合2種方法和統計檢驗，篩選到6個病害響應基因。基因表達受植物生長發育與環境脅迫響應的動態平衡調控，一方面通過R基因的過量表達抵御外界不良刺激，另一方面通過非編碼RNA調控和可變剪接等降低基因表達量[27]。然而，本研究也發現多個感病基因(S基因)，與R基因堆疊相比，破壞單個S基因是實現作物廣譜和持久抗性的更有效策略[28]。

3.2 抗病相關分子機制初步解析

上述這些物質有些可能對抗性有直接影響，有些則可能進一步形成抗性物質。其中植物-病原互作代謝通路包括PTI(Pattern-triggered immunity)和ETI(Effector-triggered immunity)，當植物受到胞外病原相關分子模式(Pathogen-associated molecular pattern，PAMP)刺激時，Ca2+濃度增加，激活鈣離子蛋白激酶(CDPK)和鈣調蛋白(CaM/CML)，將信號傳導至NADPH氧化酶活性氧中間體(RBOH)和一氧化氮合成酶(NOS)，激活活性氧(ROS)和一氧化氮(NO)的產生，誘導超敏反應(Hypersensitive response，HR)、細胞壁加固和氣孔關閉[29]。當細菌鞭毛蛋白被受體激酶(FLS2)識別，通過胞吞作用將信號傳遞到細胞內，激活下游的絲裂原活化蛋白激酶(MEKK)MKK1/MKK2→MPK4或MKK4/MKK5→MPK3MPK6。信號傳入細胞核后，WRKY22/29激活防御基因或WRKY25/33抑制相關基因[30]。本研究中，環腺苷酸門控通道(CNGCs)基因Solyc02g086980.3在抗病番茄接種后上調表達，Solyc05g050380.4在抗感番茄中均上調，Solyc09g007840.3在感病番茄中下調表達，Ca2+濃度發生變化，激活CaM/CML基因Solyc04g058170.1在感病番茄中的表達。青枯菌侵染后，FLS2基因Solyc02g031790.3和Solyc02g070890.3在抗病番茄中均上調表達，Solyc02g068820.3在抗感番茄中均上調表達，激活下游WRKY25/33基因Solyc09g014990.4在抗病番茄中的表達。植物通過ETI防御機制特異的抵抗病原菌入侵，該通路中RPM1和RPS2是2個能與RIN4基因結合的R基因。當RIN4被磷酸化時，與R基因結合，觸發HR反應。蛋白激酶RPS5與PBS1結合，也可觸發免疫反應[29]。本研究中，感病番茄接種青枯菌后，RPM1基因Solyc01g073985.1和PIK1基因Solyc06g075550.3上調表達。本研究初步探索了抗感番茄在接種青枯菌后的轉錄組動態，發現相關抗性涉及一個非常復雜的網絡調控系統。

本研究轉錄因子(Transcription factor，TF)主要涉及AP2/ERF、bHLH、bZIP、Dof、MYB、NAC和WRKY等。這些TFs能激活下游生物脅迫響應相關基因的表達，從而提高茄科作物青枯病抗性[31-40]。32個差異蛋白激酶主要為受體類激酶，其能通過磷酸化在細胞信號轉導中發揮作用，從而提高茄科植物抗病性，如AhRLK1、CaCDPK15、CaLRR51和CaRop1等能夠與其他基因互作正/負調控植物青枯病抗性[41-43]。4個差異NBS基因分正/負調控。Deslandes等[44]克隆獲得了首個抗青枯病基因RRS1。另外，在擬南芥突變體SLH1中發現了另一個編碼TIR-NBS-LRR-WRKY類抗病蛋白的基因RPS4。NBS類抗病基因AhRRS5通過多種信號途徑交聯參與煙草青枯病防御反應[45]。Ben等[46]在苜蓿中發現青枯病抗性位點MtQRRS1，包括15個TIR/NBS類抗病基因。面對病原菌入侵，植物會通過激素信號分子激活抗病基因的表達。本研究植物激素信號轉導涉及生長素、細胞分裂素(CTK)、脫落酸(ABA)、水楊酸(SA)和茉莉酸(JA)等通路，SA通路基因通過參與苯丙氨酸代謝過程正調控抗病性，而JA通路基因通過參與α-亞麻酸代謝負調控抗病性。Sánchez-Vallet等[47]和Lim等[48]報道ABA信號途徑參與植物抗病免疫，French等[10]也證明生長素信號途徑在番茄青枯病抗性上的重要性。

本研究其他次生代謝物和功能分子也能參與植物抗病過程。一些MYB可以通過調控靶基因的表達來調控植物體內次生代謝物質的合成，保護植物抵御生物脅迫，如AtMYB15調控擬南芥防御誘導的木質化和基礎免疫[49]；AtMYB11組成表達增強了苯丙醇生物合成途徑中基因的表達，導致煙草和番茄植株中類黃酮和綠原酸的積累，AtMYB11和AtMYB12調控番茄和煙草中黃酮類化合物和咖啡酰奎寧酸的生物合成[50]；AtMYB46、AtMYB58和AtMYB63在次生壁形成過程中能夠激活木質素生物合成關鍵基因[51]。PpNAC1在苯丙氨酸生物合成和調控上發揮重要作用[52]。本研究篩選到的許多差異蛋白基因可能是構成抗性的重要組分，有必要對其功能和作用機制進一步解析。