甘薯塊根響應長喙殼菌（Ceratocystis fimbriata）侵染的轉錄組測序分析

楊冬靜，馬居奎，唐 偉，陳晶偉，高方園，張成玲，孫厚俊，謝逸萍

（江蘇徐淮地區徐州農業科學研究所/農業農村部 甘薯生物學與遺傳育種重點實驗室，江蘇 徐州 221131）

由甘薯長喙殼菌（Ceratocystis fimbriata Ellis et Halsted）引起的甘薯黑斑病（black rot）是甘薯的主要病害之一，在甘薯栽培和貯藏期均可發生；甘薯長喙殼菌分布廣泛，種群分化復雜，可侵染甘薯幼苗莖基部和塊根，導致死苗、爛床以及爛窖［1］。我國每年由甘薯黑斑病造成的產量損失率達5%～10%，發病嚴重時可達20%～50%［2］。目前在生產上可采用溫湯浸種來防治甘薯黑斑病，但是該方式操作復雜，不適宜推廣應用。本課題組歷年對全國各地的甘薯品種品系進行黑斑病抗性鑒定工作，迄今尚未發現對甘薯黑斑病免疫且兼具高產、高干率特征的甘薯品種，這類甘薯品種的缺乏也成為疫區甘薯高產優質的主要限制性因素之一。在化學藥劑防治方面，張德勝等［3］研究發現粉唑醇、百菌清等藥劑對貯藏期薯塊黑斑病的防治效果較好，但是粉唑醇殘留量高于歐盟規定的最大殘留限量標準。

植物病原菌侵染寄主是一個復雜的動態過程，涉及該過程的關鍵基因及其差異表達引起的重要代謝途徑變化從而觸發植物的免疫反應一直是人們關注的焦點［4-8］。采用轉錄組測序分析方法研究植物與病原菌的互作機理，并挖掘重要的抗病基因已有較多的報道［9-13］。筆者以甘薯塊根為研究對象，基于轉錄組測序技術，研究了黑斑病菌侵染過程中甘薯塊根在轉錄組水平上的變化，深入分析了與甘薯響應黑斑病菌侵染過程密切相關的差異基因和代謝途徑，并且初步探究了甘薯與黑斑病菌互作過程的抗病分子調控網絡，以期挖掘出甘薯與黑斑病菌互作過程中的關鍵抗病調控因子，為后續研究提供理論依據。

1 材料與方法

1.1 供試材料

供試甘薯品種為相對抗黑斑病品種南京92（R）和相對感病品種徐薯18（S），選取這2個品種的健康、大小均勻的薯塊，每個品種取3塊，洗凈晾干后，采用75%的酒精噴霧殺菌，自然晾干后備用。用滅菌水洗滌在PDA平板上培養5～7 d的甘薯黑斑病菌，制備成分生孢子懸浮液，再采用滅菌紗布過濾掉洗下的菌絲，在顯微鏡下調節分生孢子懸浮液的濃度為105CFU/mL，備用。

1.2 接種處理方法

用滅菌的針頭蘸取分生孢子懸浮液，再針刺薯塊，每個薯塊均勻接種3排，每排接種5個點；分別于接種0、1、3、7 d后取少量發病部位的薯塊，將抗、感品種的各薯塊樣品分別命名為R1、R2、R3、R4，S1、S2、S3、S4；將各樣品用液氮速凍，然后置于-70 ℃冰箱中保存，待各樣品收集完成后用干冰寄送至廣州基迪奧生物科技有限公司進行測序。

1.3 RNA的提取、檢測和文庫制備

采用TIANGEN多糖多酚植物總RNA提取試劑盒［天根生化科技（北京）有限公司，貨號DP441］提取樣品的RNA，具體操作步驟參考該試劑盒的說明書。為了保證測序的質量，通過以下嚴格的質控檢測控制RNA的提取質量：（1）瓊脂糖凝膠電泳檢測，分析樣品RNA的完整性、是否降解，以及其是否存在DNA污染；（2） NanoPhotometer spectrophotometer檢測，檢測RNA的純度（OD260/OD280及OD260/OD230比值）；（3）Agilent 2100 bioanalyzer檢測，利用測定得到的RNA完整值［RNA integrity number （RIN），即RNA分子完整數］來確定RNA的質量，RIN的取值范圍為0～10，其值越大表示RNA 的質量越好，越完整。

