松樹內生細菌GD2對松材線蟲入侵寄主時轉錄組的影響

陳陽雪，趙曉佳，談家金

南京林業大學林學院，南京 210018

松材線蟲是一種危險性極高的外來有害生物，由它引起的松材線蟲病是一種毀滅性的森林病害[1]。松材線蟲病可使易感病松樹在感染松材線蟲后2～3個月內發生枯萎死亡[2]。自1905年在日本大暴發以來，該病害在世界各地被發現，截至目前該病已在歐亞和北美等多個國家或地區被發現和報道[3-4]。亞洲是松材線蟲病流行危害最嚴重的地區，其中中國是目前遭受松材線蟲病威脅最嚴重的國家，病害已經在中國10多個省造成了大面積松樹死亡，現有6 000萬hm2松林正面臨著松材線蟲病大流行的威脅[5]。經統計，該病害目前已在我國危害馬尾松(Pimusmassoniana)、油松(P.tabulaeformis)、赤松(P.densiflora)、華山松(P.armandii)、黑松(P.thunbergii)等50多種松科松屬樹種和落葉松(Larixolgensis)、云杉(Piceaasperata)、冷杉(Abiesfabri)等10多種非松屬樹種[6]，對森林生態系統造成嚴重的破壞，也造成了嚴重的經濟損失。

松材線蟲病由松材線蟲、寄主、傳播媒介、細菌真菌和環境等多個因素構成了其復雜的病害系統，因此松材線蟲的致病機制尚未完全弄清[7]。目前，國內外學者對松材線蟲病的研究工作主要以松材線蟲為病原展開[1，8-10]。近年來，已有研究顯示松材線蟲致病作用是松材線蟲與其共生細菌共同致病[11]。Kawazu等[12-15]、洪英娣等[16-17]從松材線蟲體表分離出細菌，并通過接種試驗表明松材線蟲與其伴生細菌共同導致松材線蟲病。但是細菌在這種復雜病害中的特定參與特性仍然存在爭議。Wang等[18]、He等[19]、Xue等[20]通過轉錄組測序研究了細菌對松材線蟲基因表達的影響，發現細菌對松材線蟲的致病相關基因有一定的影響。這些對于深入研究細菌在松材線蟲病中的作用提供了重要的實驗依據。菌株GD2是課題組前期從馬尾松體內分離的內生細菌，初步研究顯示其在松材線蟲致病過程中發揮著重要作用[21]。本研究擬進一步探討在松材線蟲侵入寄主建立適應性的過程中，菌株GD2對線蟲轉錄組的影響，這對全面認識松材線蟲的致病機理和松材線蟲-細菌微生態系統具有重要的理論意義。

1 材料與方法

1.1 材料與處理

供試菌株為馬尾松內生細菌GD2。供試松材線蟲AA3，分離自安徽的感病黃山松(P.taiwanensis)。蟲株由南京林業大學森林病理實驗室提供。供試松苗為生長較一致的 2 年生馬尾松(中國，廣西)，由江蘇沂北園林苗木有限公司提供。

將松材線蟲 AA3 接種到長滿新鮮灰葡萄孢的 PDA 平板上，置于 25 ℃ 培養6 d 后，采用貝爾曼漏斗法分離。將線蟲液收集到無菌的 15 mL 離心管中，3 000 r/min 離心 3 min,棄上清，用 0.1%硫酸鏈霉素處理 15 min，用無菌水沖洗離心3次。將松材線蟲調成濃度為 15 000 條/mL 的蟲懸液。

用接種環挑取1環菌株GD2的新鮮單菌落，接入含有20 mL NB培養液的錐形瓶中，27 ℃、200 r/min振蕩搖培12 h，置入無菌的 50 mL離心管中，6 500 r/min 離心 3 min，棄上清，無菌去離子水洗 3 次。將菌體用無菌水重懸制成菌懸液(濃度為5×107CFU/mL)。

