陳陽(yáng)雪,趙曉佳,談家金
南京林業(yè)大學(xué)林學(xué)院,南京 210018
松材線蟲(chóng)是一種危險(xiǎn)性極高的外來(lái)有害生物,由它引起的松材線蟲(chóng)病是一種毀滅性的森林病害[1]。松材線蟲(chóng)病可使易感病松樹(shù)在感染松材線蟲(chóng)后2~3個(gè)月內(nèi)發(fā)生枯萎死亡[2]。自1905年在日本大暴發(fā)以來(lái),該病害在世界各地被發(fā)現(xiàn),截至目前該病已在歐亞和北美等多個(gè)國(guó)家或地區(qū)被發(fā)現(xiàn)和報(bào)道[3-4]。亞洲是松材線蟲(chóng)病流行危害最嚴(yán)重的地區(qū),其中中國(guó)是目前遭受松材線蟲(chóng)病威脅最嚴(yán)重的國(guó)家,病害已經(jīng)在中國(guó)10多個(gè)省造成了大面積松樹(shù)死亡,現(xiàn)有6 000萬(wàn)hm2松林正面臨著松材線蟲(chóng)病大流行的威脅[5]。經(jīng)統(tǒng)計(jì),該病害目前已在我國(guó)危害馬尾松(Pimusmassoniana)、油松(P.tabulaeformis)、赤松(P.densiflora)、華山松(P.armandii)、黑松(P.thunbergii)等50多種松科松屬樹(shù)種和落葉松(Larixolgensis)、云杉(Piceaasperata)、冷杉(Abiesfabri)等10多種非松屬樹(shù)種[6],對(duì)森林生態(tài)系統(tǒng)造成嚴(yán)重的破壞,也造成了嚴(yán)重的經(jīng)濟(jì)損失。
松材線蟲(chóng)病由松材線蟲(chóng)、寄主、傳播媒介、細(xì)菌真菌和環(huán)境等多個(gè)因素構(gòu)成了其復(fù)雜的病害系統(tǒng),因此松材線蟲(chóng)的致病機(jī)制尚未完全弄清[7]。目前,國(guó)內(nèi)外學(xué)者對(duì)松材線蟲(chóng)病的研究工作主要以松材線蟲(chóng)為病原展開(kāi)[1,8-10]。近年來(lái),已有研究顯示松材線蟲(chóng)致病作用是松材線蟲(chóng)與其共生細(xì)菌共同致病[11]。Kawazu等[12-15]、洪英娣等[16-17]從松材線蟲(chóng)體表分離出細(xì)菌,并通過(guò)接種試驗(yàn)表明松材線蟲(chóng)與其伴生細(xì)菌共同導(dǎo)致松材線蟲(chóng)病。但是細(xì)菌在這種復(fù)雜病害中的特定參與特性仍然存在爭(zhēng)議。Wang等[18]、He等[19]、Xue等[20]通過(guò)轉(zhuǎn)錄組測(cè)序研究了細(xì)菌對(duì)松材線蟲(chóng)基因表達(dá)的影響,發(fā)現(xiàn)細(xì)菌對(duì)松材線蟲(chóng)的致病相關(guān)基因有一定的影響。這些對(duì)于深入研究細(xì)菌在松材線蟲(chóng)病中的作用提供了重要的實(shí)驗(yàn)依據(jù)。菌株GD2是課題組前期從馬尾松體內(nèi)分離的內(nèi)生細(xì)菌,初步研究顯示其在松材線蟲(chóng)致病過(guò)程中發(fā)揮著重要作用[21]。本研究擬進(jìn)一步探討在松材線蟲(chóng)侵入寄主建立適應(yīng)性的過(guò)程中,菌株GD2對(duì)線蟲(chóng)轉(zhuǎn)錄組的影響,這對(duì)全面認(rèn)識(shí)松材線蟲(chóng)的致病機(jī)理和松材線蟲(chóng)-細(xì)菌微生態(tài)系統(tǒng)具有重要的理論意義。
供試菌株為馬尾松內(nèi)生細(xì)菌GD2。供試松材線蟲(chóng)AA3,分離自安徽的感病黃山松(P.taiwanensis)。蟲(chóng)株由南京林業(yè)大學(xué)森林病理實(shí)驗(yàn)室提供。供試松苗為生長(zhǎng)較一致的 2 年生馬尾松(中國(guó),廣西),由江蘇沂北園林苗木有限公司提供。
將松材線蟲(chóng) AA3 接種到長(zhǎng)滿新鮮灰葡萄孢的 PDA 平板上,置于 25 ℃ 培養(yǎng)6 d 后,采用貝爾曼漏斗法分離。將線蟲(chóng)液收集到無(wú)菌的 15 mL 離心管中,3 000 r/min 離心 3 min,棄上清,用 0.1%硫酸鏈霉素處理 15 min,用無(wú)菌水沖洗離心3次。將松材線蟲(chóng)調(diào)成濃度為 15 000 條/mL 的蟲(chóng)懸液。
用接種環(huán)挑取1環(huán)菌株GD2的新鮮單菌落,接入含有20 mL NB培養(yǎng)液的錐形瓶中,27 ℃、200 r/min振蕩搖培12 h,置入無(wú)菌的 50 mL離心管中,6 500 r/min 離心 3 min,棄上清,無(wú)菌去離子水洗 3 次。將菌體用無(wú)菌水重懸制成菌懸液(濃度為5×107CFU/mL)。
線蟲(chóng) AA3 與菌株 GD2 預(yù)處理后,采用皮接法[22]接種到2年生馬尾松。設(shè)3個(gè)處理,1個(gè)對(duì)照,每處理接種13棵馬尾松:(1)AA3+GD2(T):接種 1 mL 蟲(chóng)懸液和 1 mL 菌懸液;(2)AA3(CK):接種 1 mL 蟲(chóng)懸液和 1 mL無(wú)菌去離子水;(3)菌株 GD2:接種1 mL 菌懸液和 1 mL無(wú)菌去離子水;(4)AG:接無(wú)菌去離子水共 2 mL。放于溫室 25 ℃培養(yǎng),24 h后將接種無(wú)菌水(AG)、AA3(CK)、AA3+GD2(T)的馬尾松接種點(diǎn)部位剪下,-80 ℃保存?zhèn)溆?。每個(gè)處理各3個(gè)生物學(xué)重復(fù)。剩余的馬尾松每 7 d 觀察接種寄主的感病癥狀,統(tǒng)計(jì)感病率和感病指數(shù)。根據(jù)病害分級(jí)標(biāo)準(zhǔn)和參照Yu等[23]方法計(jì)算感病指數(shù)。馬尾松的感病嚴(yán)重程度分為5個(gè)級(jí)別,Ⅰ級(jí):植物針葉綠色,健康生長(zhǎng),代表值0;Ⅱ級(jí):0~25%針葉褪綠變黃,代表值1;Ⅲ級(jí):25%~50% 針葉褪綠變黃,枝頭彎曲,代表值2;Ⅳ級(jí):50%~70% 針葉褪綠變黃,枝頭下垂,代表值3;Ⅴ級(jí):75%~100% 針葉失綠褐變,植株枯萎死亡,代表值4。計(jì)算公式如下:
感病率=(感病株數(shù)/接種株樹(shù)) × 100% 感病指數(shù)= [∑(各級(jí)病樹(shù)株數(shù) × 代表值)/ (總株數(shù) × 發(fā)病最高級(jí)代表值)] × 100
將前面-80 ℃保存?zhèn)溆玫臉悠凡捎肙minPlant RNA kit(DNase I)(CWBIO)法提取各處理組的總RNA,用 1% 的瓊脂糖凝膠電泳檢測(cè),2100 Bioanalyser (Agilent) 、ND-2000 (NanoDrop Technologies)方法檢測(cè)RNA樣品的質(zhì)量。RNA文庫(kù)的建立采用TruSeqTM RNA sample preparation Kit( Illumina,San Diego,CA)試劑盒。使用IlluminaNovaSeq 6000測(cè)序平臺(tái)進(jìn)行高通量測(cè)序(上海美吉生物醫(yī)藥科技有限公司)。數(shù)據(jù)已提交NCBI,登錄號(hào)為PRJNA690180。
首先將測(cè)得馬尾松(AG組)的Raw reads先進(jìn)行過(guò)濾處理得到 Clean reads,再用CK組和T組的Raw reads直接比對(duì)馬尾松的Clean reads,去除比對(duì)上的reads,將剩下的reads經(jīng)過(guò)濾后比對(duì)到參考基因組松材線蟲(chóng)基因組(參考基因組版本:GCA_000231135.1;參考基因組來(lái)源:https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)進(jìn)行序列比對(duì)。
使用FPKM(Fragments Per Kilobases per Millionreads)軟件進(jìn)行定量計(jì)算獲得基因的表達(dá)量,再用軟件DESeq2進(jìn)行差異基因分析,并用BH(fdr correction with Benjamini/Hochberg)進(jìn)行多重檢驗(yàn)。將P<0.05 且 |log2FC|≥0.585設(shè)為條件篩選差異表達(dá)基因并進(jìn)行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/)分析。
選取8個(gè)差異表達(dá)基因(5個(gè)上調(diào)基因、3個(gè)下調(diào)基因),利用SIGMA在線設(shè)計(jì)軟件(http://www.oligoarchitect.com/LoginServlet)為 8個(gè)差異表達(dá)基因設(shè)計(jì)引物(表1),以松材線蟲(chóng)肌動(dòng)蛋白基因Actin(GenBank EU100952)為內(nèi)參基因,每個(gè)樣品3個(gè)技術(shù)學(xué)重復(fù),3個(gè)生物學(xué)重復(fù),并采用 2-ΔΔCt算法分析數(shù)據(jù)。
溫室接種試驗(yàn)表明,菌株GD2增強(qiáng)了松材線蟲(chóng)在馬尾松中的致病力。AA3 處理組、GD2+AA3混合處理組都使 2 年生馬尾松盆栽苗出現(xiàn)感病癥狀,并且所致萎蔫癥狀一致,均為針葉先褪色、黃化,進(jìn)而枯萎變成紅褐色;而 GD2 處理組和對(duì)照無(wú)菌水處理組均保持健康,馬尾松松針綠色。接種后21 d,接種 AA3 的馬尾松感病率為20%,感病指數(shù)為5,而混合接種AA3和菌株GD2的馬尾松感病率達(dá)到40%,感病指數(shù)為10;接種后 42 d, AA3處理組的感病率為 60%,感病指數(shù) 52.5,而GD2+AA3 處理組的感病率為 70%,感病指數(shù)為55(表2)。
表1 松材線蟲(chóng) 8 個(gè)差異表達(dá)基因的 qRT-PCR 引物 Table 1 The qRT-PCR primers for 8 differentially expressed genes of Bursaphelenchus xylophilus
表2 接種不同處理的松材線蟲(chóng)的馬尾松的感病率和感病指數(shù) Table 2 Disease incidence and disease index (DI) of Pinus massoniana inoculated with different treatments of Bursaphelenchus xylophilus
對(duì)含松材線蟲(chóng)的RNA樣本進(jìn)行轉(zhuǎn)錄組測(cè)序,去除比對(duì)到的馬尾松reads,剩下的 reads經(jīng)過(guò)濾后得到的clean reads質(zhì)量達(dá)到 Q20 以上的堿基比例大于 80%,達(dá)到 Q30 以上的堿基比例大于85%,每個(gè)樣本堿基 GC 含量較為一致,且測(cè)序堿基平均錯(cuò)誤率均在0.1%以下,說(shuō)明測(cè)序數(shù)據(jù)良好(表 3)。將過(guò)濾后的數(shù)據(jù)與rRNA數(shù)據(jù)庫(kù)進(jìn)行比對(duì)(表4),再利用HISAT2 軟件將質(zhì)控后的原始數(shù)據(jù),即Clean reads,與參考基因組(https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)比對(duì)(表5)。將P<0.05 且|log2FC|≥0.585設(shè)為條件共篩選到143個(gè)差異表達(dá)基因,其中 63 個(gè)基因表達(dá)上調(diào),80個(gè)基因表達(dá)下調(diào)(圖1)。這些差異基因與內(nèi)生細(xì)菌GD2協(xié)同松材線蟲(chóng)侵染馬尾松相關(guān)。
表3 轉(zhuǎn)錄組序列質(zhì)量統(tǒng)計(jì)表 Table 3 The data statistics for RNA-sequencing
表4 過(guò)濾后數(shù)據(jù)與rRNA數(shù)據(jù)庫(kù)比對(duì)統(tǒng)計(jì)表 Table 4 Results of comparison with rRNA database
表5 序列與松材線蟲(chóng)參考基因組比對(duì)結(jié)果統(tǒng)計(jì) Table 5 Results of comparison with the reference genomeof Bursaphelenchus xylophilus
圖1 差異基因火山圖Fig.1 Differential gene volcano
對(duì)差異表達(dá)基因進(jìn)行 GO 功能注釋與分類統(tǒng)計(jì)。差異基因被歸類到GO 3個(gè)組分的38個(gè)功能分類中(圖2)。細(xì)胞組成(cellular component)中細(xì)胞組分(cell part)和膜組分(membrane part)最多,分別有52個(gè)基因和44個(gè)基因;分子功能(molecular functions)中結(jié)合(binding)和催化活性(catalytic activity)最多,分別有60個(gè)基因和51個(gè)基因,生物代謝過(guò)程(biological processes)中細(xì)胞進(jìn)程(cellular process)和代謝過(guò)程(metabolic process)最多,分別有48個(gè)基因和41個(gè)基因。此外,還檢測(cè)到一些涉及先天免疫的生物過(guò)程,如免疫系統(tǒng)過(guò)程、行為、對(duì)刺激的反應(yīng)等。
將差異表達(dá)基因比對(duì)到 COG 數(shù)據(jù)庫(kù)當(dāng)中,預(yù)測(cè)基因的功能并進(jìn)行分類和統(tǒng)計(jì),可分為16類(圖3)。除了未知功能,所占數(shù)量最多的功能是蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn),分子伴侶(posttranslational modification,protein turnover,chaperones),注釋到此類的差異基因有11個(gè)。其次是信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(signal transduction mechanisms)、碳水化合物運(yùn)輸與代謝(carbohydrate transport and metabolism)和細(xì)胞內(nèi)運(yùn)輸、分泌和囊泡運(yùn)輸(intracellular trafficking,secretion,and vesicular transport)。
圖2 差異基因GO功能分析Fig.2 GO functional analysis of differential gene
將差異表達(dá)基因進(jìn)行KEGG生物通路富集分析,為挖掘細(xì)菌GD2增強(qiáng)松材線蟲(chóng)寄主適應(yīng)性相關(guān)基因提供依據(jù)。全部差異表達(dá)基因共映射到106個(gè)代謝通路,圖4顯示了富集程度在前50的通路,其中新陳代謝(metabolism)、有機(jī)系統(tǒng)(organismal systems)、人類疾病(human diseases)涉及的通路最多。ECM-receptor互作(ECM-receptor interaction)、細(xì)胞凋亡(cellular senescence)、腎素-血管緊張素系統(tǒng)(renin-angiotensin system)和PPAR信號(hào)通路(PPAR signaling pathway)、AMPK 信號(hào)通路(AMPK signaling pathway)、谷胱甘肽代謝(glutathione metabolism)等通路主要涉及機(jī)體的免疫調(diào)節(jié)、延長(zhǎng)壽命和生長(zhǎng)代謝。
B:染色質(zhì)結(jié)構(gòu)和變化 Chromatin structure and dynamics; C:能量生產(chǎn)和轉(zhuǎn)化 Energy production and conversion; E:氨基酸運(yùn)輸和代謝 Amino acid transport and metabolism; F:核苷酸運(yùn)輸和代謝 Nucleotide transport and metabolism; G:碳水化合物運(yùn)輸和代謝 Carbohydrate transport and metabolism; I:脂質(zhì)運(yùn)輸和代謝 Lipid transport and metabolism; K:轉(zhuǎn)錄 Transcription; M:細(xì)胞壁/膜/包膜生物發(fā)生 Cell wall/membrane/envelope biogenesis; N:細(xì)胞運(yùn)動(dòng) Cell motility; O:蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn),分子伴侶 Posttranslational modification,protein turnover,chaperones; P:無(wú)機(jī)離子轉(zhuǎn)運(yùn)與代謝 Inorganic ion transport and metabolism; Q:次級(jí)代謝物的生物合成、運(yùn)輸和代謝 Secondary metabolites biosynthesis,transport and catabolism; S:功能未知 Function unknown; T:信號(hào)轉(zhuǎn)導(dǎo)機(jī)制 Signal transduction mechanisms; U:細(xì)胞內(nèi)運(yùn)輸、分泌和囊泡運(yùn)輸 Intracellular trafficking,secretion,and vesicular transport; Z:細(xì)胞骨架 Cytoskeleton.
圖4 差異基因 KEGG 功能分類Fig.4 KEGG functional classification of differential gene
通過(guò)RT-qPCR對(duì)8個(gè)表達(dá)豐度不同的轉(zhuǎn)錄本(5個(gè)上調(diào)和3個(gè)下調(diào))進(jìn)行驗(yàn)證,發(fā)現(xiàn)基因的qPCR 驗(yàn)證結(jié)果與轉(zhuǎn)錄組中的趨勢(shì)一致,說(shuō)明轉(zhuǎn)錄組測(cè)序數(shù)據(jù)較為準(zhǔn)確,可用于下一步研究(圖5)。
A.轉(zhuǎn)錄組所測(cè)8 個(gè)差異表達(dá)基因 Read count 值; B. RT-qPCR 檢測(cè)所得8 個(gè)差異表達(dá)基因的相對(duì)表達(dá)量。*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.
RNA測(cè)序(RNA-Seq)是一種用高通量測(cè)序技術(shù)進(jìn)行測(cè)序分析的轉(zhuǎn)錄組測(cè)序技術(shù),目前已成為所有物種全基因組分析的第一方法,而鑒定比較組間的差異表達(dá)基因是RNA-Seq 最典型的應(yīng)用[24-25]。盡管RNA-Seq是一項(xiàng)正在積極開(kāi)發(fā)的技術(shù),但它提供的基因關(guān)鍵功能仍使之成為目前最為廣泛使用的測(cè)序技術(shù)[26]。本研究通過(guò)轉(zhuǎn)錄組測(cè)序分析了細(xì)菌GD2對(duì)松材線蟲(chóng)入侵寄主馬尾松24 h后的影響,但是我們發(fā)現(xiàn)松材線蟲(chóng)轉(zhuǎn)錄組序列與參考基因組的比對(duì)率比較低。原因可能是①?zèng)]有確保每棵馬尾松都準(zhǔn)確接種到15 000條松材線蟲(chóng);②在松材線蟲(chóng)接種到馬尾松24 h 后,沒(méi)有將松材線蟲(chóng)用貝爾曼漏斗法漏下線蟲(chóng)進(jìn)行RNA 提取,而是直接將接種部位取下,這時(shí)候已經(jīng)有部分線蟲(chóng)已經(jīng)擴(kuò)散到松樹(shù)其他地方,再加上馬尾松和松材線蟲(chóng)混樣提取RNA,所以導(dǎo)致提取的RNA 里馬尾松的RNA 占比很大,只有較少的為松材線蟲(chóng)的 RNA;③由于混樣提取RNA,為了剔除馬尾松的同源序列,我們將測(cè)得的松材線蟲(chóng)的reads 與馬尾松的Clean reads 進(jìn)行了比對(duì),去除了比對(duì)上的reads,將剩下的reads 進(jìn)行過(guò)濾后得到的松材線蟲(chóng)的Clean reads 再與參考基因組比對(duì)。所以可能是松材線蟲(chóng)的Clean reads 的基因覆蓋度比較低,導(dǎo)致了松材線蟲(chóng)轉(zhuǎn)錄組序列與參考基因組的比對(duì)率比較低。
通過(guò)本次轉(zhuǎn)錄測(cè)序,從中篩選出143個(gè)差異表達(dá)顯著的基因,并對(duì)這些差異基因進(jìn)行 GO、COG以及 KEGG 分析。
從GO功能注釋可以看出,差異表達(dá)基因主要集中在代謝過(guò)程、結(jié)合、催化活性、膜等與藥物的吸收、運(yùn)輸、代謝有關(guān)的功能,這些基因可能參與了外源有害物質(zhì)(馬尾松抗病反應(yīng)中產(chǎn)生的植保素)的解毒作用。此外,涉及的GO功能免疫和行為是動(dòng)物保護(hù)自己免受傷害的2種機(jī)制[27]。這可能是松材線蟲(chóng)應(yīng)對(duì)寄主的防御反應(yīng)所產(chǎn)生的。
COG功能注釋結(jié)果表明,細(xì)菌GD2與松材線蟲(chóng)混合接種馬尾松后,GD2幫助線蟲(chóng)增強(qiáng)信號(hào)轉(zhuǎn)導(dǎo)機(jī)制來(lái)感受寄主的防御反應(yīng),然后線蟲(chóng)通過(guò)改變自身的能量運(yùn)輸與代謝以及蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)等方面來(lái)加強(qiáng)對(duì)寄主的適應(yīng)性。
對(duì)差異基因富集的KEGG 通路進(jìn)行分析,發(fā)現(xiàn)差異基因集中于視黃醇新陳代謝、酪氨酸代謝、果糖和甘露糖的代謝、鞘脂類代謝等代謝通路。可見(jiàn)細(xì)菌GD2對(duì)松材線蟲(chóng)的生長(zhǎng)代謝產(chǎn)生了影響。此外,一些差異基因如BXY_0706100涉及的PI3K-Akt信號(hào)通路在與細(xì)菌互作及對(duì)抗逆境脅迫中發(fā)揮著重要作用[28],涉及的胰島素信號(hào)通路調(diào)節(jié)固有免疫[29];BXY_1558900涉及的ECM-receptor互作通路被認(rèn)為是免疫相關(guān)通路[30];BXY_1767700涉及的細(xì)胞色素P450對(duì)外源物質(zhì)的代謝通路可催化殺蟲(chóng)劑、植物次生代謝物等外源物質(zhì)的代謝,使生物體適應(yīng)其生存環(huán)境或產(chǎn)生抗逆性[31]。
本研究通過(guò) RNA-Seq,初步揭示了接種馬尾松后,馬尾松內(nèi)生細(xì)菌GD2對(duì)松材線蟲(chóng)轉(zhuǎn)錄組的影響,發(fā)現(xiàn)菌株GD2影響著松材線蟲(chóng)的143個(gè)基因,這些DEGs顯著富集于一些與代謝、免疫相關(guān)的GO、COG類別和KEGG途徑,這些與松材線蟲(chóng)的生命活動(dòng)密切相關(guān)。推測(cè)菌株GD2在松材線蟲(chóng)抵御寄主的防御反應(yīng)時(shí)增強(qiáng)了線蟲(chóng)的信號(hào)傳導(dǎo)機(jī)制,刺激了線蟲(chóng)的免疫相關(guān)反應(yīng),提高了線蟲(chóng)的能量運(yùn)輸和對(duì)寄主所產(chǎn)生的次生代謝物代謝的能力,從而加重了松材線蟲(chóng)病的發(fā)生。這為進(jìn)一步從分子層面弄清細(xì)菌增強(qiáng)松材線蟲(chóng)對(duì)寄主的適應(yīng)性機(jī)制奠定了基礎(chǔ)。