文庫制備方法：用帶有Oligo（dT）的磁珠富集具有polyA尾巴的真核mRNA后，再用超聲波把mRNA打斷。以片段化的mRNA為模板，以隨機寡核苷酸為引物，在M-MuLV逆轉錄酶體系中合成cDNA第一條鏈；隨后用RNaseH降解RNA鏈，并在DNA polymerase I體系下，以dNTPs為原料合成cDNA第二條鏈。純化后的雙鏈cDNA經過末端修復、加A尾并連接測序接頭，用AMPure XP beads篩選200 bp左右的cDNA，進行PCR擴增，并再次使用AMPure XP beads純化PCR產物，最終獲得文庫。

1.4 數據質控

為了保證數據質量，在信息分析前對原始數據進行過濾，從而減少無效數據帶來的分析干擾。首先，對下機的raw reads利用fastp［14］進行質控，過濾掉低質量數據，得到clean reads，reads過濾的步驟如下：去除含adapter的reads；去除含N比例大于10%的reads；去除全部都是A堿基的reads；去除低質量reads（質量值Q≤20的堿基數占整條read的50%以上）。原始數據經過過濾后，分析其堿基的組成及質量分布，以直觀展示數據的質量情況。堿基組成越平衡，則質量越高，后續分析越準確。

1.5 序列比對分析

1.5.1 核糖體比對 在正常生物體內，mRNA只占總RNA的很少部分，其他絕大多數都是核糖體RNA（rRNA）。用帶有Oligo（dT）的磁珠富集帶有polyA尾的mRNA過程中，由于受樣品質量和物種的影響，仍可能會有部分rRNA殘留。筆者使用短reads比對工具bowtie2［15］將clean reads比對到該物種的核糖體數據庫，在不允許錯配情況下去除比對上核糖體的reads，將保留下來的unmapped reads用于后續轉錄組分析。

1.5.2 參考比對 本研究的比對選用甘薯品種泰中6號的參考基因組（http：//public-genomes-ngs.molgen.mpg.de/SweetPotato/DOWNLOADS/）。由 于選擇比對到基因組，因此所選用的比對軟件要求能夠進行splice，否則會丟失很多有效的junction reads，這是因為可變剪接和polyA尾在真核生物中普遍存在，跨越可變剪接位點或ployA尾邊緣的reads普遍存在。HISAT2軟件［16］采用全局和局部搜索的方法能夠有效地比對到RNASeq測序數據中的spliced reads，是目前比對率最高且最準確的比對軟件。因此本研究采用HISAT2軟件開展基于參考基因組的比對分析。

1.6 表達量統計

基因表達量的準確性非常依賴轉錄本重構結果的完善程度。Stringtie軟件［17］能夠有效利用轉錄本的豐度信息，組裝出更多的轉錄本，組裝結果也更為準確。根據HISAT2的比對結果，利用Stringtie重構轉錄本，并計算每個樣本中所有基因的表達量。表達量以原始reads count和FPKM展示，其中原始reads count表示該轉錄本所包含的reads數目，但受測序量和基因長度的影響，原始reads count不利于樣本間的差異基因比較。為了保證后續分析的準確性，本文先對測序深度進行校正，再對基因或轉錄本的長度進行校正，在獲得基因的FPKM值后，再進行后續分析。

1.7 樣本關系的主成分分析

基于表達量信息，利用R（http：//www.r-proje-ct.org/）開展主成分分析（Principal Component Analysis，PCA），從而利用降維的思路研究樣本間的距離關系。在所得的分析結果中，樣品組成越相似，則在PCA圖中的距離越近；來自不同有效處理的樣品往往表現出各自聚集的分布。