線蟲 AA3 與菌株 GD2 預處理后，采用皮接法[22]接種到2年生馬尾松。設3個處理，1個對照，每處理接種13棵馬尾松：(1)AA3+GD2(T)：接種 1 mL 蟲懸液和 1 mL 菌懸液；(2)AA3(CK)：接種 1 mL 蟲懸液和 1 mL無菌去離子水；(3)菌株 GD2：接種1 mL 菌懸液和 1 mL無菌去離子水；(4)AG：接無菌去離子水共 2 mL。放于溫室 25 ℃培養，24 h后將接種無菌水(AG)、AA3(CK)、AA3+GD2(T)的馬尾松接種點部位剪下，-80 ℃保存備用。每個處理各3個生物學重復。剩余的馬尾松每 7 d 觀察接種寄主的感病癥狀，統計感病率和感病指數。根據病害分級標準和參照Yu等[23]方法計算感病指數。馬尾松的感病嚴重程度分為5個級別，Ⅰ級：植物針葉綠色，健康生長，代表值0；Ⅱ級：0～25%針葉褪綠變黃，代表值1；Ⅲ級：25%～50% 針葉褪綠變黃，枝頭彎曲，代表值2；Ⅳ級：50%～70% 針葉褪綠變黃，枝頭下垂，代表值3；Ⅴ級：75%～100% 針葉失綠褐變，植株枯萎死亡，代表值4。計算公式如下：

感病率=(感病株數/接種株樹) × 100% 感病指數= [∑(各級病樹株數 × 代表值)/ (總株數 × 發病最高級代表值)] × 100

1.2 總RNA提取、文庫構建及測序

將前面-80 ℃保存備用的樣品采用OminPlant RNA kit(DNase I)(CWBIO)法提取各處理組的總RNA，用 1% 的瓊脂糖凝膠電泳檢測，2100 Bioanalyser (Agilent) 、ND-2000 (NanoDrop Technologies)方法檢測RNA樣品的質量。RNA文庫的建立采用TruSeqTM RNA sample preparation Kit( Illumina,San Diego,CA)試劑盒。使用IlluminaNovaSeq 6000測序平臺進行高通量測序(上海美吉生物醫藥科技有限公司)。數據已提交NCBI，登錄號為PRJNA690180。

1.3 數據處理

首先將測得馬尾松(AG組)的Raw reads先進行過濾處理得到 Clean reads，再用CK組和T組的Raw reads直接比對馬尾松的Clean reads，去除比對上的reads，將剩下的reads經過濾后比對到參考基因組松材線蟲基因組(參考基因組版本：GCA_000231135.1;參考基因組來源：https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)進行序列比對。

1.4 生物信息學分析