1.8 差異表達基因（DEGs）分析

基因差異表達分析的輸入數據為基因表達水平分析中得到的reads count數據，使用DESeq2軟件［18］進行差異表達基因的分析，主要分為3個步驟：對reads count數據進行標準化（normalization）；根據模型進行假設檢驗概率（P value）的計算；最后進行多重假設檢驗校正，得到FDR值（錯誤發現率）。本研究以FDR＜0.05且|log2FC|＞1為標準篩選差異表達基因，然后進行統計分析。

1.8.1 GO富 集 分 析 Gene Ontology（簡 稱GO）是一個國際標準化的基因功能分類體系，其提供了一套動態更新的標準詞匯表（controlled vocabulary）來全面描述生物體中基因和基因產物的屬性。GO總共有3個ontology（本體），分別描述基因的分子功能（molecular function）、細胞組分（cellular comp-onent）、參與的生物過程（biological process）。GO的基本單位是term（詞條、節點），每個term都對應1個屬性。GO功能分析一方面給出基因的GO功能分類注釋，另一方面給出基因的GO功能顯著性富集分析。在進行GO富集分析時，首先將基因向GO數據庫（http：//www.geneontology.org/）的各個term映射，并計算每個term的基因數，從而得到具有某個GO功能的基因列表及基因數目；然后應用超幾何檢驗，找出與整個基因組背景相比在基因中顯著富集的GO條目。

1.8.2 KEGG富集分析 在生物體內，不同的基因相互協調，行使其生物學功能，基于Pathway的分析有助于進一步了解基因的生物學功能。KEGG是有關Pathway的主要公共數據庫。Pathway顯著性富集分析以KEGG Pathway為單位，應用超幾何檢驗，找出與整個基因組背景相比在基因中顯著富集的Pathway。通過Pathway顯著富集分析能確定基因參與的最主要生化代謝途徑和信號轉導途徑［19］。

1.9 植物抗病相關基因（R-gene）分析

將預測的蛋白序列與相應的R-Gene數據庫（PRGdb v3.0，http：//prgdb.org）進行hmmscan 比對，篩選與植物抗病相關的基因［20］。

2 結果與分析

2.1 抗、感甘薯品種接種黑斑病菌的表型差異

抗、感甘薯品種在接種甘薯黑斑病菌分生孢子懸浮液后，分別于第1、3、7天取出薯塊，進行表皮和剖面觀察，并拍照記錄發病表型差異。從圖1可以看出：薯塊被侵染的第1天，抗、感品種的表皮和剖面均未出現明顯的發病癥狀；薯塊被侵染的第3天，抗、感品種的表皮和剖面均表現出明顯的黑斑癥狀，但在表型上并無顯著差異；薯塊被侵染的第7天，抗、感品種的表皮均表現出明顯的黑斑癥狀，病斑直徑差異不顯著，但從剖面來看，感病品種被侵染的深度明顯比抗病品種深，且感病品種與抗病品種的發病體積呈顯著差異。

圖1 2個甘薯品種塊根接種黑斑病菌后的表型

2.2 總RNA提取和質量分析

基于抗、感品種接種黑斑病菌后的表型差異，分別在不同侵染時間點取樣并提取RNA，瓊脂糖凝膠電泳檢測結果顯示，本研究中24個樣品的RNA完整性較好，未發現明顯降解，且不存在DNA污染。此外，采用RNA純度和濃度檢測的結果也顯示所提取的RNA質量完全符合建庫和轉錄組測序的要求。

2.3 測序質控

對原始數據進行質量控制以后，獲得的clean reads數目統計結果如表1所示，各樣品的clean reads數目均在62529524～88661116。堿基測序錯誤率分析顯示：Q20堿基百分比（錯誤率＜1.0%）分布介于97.51%～98.33%，Q30堿基百分比（錯誤率＜0.1%）分布介于92.87%～95.08%，GC含量介于56.40%～77.77%。以上結果均表明測序數據的質量較高，過濾后的有效序列reads可用于后續生物信息學分析。

表1 轉錄組測序數據統計分析結果

2.4 主成分分析

各樣本基因表達值（FPKM）的主成分分析結果如圖2所示，在相同時間點組內樣品的表達基因在每個生物學重復中都比較接近（S4除外），抗病品種的第0天對照組R1與第1、3和7天的處理組R2、R3、R4明顯比較遠，且處理時間越長，各處理組與R1的距離越遠。抗病處理組R2雖然在癥狀上表現不明顯，但是隨著黑斑病菌的侵染其差異表達基因的類別和表達量都有了較大的變化，故而與對照組R1差異較大。隨著侵染時間的增加，差異表達基因的類別和表達量變化更多，但與R1相比R4、R3已經很接近，而侵染7 d時的癥狀表型差異更明顯，這表明差異基因的表達和系列防衛反應發生在侵染早期。感病品種4個組（S1～S4）的基因表達規律與抗病品種相似，不同的是感病品種的對照組S1與處理組S2、S3、S4之間的距離更近，表明隨著病原菌侵染時間的延長，感病品種的表達基因類別和表達量也有變化，但是變化幅度小于抗病品種。從圖3還可以看出，感病品種S1與抗病品種R1之間的距離較大，表明在兩者之間基因類別和表達量差異也非常明顯。

圖2 各樣本的主成分分析結果

2.5 抗、感品種中差異表達基因的統計分析

基于差異分析的結果，篩選FDR＜0.05并且|log2FC|＞1的基因為顯著差異基因。由圖3可見，與抗病品種對照組R1相比，處理組R2有4771個差異基因上調表達，2423個基因下調表達；處理組R3有6405個差異基因上調表達，7725個基因下調表達；處理組R4的上調表達差異基因數比R3略少，下調表達的差異基因數比R3增多；總的來看，隨著侵染時間的增加，差異基因數呈上升趨勢。感病品種中差異基因上調和下調表達的變化規律與抗病品種相似，也是隨著侵染時間的增加，差異表達的基因數增多，且均在侵染第3天的增幅最大。比較同一取樣時間點抗病品種和感病品種的基因表達情況，發現隨著侵染時間的增加，抗、感品種的差異表達基因數逐漸減少，差異表達基因數從0 d的12124個逐漸下降到7 d的3468個。

圖3 不同處理組間差異表達基因統計結果

抗病品種3個處理組與對照組對比的維恩圖（圖4A）表明，在第1、3、7天均出現差異表達的基因有4955個。感病品種3個處理組與對照組對比的維恩圖（圖4B）顯示，在第1、3、7天均出現差異表達的基因有5644個。從圖4C可以看出，參與響應黑斑病菌侵染的核心差異基因有2524個。從相同侵染時間點來分析抗、感品種之間的差異，發現在R1vs S1中，R1上調表達的基因數為5290，下調表達的基因數為6834，表明在侵染之前抗、感品種之間就存在了差異，且隨著侵染時間的增加，抗、感品種之間上調和下調表達的基因數均逐漸減少；在R1vs S1、R2vs S2、R3vs S3和R4vs S4這4個對比組中有1057個共有的差異表達基因（圖4D），這些共有的DEGs與侵染時間點沒有關系，是抗、感品種之間本身就存在的差異，與黑斑病菌侵染與否無關。

圖4 抗、感甘薯品種差異表達基因的維恩分析結果

2.6 差異基因GO功能的顯著性富集分析