使用FPKM(Fragments Per Kilobases per Millionreads)軟件進行定量計算獲得基因的表達量，再用軟件DESeq2進行差異基因分析，并用BH(fdr correction with Benjamini/Hochberg)進行多重檢驗。將P<0.05 且 |log2FC|≥0.585設為條件篩選差異表達基因并進行GO (Gene Ontology,http://www.geneontology.org/)、COG (clusters of orthologous groups,http://www.ncbi.nlm.nih.gov/research/cog-project/)以及KEGG (Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)分析。

1.5 實時熒光定量 PCR 驗證差異基因

選取8個差異表達基因(5個上調基因、3個下調基因)，利用SIGMA在線設計軟件(http://www.oligoarchitect.com/LoginServlet)為 8個差異表達基因設計引物(表1)，以松材線蟲肌動蛋白基因Actin(GenBank EU100952)為內參基因，每個樣品3個技術學重復，3個生物學重復，并采用 2-ΔΔCt算法分析數據。

2 結果與分析

2.1 菌株GD2影響松材線蟲病的發生

溫室接種試驗表明，菌株GD2增強了松材線蟲在馬尾松中的致病力。AA3 處理組、GD2+AA3混合處理組都使 2 年生馬尾松盆栽苗出現感病癥狀，并且所致萎蔫癥狀一致，均為針葉先褪色、黃化，進而枯萎變成紅褐色；而 GD2 處理組和對照無菌水處理組均保持健康，馬尾松松針綠色。接種后21 d，接種 AA3 的馬尾松感病率為20%，感病指數為5，而混合接種AA3和菌株GD2的馬尾松感病率達到40%，感病指數為10；接種后 42 d， AA3處理組的感病率為 60%，感病指數 52.5，而GD2+AA3 處理組的感病率為 70%，感病指數為55(表2)。

表1 松材線蟲 8 個差異表達基因的 qRT-PCR 引物 Table 1 The qRT-PCR primers for 8 differentially expressed genes of Bursaphelenchus xylophilus

表2 接種不同處理的松材線蟲的馬尾松的感病率和感病指數 Table 2 Disease incidence and disease index (DI) of Pinus massoniana inoculated with different treatments of Bursaphelenchus xylophilus

2.2 差異表達基因的篩選

對含松材線蟲的RNA樣本進行轉錄組測序，去除比對到的馬尾松reads,剩下的 reads經過濾后得到的clean reads質量達到 Q20 以上的堿基比例大于 80%，達到 Q30 以上的堿基比例大于85%，每個樣本堿基 GC 含量較為一致，且測序堿基平均錯誤率均在0.1%以下，說明測序數據良好(表 3)。將過濾后的數據與rRNA數據庫進行比對(表4)，再利用HISAT2 軟件將質控后的原始數據，即Clean reads，與參考基因組(https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)比對(表5)。將P<0.05 且|log2FC|≥0.585設為條件共篩選到143個差異表達基因，其中 63 個基因表達上調，80個基因表達下調(圖1)。這些差異基因與內生細菌GD2協同松材線蟲侵染馬尾松相關。

表3 轉錄組序列質量統計表 Table 3 The data statistics for RNA-sequencing

表4 過濾后數據與rRNA數據庫比對統計表 Table 4 Results of comparison with rRNA database

表5 序列與松材線蟲參考基因組比對結果統計 Table 5 Results of comparison with the reference genomeof Bursaphelenchus xylophilus

圖1 差異基因火山圖Fig.1 Differential gene volcano

2.3 差異基因 GO 分析

對差異表達基因進行 GO 功能注釋與分類統計。差異基因被歸類到GO 3個組分的38個功能分類中(圖2)。細胞組成(cellular component)中細胞組分(cell part)和膜組分(membrane part)最多，分別有52個基因和44個基因；分子功能(molecular functions)中結合(binding)和催化活性(catalytic activity)最多，分別有60個基因和51個基因，生物代謝過程(biological processes)中細胞進程(cellular process)和代謝過程(metabolic process)最多，分別有48個基因和41個基因。此外，還檢測到一些涉及先天免疫的生物過程，如免疫系統過程、行為、對刺激的反應等。

2.4 差異表達基因 COG 數據庫注釋分析

將差異表達基因比對到 COG 數據庫當中，預測基因的功能并進行分類和統計，可分為16類(圖3)。除了未知功能，所占數量最多的功能是蛋白質翻譯后修飾與轉運，分子伴侶(posttranslational modification,protein turnover,chaperones)，注釋到此類的差異基因有11個。其次是信號轉導機制(signal transduction mechanisms)、碳水化合物運輸與代謝(carbohydrate transport and metabolism)和細胞內運輸、分泌和囊泡運輸(intracellular trafficking,secretion,and vesicular transport)。

圖2 差異基因GO功能分析Fig.2 GO functional analysis of differential gene

2.5 差異基因 KEGG 功能分析

將差異表達基因進行KEGG生物通路富集分析，為挖掘細菌GD2增強松材線蟲寄主適應性相關基因提供依據。全部差異表達基因共映射到106個代謝通路，圖4顯示了富集程度在前50的通路，其中新陳代謝(metabolism)、有機系統(organismal systems)、人類疾病(human diseases)涉及的通路最多。ECM-receptor互作(ECM-receptor interaction)、細胞凋亡(cellular senescence)、腎素-血管緊張素系統(renin-angiotensin system)和PPAR信號通路(PPAR signaling pathway)、AMPK 信號通路(AMPK signaling pathway)、谷胱甘肽代謝(glutathione metabolism)等通路主要涉及機體的免疫調節、延長壽命和生長代謝。

B:染色質結構和變化 Chromatin structure and dynamics； C:能量生產和轉化 Energy production and conversion； E:氨基酸運輸和代謝 Amino acid transport and metabolism； F:核苷酸運輸和代謝 Nucleotide transport and metabolism； G:碳水化合物運輸和代謝 Carbohydrate transport and metabolism； I:脂質運輸和代謝 Lipid transport and metabolism； K:轉錄 Transcription； M:細胞壁/膜/包膜生物發生 Cell wall/membrane/envelope biogenesis； N:細胞運動 Cell motility； O:蛋白質翻譯后修飾與轉運,分子伴侶 Posttranslational modification,protein turnover,chaperones； P:無機離子轉運與代謝 Inorganic ion transport and metabolism； Q:次級代謝物的生物合成、運輸和代謝 Secondary metabolites biosynthesis,transport and catabolism； S:功能未知 Function unknown； T:信號轉導機制 Signal transduction mechanisms； U:細胞內運輸、分泌和囊泡運輸 Intracellular trafficking,secretion,and vesicular transport； Z:細胞骨架 Cytoskeleton.

圖4 差異基因 KEGG 功能分類Fig.4 KEGG functional classification of differential gene

2.6 qPCR驗證

通過RT-qPCR對8個表達豐度不同的轉錄本(5個上調和3個下調)進行驗證，發現基因的qPCR 驗證結果與轉錄組中的趨勢一致，說明轉錄組測序數據較為準確，可用于下一步研究(圖5)。

A.轉錄組所測8 個差異表達基因 Read count 值； B. RT-qPCR 檢測所得8 個差異表達基因的相對表達量。*P<0.05，**P<0.01，***P<0.001，****P<0.000 1。 A. Read count of the 8 differentially expressed genes measured in the transcriptome;B. The relative expression levels of the 8 differentially expressed genes detected by RT-qPCR. *P<0.05,**P<0.01,*** P<0.001,**** P<0.000 1.

3 討 論

RNA測序(RNA-Seq)是一種用高通量測序技術進行測序分析的轉錄組測序技術，目前已成為所有物種全基因組分析的第一方法，而鑒定比較組間的差異表達基因是RNA-Seq 最典型的應用[24-25]。盡管RNA-Seq是一項正在積極開發的技術，但它提供的基因關鍵功能仍使之成為目前最為廣泛使用的測序技術[26]。本研究通過轉錄組測序分析了細菌GD2對松材線蟲入侵寄主馬尾松24 h后的影響，但是我們發現松材線蟲轉錄組序列與參考基因組的比對率比較低。原因可能是①沒有確保每棵馬尾松都準確接種到15 000條松材線蟲；②在松材線蟲接種到馬尾松24 h 后，沒有將松材線蟲用貝爾曼漏斗法漏下線蟲進行RNA 提取，而是直接將接種部位取下，這時候已經有部分線蟲已經擴散到松樹其他地方，再加上馬尾松和松材線蟲混樣提取RNA，所以導致提取的RNA 里馬尾松的RNA 占比很大，只有較少的為松材線蟲的 RNA；③由于混樣提取RNA，為了剔除馬尾松的同源序列，我們將測得的松材線蟲的reads 與馬尾松的Clean reads 進行了比對，去除了比對上的reads，將剩下的reads 進行過濾后得到的松材線蟲的Clean reads 再與參考基因組比對。所以可能是松材線蟲的Clean reads 的基因覆蓋度比較低，導致了松材線蟲轉錄組序列與參考基因組的比對率比較低。

通過本次轉錄測序，從中篩選出143個差異表達顯著的基因，并對這些差異基因進行 GO、COG以及 KEGG 分析。

從GO功能注釋可以看出，差異表達基因主要集中在代謝過程、結合、催化活性、膜等與藥物的吸收、運輸、代謝有關的功能，這些基因可能參與了外源有害物質(馬尾松抗病反應中產生的植保素)的解毒作用。此外，涉及的GO功能免疫和行為是動物保護自己免受傷害的2種機制[27]。這可能是松材線蟲應對寄主的防御反應所產生的。

COG功能注釋結果表明，細菌GD2與松材線蟲混合接種馬尾松后，GD2幫助線蟲增強信號轉導機制來感受寄主的防御反應，然后線蟲通過改變自身的能量運輸與代謝以及蛋白質翻譯后修飾與轉運等方面來加強對寄主的適應性。

對差異基因富集的KEGG 通路進行分析，發現差異基因集中于視黃醇新陳代謝、酪氨酸代謝、果糖和甘露糖的代謝、鞘脂類代謝等代謝通路。可見細菌GD2對松材線蟲的生長代謝產生了影響。此外，一些差異基因如BXY_0706100涉及的PI3K-Akt信號通路在與細菌互作及對抗逆境脅迫中發揮著重要作用[28]，涉及的胰島素信號通路調節固有免疫[29]；BXY_1558900涉及的ECM-receptor互作通路被認為是免疫相關通路[30]；BXY_1767700涉及的細胞色素P450對外源物質的代謝通路可催化殺蟲劑、植物次生代謝物等外源物質的代謝，使生物體適應其生存環境或產生抗逆性[31]。

本研究通過 RNA-Seq，初步揭示了接種馬尾松后，馬尾松內生細菌GD2對松材線蟲轉錄組的影響，發現菌株GD2影響著松材線蟲的143個基因，這些DEGs顯著富集于一些與代謝、免疫相關的GO、COG類別和KEGG途徑，這些與松材線蟲的生命活動密切相關。推測菌株GD2在松材線蟲抵御寄主的防御反應時增強了線蟲的信號傳導機制，刺激了線蟲的免疫相關反應，提高了線蟲的能量運輸和對寄主所產生的次生代謝物代謝的能力，從而加重了松材線蟲病的發生。這為進一步從分子層面弄清細菌增強松材線蟲對寄主的適應性機制奠定了基礎。