選取侵染3 d的處理組與對照組相比,進行差異基因GO（Gene Ontology）功能的顯著性富集分析，篩選了顯著富集前20的GO條目，GO富集圈結果如圖5所示。抗病品種被侵染3 d的差異基因中，可明顯觀察到隸屬于生物過程的是單一生物代謝過程（GO：0044710）、前體代謝產物和能量的產生（GO：0006091）和對金屬離子的響應（GO：0010038）；分子功能的GO富集分析顯示，主要是氧化還原酶活性（GO：0016491）和四吡咯結合（GO：0046906）；細胞組分的GO富集分析顯示，主要包括細胞膜（GO：0016020）、膜的固有成分（GO：0031224）和質體（GO：0009536）等15個條目（圖5A）。對感病品種侵染3 d的差異基因進行GO富集分析，結果（圖5B）顯示，在被富集到的前20的GO條目中，有4個條目（GO：0044710等）是參與生物過程的；有3個條目（GO：0016491等）是參與分子功能的；有13個條目（GO：0009536等）均是參與細胞組分相關功能的。對比圖5A和圖5B可以看出，隸屬于生物過程的單一生物代謝過程（GO：0044710）、隸屬于分子功能的氧化還原酶活性（GO：0016491）和四吡咯結合（GO：0046906），以及隸屬于細胞組分的質體（GO：0009536）和細胞外圍（GO：00071944）等條目在抗、感品種中均被富集到，表明當長喙殼菌侵染甘薯塊根時抗、感品種在這些生物過程中均有響應。

圖5 抗病及感病品種中對照組和處理組的差異GO富集圈分析結果

此外，還有一些GO條目是在抗病品種中富集到而在感病品種中沒有富集到的，如隸屬于細胞組分的膜的固有成分（GO：0031224）和細胞膜部件（GO：0044425），這些過程可能是抗病品種與感病品種抵御長喙殼菌侵染的差異所在。

2.7 差異基因的KEGG通路富集分析

圖6A和圖6B分別為在抗病甘薯塊根R1vs R3對比組和感病甘薯塊根S1vs S3對比組中，差異基因的功能主要富集到的前20的KEGG通路。從中可以看出，抗、感品種中的差異基因均被富集到的KEGG通路包括代謝途徑、次級代謝產物的生物合成、碳代謝、苯丙烷生物合成、檸檬酸循環/三羧酸循環、光合作用、氧化磷酸化、谷胱甘肽代謝、氨基酸生物合成、糖酵解/糖異生、光合作用天線蛋白、光合生物的固碳作用、托烷、哌啶和吡啶生物堿的生物合成以及果糖和甘露糖代謝等14個通路；與感病品種不同，抗病品種中的差異基因還被富集到植物激素信號轉導、淀粉和糖代謝、萜類主鏈生物合成、氨基糖和核苷酸糖代謝、戊糖磷酸途徑，以及ABC轉運蛋白等信號通路；與抗病品種不同，感病品種中的差異基因還被富集到植物MAPK信號通路、半胱氨酸和蛋氨酸代謝、乙醛酸和二羧酸代謝、甘氨酸、絲氨酸和蘇氨酸代謝、二苯乙烯類、二芳基庚烷類和姜辣素生物合成以及苯丙氨酸代謝途徑。

圖6 差異表達基因顯著富集的前20的KEGG通路

進一步對抗、感品種間4個對比組中1057個共有的差異表達基因進行KEGG通路富集分析，結果如表2所示，從中可以看出，被富集到的前20的KEGG通路包括單萜生物合成、亞油酸代謝、氨基糖和核苷酸糖代謝、剪接體、異喹啉生物堿生物合成、核苷酸切除修復、其他聚糖降解、酪氨酸代謝、錯配修復、泛素介導的蛋白質水解、自噬-其他真核生物、硒化合物代謝、植物激素信號轉導、DNA復制、卟啉與葉綠素代謝、RNA降解、核糖體代謝、果糖和甘露糖代謝、氨基酸生物合成以及同源重組等通路，共涉及158個關鍵抗性相關基因。

表2 基于各時間點抗、感品種共同響應長喙殼菌侵染的差異基因的前20的KEGG通路

3 討論和結論

近年來，植物-病原菌互作研究逐漸成為熱點，以多種植物為寄主的長喙殼菌以及甘薯栽培品種泰中6號的全基因組序列均已公布，以全基因組序列為依據開展植物抗病機制和病原菌致病機制的研究越來越多，但迄今在國內外尚未見有關長喙殼菌侵染甘薯的轉錄組測序和分析的研究報道。筆者選擇對黑斑病抗性表現不同的2個甘薯品種，研究了長喙殼菌侵染的轉錄組測序，分析了抗、感品種間本身存在的差異基因及其參與抗病相關途徑的差異，所得研究結果有助于揭示不同抗性甘薯品種抗性差異的原因。本文還分析了長喙殼菌侵染過程中抗病和感病品種響應病原菌的差異，揭示了抗、感品種抵御病原菌侵染過程中存在的差異變化。在本研究設置的侵染時間1、3和7 d下，接種的甘薯薯塊分別未見明顯癥狀、已初步顯示病癥、黑斑病癥狀非常明顯，這3個時間點分別是長喙殼菌侵染甘薯的早期、中期和后期。因此，本研究在上述3個時間點取樣進行轉錄組測序，測序結果能全面反映長喙殼菌侵染狀態下抗、感甘薯品種的轉錄本差異，可以為黑斑病抗病調控網絡和抗性基因的挖掘奠定基礎。在本研究中，RNA-seq數據質控的堿基質量、測序飽和度、組內樣品重復性等質量評估結果都比較理想。在與參考基因組的比對中，各處理的clean reads數均在9340340097 bp以上，因此基于該轉錄組測序數據的后續分析和實驗結果均具有較好的真實性和可靠性。

植物響應病原菌侵染及其對病原菌的防御機制是復雜的新陳代謝過程，對病原物侵入及擴展的生理反應是以酶的催化活動來實現的［22］。超氧化物歧化酶、過氧化物酶及過氧化氫酶等是植物體內重要的保護酶類，在保護植物細胞膜免受自由基傷害的過程中起非常重要的作用［23］。本研究通過對南京92和徐薯18轉錄組測序分析發現，在甘薯長喙殼菌侵染條件下，2個材料中均有大量差異基因被富集到屬于分子功能的與氧化還原酶活性相關的GO條目（GO：0016491）中，南京92和徐薯18中分別有451和482個DEGs，表明在長喙殼菌侵染下抗、感品種體內的保護酶系統平衡狀態均被打破，從而抵御長喙殼菌的侵染，維持機體的正常運轉；同時，長喙殼菌侵染甘薯是復雜的生物過程，在此過程中差異基因被富集到屬于生物過程的單一生物代謝GO條目（GO：0044710），抗、感品種中分別有1123和1281個DEGs，從富集的DEGs數量來看，南京92受長喙殼菌侵染的影響比徐薯18小。

KEGG富集分析結果顯示，抗、感品種中的差異基因均被富集到的前20的KEGG通路包括代謝途徑、次級代謝產物的生物合成、碳代謝、苯丙烷生物合成、檸檬酸循環/三羧酸循環、光合作用、氧化磷酸化、谷胱甘肽代謝等，這表明抗、感品種在抵御長喙殼菌侵染過程中有相似的反應。另外，在抗、感品種間的R1vs S1、R2vs S2、R3vs S3和R4vs S4四個對比組中發現了1057個共有的差異表達基因，這些共有的DEGs與黑斑病菌侵染沒有關系。此外，通過R-gene分析，鑒定到221個抗性相關基因，這些抗病基因主要為激酶類基因（圖7）。

圖7 響應長喙殼菌侵染的甘薯R-gene分析結果

本研究通過對甘薯長喙殼菌侵染甘薯塊根的轉錄組測序和分析，發現了甘薯在長喙殼菌侵染條件下的大量DEGs；通過GO和KEGG富集分析，初步明確了在長喙殼菌侵染條件下甘薯的生物學過程變化以及差異基因被富集到的主要通路；此外還鑒定出了一些抗性相關基因。后續擬結合基因注釋和富集分析結果，對重點關注的相關基因進行克隆和分析，并對這些基因表達的蛋白質功能和調控網絡進行深入研究，以進一步解析甘薯響應長喙殼菌侵染的分子機制。