油茶根系與內生細菌枯草芽孢桿菌互作早期的轉錄組分析*

李梓楊 陳曉琳 李麗麗 許詩萍 何苑皞

(1.南方人工林病蟲害防治國家林業和草原局重點實驗室 森林有害生物防控湖南省重點實驗室 經濟類培育與保護教育部重點實驗室 中南林業科技大學 長沙 410004)

內生細菌廣泛存在于健康植物體內，具有促進植物生長、提高作物產量、提高植物存活率、抑制病原菌生長、修復環境污染、溶解磷酸鹽、促進氮吸收等作用(Bulgarellietal., 2012),通過長期的進化與寄主植物形成了互利互惠和的關系。內生細菌在互作早期，利用自身的趨化性、運動性、黏附作用以及分泌一些小分子物質與植物建立互作關系實現定殖并發揮作用；植物通過激素應答、免疫反應啟動、次生代謝物合成、碳代謝等一系列的變化來響應內生細菌的定殖。植物根系或根系分泌物接觸后內生菌轉錄組變化明顯，其中氨基酸轉運、抗壓、IV型菌毛、粘附、植物來源化合物的利用等相關基因的表達變化顯著(Shidoreetal., 2014; Yietal., 2014; Pankieviczetal., 2014; Sheibani-Tezerjietal., 2015)。植物根系在接觸內生菌后，受體激酶、植物激素信號通路、木質素積累、細胞壁合成相關基因的表達變化明顯(Drogueetal., 2014; Irizarry, 2018; Galambosetal., 2020)。目前研究主要以互作后期為主，而互作早期的研究較少，而互作早期植物各項生理變化較為顯著，是互作的關鍵時間，以互作早期的植物根系為研究對象。目前，內生細菌與植物互作研究主要集中在擬南芥(Arabidopsisthaliana)、水稻(Oryzasativa)、棉花(Anemonevitifolia)等模式生物上，對于木本植物則關注較少，而內生細菌與油茶(Camelliaolefolia)互作的研究還未見報道。

筆者課題組前期從油茶內分離出1株內生細菌——枯草芽孢桿菌(Bacillussubtilis1-L-29gfpr)，該菌株可成功在油茶根系內定殖，并具有促進生長、抗病作用(Xuetal., 2020)。在本研究將該菌株接種油茶根系構建油茶-內生細菌活體互作體系，利用轉錄組測序技術研究不同時間(0、6、12、24 h)油茶根系轉錄組的變化，旨在解析油茶根系響應內生細菌侵染的變化規律，為后期深入研究內生細菌與油茶的互作機制提供參考。

1 材料與方法

1.1 試驗材料

將油茶種子表面消毒后，放入無菌組培瓶中，17 ℃保濕培養，直至幼莖5 cm左右。內生細菌B.subtilis1-L-29gfpr為中南林業科技大學微生物菌種保藏室提供。

1.2 試驗方法

1.2.1 油茶接種內生菌 將培養18 h的內生細菌1-L-29gfpr離心后重懸于PBS緩沖液中，濃度調至1×108cfu·mL-1，用注射器吸取20 mL緩沖液注射到組培瓶中，將油茶根系浸泡于緩沖液中， 20 min后將根系取出以無菌水沖洗3遍后置于無菌組培瓶中繼續保濕培養。分別于0 h(不接種對照)、接種后6 h(6 hpi)、12 h(12 hpi)、24 h(24 hpi)截取油茶根系根尖部位。取樣后將樣品液氮速凍，于-80 ℃保存。每個處理3個生物學重復。

1.2.2 RNA提取、文庫構建及測序 采用TRIZOL試劑盒(Invitrogen，美國) 分別提取0 、6 、12 、24 h共 12 個樣品的總RNA，按照試劑盒說明書進行RNA提取。利用 Agilent 2100 Bioanalyzer (Agilent，美國)檢測RNA質量，用帶有 Oligo(dT) 的磁珠富集真核生物Mrna； 使用緩沖液將富集的mRNA片段化為短片段，并使用隨機引物將其逆轉錄為cDNA； 通過DNA聚合酶I，RNase H，dNTP和緩沖液合成第二鏈cDNA； 用QiaQuick PCR提取試劑盒(Qiagen，荷蘭)純化cDNA片段，黏性末端修復，添加poly(A)，并連接測序接頭。通過瓊脂糖凝膠電泳選擇連接產物的大小，PCR擴增，委托廣州基迪奧生物科技有限公司使用Illumina Novaseq 6000進行測序。

1.2.3 差異表達基因 qRT-PCR 驗證 測序獲得的原始序列去除低質量序列得到Clean reads。使用對比軟件bowtie2將測序數據與枯草芽孢桿菌ASM211716V1菌株的基因組進行mapping(Langmeadetal., 2012)，mapping后的枯草芽孢桿菌reads剔除掉，將所有油茶樣本剩余的reads進行混合組裝，得到油茶的參考序列信息。然后，使用序列聚類軟件 TGICL 做進一步序列拼接和去冗余處理(Perteaetal., 2003)，得到Unigene。FPKM (fragments per kilobases per millionreads)值用于評估每個基因的表達水平。差異表達基因的篩選條件為P-adjust<0.05 且∣log2FC∣>=1。采用 NR、Swiss-Prot、Pfam、EggNOG、GO 和 KEGG 等數據庫進行基因功能注釋。通過 GO 和 KEGG 富集分析，確定差異表達基因主要富集的通路及其可能的生物學功能。

1.2.4 差異表達基因 qRT-PCR 驗證 為了驗證RNA-Seq 數據的可靠性，隨機選取5個差異表達的基因進行 qRT-PCR表達量驗證。利用軟件Primer5.0設計引物，內參基因選用EF1α，引物列表見表1，隨后采用 ABI 7000 系統進行 PCR擴增和分析。每個基因設 3 個重復。采用2-ΔΔCt法計算不同侵染時間點差異基因相對表達量。

表1 差異表達基因 qRT-PCR 驗證使用的引物Tab.1 Primers used for qRT-PCR gene expression analysis

2 結果與分析

2.1 測序數據分析

油茶和內生細菌互作的文庫0 (不接種對照)、6 、12 、24 h共產生序列52 958 922條，約46.1 Gb數據，差異表達基因12 561條，GC含量變化為43.14%～44.47%。過濾后Q30比例為91.94%以上。組裝的Unigene序列N50為1 175，Unigene平均長度為866 bp，數據的比對率為71.13%。通過blastx將Unigene序列比對到蛋白數據庫nr、SwissProt、KEGG 和 COG/KOG(evalue<0.000 01)發現有67 700獲得了注釋(44.5%)(表2)。對樣品的生物學重復和樣品組間的差異表達情況進行評估，相關系數均大于0.8，這表明重復樣品間的相關性較強。

表2 內生細菌接種油茶不同時間點油茶轉錄組測序數據Tab.2 Summary of De novo data generated for C. oleifera transcriptomes inoculated with endophytic bacteria at different time points

2.2 差異表達基因分析

通過對轉錄組測序后的基因進行注釋和表達量(RPKM 值)計算可知，3個時間點共有相同差異表達基因84個，內生菌接種油茶根部6 h時的差異表達基因數量最多為5 715個，其次為24 h，差異表達基因數量為1 871個(圖1)。

圖1 3個不同侵染時間點的油茶差異表達基因分析Fig. 1 DEGs of three different infection time points of C. oleiferaA： 不同時間點基因表達差異； B： 特異差異基因和共有差異基因的維恩圖分析。 A： Differences in gene expression at different time points； B： Venn diagrams representing the numbers of DEGs and the overlaps of sets obtained across four comparisons.

2.3 差異表達基因的 GO 富集分析

以 correctedP-value ≤0.05 為篩選閾值，對不同時間獲得的差異表達基因進行 GO 注釋分析，這些差異表達基因被富集注釋到3個 GO 分類，分別是分子功能(molecular function)、細胞組成(cellular component)、生物過程 (biological process)。差異表達基因生物學過程中有2 188 個差異表達基因；細胞組分有3 489 個； 分子功能有2 490 個。

由圖2可知，油茶根系接種內生細菌6 h時生物學過程中磷酸代謝途徑、含磷化合物代謝途徑、響應激素途徑、茉莉酸代謝途徑、響應細菌途徑等得到富集，這些途徑均與植物響應細菌刺激有關； 分子功能中的萜烯合酶活性，雙加氧酶活性等得到富集； 細胞組分中膜、膜固有成分等得到富集。接種12 h時,富集的生物學過程主要有單個有機體代謝過程、小分子代謝過程、羧酸代謝過程、酮酸代謝過程、有機酸代謝過程； 富集的分子功能主要有氧化還原酶活性、四吡咯綁定等； 富集的細胞組成有膜結合細胞器、胞內膜細胞器等。接種24 h時，代謝途徑、單有機體代謝途徑、小分子代謝途徑等生物學過程得到富集； 富集的分子功能主要有催化活性； 細胞組分中膜結合細胞器、細胞等得到富集(圖2)。

由GO富集分析中可知，油茶根系接種枯草芽孢桿菌1-L-29gfpr6 h后，生物學過程分子功能變化最明顯。

圖2 油茶-內生細菌互作早期差異表達基因GO富集分析Fig. 2 GO enrichment analysis of DEGs at early stage of C. olefolia-endophytic bacteria

2.4 差異表達基因的 KEGG 代謝通路富集分析

對不同時間點的差異表達基因進行KEGG pathways富集分析表明： 差異表達基因被富集到碳水化合物代謝、氨基酸代謝、信號轉導、環境適應性等通路。

油茶根系接種內生細菌1-L-29gfpr6 h，植物-病原互作通路、植物激素信號轉導、MAPK信號途徑、芪類化合物、二芳基庚烷類化合物及姜辣素生物合成通路都有一定程度的提高，這些途徑均參與植物對病原菌的響應。接種12 h，玉米素生物合成、核黃素代謝、雙萜生物合成有一定程度的提高。接種24 h，酪氨酸代謝、硫代謝等有一定程度的提高。這些代謝途徑可能參與植物與內生菌的互作(圖3)。由此可見，互作過程中6 h 富集的代謝通路與12 h 和24 h的差別較大。

圖3 油茶-內生細菌互作早期差異表達基因的 KEGG 富集分析Fig. 3 KEGG enrichment analysis of DEGs at early stage of C. olefolia-endophytic bacteria

2.5 植物激素信號轉導相關的差異基因

調節植物激素水平是內生細菌重要的特點之一。在內生細菌與宿主植物互作過程中生長素、赤霉素、脫落酸、茉莉酸、水楊酸等植物激素信號轉導相關基因均發生變化(Pinskietal.， 2019)。本研究中植物激素轉導途徑差異表達基因占總差異表達基因的比例隨時間遞減，6 、12 、24 h分別為6.41%、 3.96%、1.98%。由圖4A和表3可知，生長素響應基因(SAUR15/32/50)、生長素早期響應基因(IAA)表達量在6 h或12 h都有一定程度的提高，表明內生細菌1-L-29gfpr在一定程度可以刺激油茶根系的生長。木葡聚糖內糖基轉移酶/水解蛋白酶基因(XTH22/23/25)、TIFY家族基因(TIFY9/10B)、轉錄因子AP2/ERF(ERF1B)表達量在6 h顯著上調隨后下降。轉錄因子PIF(PIF3)、轉錄因子bZIP、DELLA蛋白編碼基因(GAI)表達量呈現低-高-低的狀態，在12 h最大。生長素響應蛋白編碼基因(GH3.1)表達量則在24 h達到最大值。

2.6 植物病原互作相關差異基因

油茶根部接種內生細菌1-L-29gfpr后，植物病原互作通路檢測到差異基因100個，其差異基因占總差異基因比例6、12、24 h分別為9.61%、3.96%、3.06%。油茶根部多個參與植物病原互作的基因均有上調，但到達峰值時間不同。由圖4B和表4可知，鈣調蛋白編碼基因(CAM2)、呼吸爆發氧化酶同源基因(RBOH)以及RPM1互作蛋白編碼基因(RIN4)表達量在6 h達到峰值，隨后一直下降。轉錄因子MYB、MAPK激酶(NTF4)表達量呈現高-低-高的狀態，在12 h達到最低值。抗病蛋白相關基因(RPM1、RPS5)表達量則在12 h達到最大值。由此可見，內生菌接種12 h后植物啟動抗病蛋白相關基因，LysM類受體激酶基因(CERK1)則持續性的下調。

2.7 苯丙素生物合成相關差異基因

苯丙素參與植物對生物和非生物刺激的的反應(Vogt， 2010)。接種內生細菌后6 、12 、24 h，苯丙素生物合成途徑差異基因占總差異基因比例分別為5.74%、3.96%、2.52%。在接種內生細菌后的不同時間點，檢測到差異基因51個，6 h上調基因35個，下調基因16個； 12 h上調基因33個，下調基因18個； 24 h上調基因21個，下調基因30個。由圖4C和表5可知，過氧化物酶基因(pod，PER5/44/58，PNC2)表達量在6 h達到峰值。糖苷水解酶(BGLU12、AA7GT、BGLU11)表達量在12 h達到峰值。肉桂醇脫氫酶(CAD1)表達量在24 h達到峰值。其中，大部分的基因都參與羥基苯基木質素、愈創木基木質素、紫丁香木基木質素以及香豆素的合成。

2.8 差異表達基因 qRT-PCR 驗證

qRT-PCR 結果表明，隨機挑取的 5個基因在不同時間點的差異表達量變化與轉錄組測序結果的趨勢一致(圖5)，說明本次轉錄組測序的結果可靠。

3 討論

3.1 植物激素信號轉導相關的差異基因

內生菌或病原菌的定殖都可引起植物體內激素的變化，植物響應內生菌或病原菌最為重要的是三類激素：水楊酸(SA)、茉莉酸(JA)和乙烯(ET)(Dominiketal., 2012)。本研究中植物激素轉導途徑差異基因占總差異基因比例從6.41%下降到1.98%，水楊酸、茉莉酸途徑產生了變化，但乙烯未發生明顯變化。在JA信號轉導過程中，ZIM(TIFY)結構域可介導JAZ蛋白與其共抑制因子NINJA相互作用，進而共同抑制JA信號轉導(段龍飛等， 2013)。本研究中JAZ編碼基因TIFY9、TIFY108在接種6 h表達量上升，隨后下降，這表明在內生菌與植物互作開始階段JA信號轉導受到抑制，12 h抑制解除。這與King等(2019)研究相似，接種內生菌后JA的響應會存在延時效應。TGA 基因家族是 bZIP 轉錄因子家族中十分重要的一組，它可與水楊酸信號途徑關鍵調節因子 NPR1 相互作用，促進 TGA 對下游PR1啟動子的結合并上調其表達，提高植株抗病性(Johnsonetal., 2003)。本研究中TGA表達量在接種6 h時下降，但24 h時PR1是有明顯的上調表達。由此可見，內生菌接種后，油茶根部的水楊酸(SA)、茉莉酸(JA)開始時受到抑制，但隨后抑制解除。這可能是由于內生菌為了定殖成功，在定殖初期(6 h)迫使植物降低了免疫反應，定殖成功后內生菌解除了免疫抑制，從而提高了植物的免疫力，這也證明內生菌可提高植物抗性。關于內生菌是如何迫使植物降低免疫反應還有待進一步研究。

圖4 差異基因表達熱圖Fig. 4 Heat map diagram of the expression levels of DEGsA： 植物激素信號轉導相關差異基因表達熱圖； B： 植物病原互作相關差異基因表達熱圖； C： 苯丙素生物合成相關差異基因表達熱圖。A:Heat map diagram of the expression levels of DEGs which are involved in plant hormone signal transduction. B:Heat map diagram of the expression levels of DEGs which are involved in plant-pathogen interaction. C:Heat map diagram of theexpression levels of DEGs which are involved in plant phenylpropanoid biosynthesis.

圖5 qRT-PCR表達量分析和轉錄組數據的驗證Fig. 5 Expression analysis of qRT-PCR and verification of transcriptome dataA： TCH4基因； B:RBOHC基因； C:POD基因； D:TIFY基因； E:XTH基因。 A： TCH4 gene； B:RBOHC gene； C:POD gene； D:TIFY gene； E:XTH gene.

表3 植物激素信號轉導通路差異基因Tab.3 DEGs of plant hormone signal transduction pathway

內生細菌通過調控植物激素信號傳遞，起到促進植物生長的作用(Straubetal., 2013)。赤霉素(GA)在植物株高等生長發育過程中起著重要的調控作用(張佳琦等， 2019)。GA 信號通路中的關鍵組分 DELLA 蛋白，主要起阻遏作用，GA可以解除DELLA 蛋白的阻遏作用。本研究中DELLA 蛋白基因表達量最后降低，表明內生菌定殖后期可以促進油茶生長。同時本研究也發現，DELLA 蛋白基因表達量與JA的相關基因表達量呈現相反的變化趨勢，這與Yimer等(2018)的研究結果一致，即GA和JA之間存在對立關系。推測其原因可能是植物將DELLA 蛋白作為調節生長與防御之間平衡的手段之一。

本研究中，接種內生細菌6 h后油茶根系生長素Aux/IAA基因表達量明顯上升，但隨后開始出現下降，這與Irizarry等(2018)研究結果相反，這可能是由于本研究選取的時間點為侵染前期，而Irizarry的研究選取的是接種10天后的樣品。

表4 植物病原互作通路差異基因Tab.4 DEGs of plant-pathogen interaction pathway

表5 苯丙素生物合成通路差異基因Tab.5 DEGs of phenylpropanoid biosynthesis pathway

3.2 植物病原互作相關差異基因

植物與微生物共同進化過程中形成了復雜的免疫防衛體系。該免疫防衛體系的正常運行保證了植物的正常生長。由于植物細胞可以通過先天免疫反應對MAMPs進行監測及響應(Machoetal., 2014)，因此內生菌或病原菌的定殖都可引起植物的免疫反應。Drogue等(2014)研究發現,內生菌可以抑制更多與植物防御有關基因的表達。在本研究中，內生菌在6 h引起了植物病原互作通路大量基因的表達，但在12 h表達量明顯下降，表明在內生菌侵染的初期油茶起動了防御反應，到后期防御反應減弱，讓內生菌更好的實現定殖。

轉錄因子在調節植物與微生物互作相關基因的過程中起著重要的作用。已有研究證明WRKY轉錄因子參與大量的植物-病原菌互作反應(Banerjeeetal., 2015; Ishihamaetal., 2012; Pandeyetal., 2009; Ryuetal., 2006)，但是對于內生菌與植物相互作用中的WRKY轉錄因子知之甚少。本研究中油茶根系的WRKY22/24/27/29/33都在6 h表達量上調，隨后出現一定的下調。Jorge對有益真菌與擬南芥互作過程中WRKY基因的表達進行了研究，發現不同的轉錄因子調節的途徑不同，其中WRKY33是受茉莉酸調節，并且在24 h的表達量是下降的(Saenz-Mataetal., 2014)。

鈣是植物中PAMP觸發的免疫力重要信號(Yuanetal., 2017)。植物環核苷酸門控通道(CNGC)蛋白介導的鈣離子內流是識別受體復合物與鈣離子依賴性免疫的一個重要關聯。Tian等(2019)研究表明，CNGC通道直接被BIK1磷酸化并響應flg22激活PTI信號傳導途徑。本研究中，內生細菌枯草芽孢桿菌也存在flg22，CNGC編碼基因表達量上調可能與此有關。同時本研究中油茶根系的鈣調蛋白基因(CAM2)在6 h上調明顯，主要是由于鈣調蛋白參與了植物免疫反應(Singhetal., 2020)。

3.3 苯丙素生物合成相關差異基因

植物的苯丙素化合物在抵抗病原體侵襲中起重要作用。苯丙素合成途徑中的木質素可以使細胞壁增厚，抵抗真菌侵染釘形成的壓力，同時木質化的細胞壁可以減弱病原菌細胞壁降解酶的作用(Bechingeretal., 1999; Naoumkinaetal., 2010)。苯丙素合成途徑在病原菌與植物互作過程中變化較大，但內生菌與植物互作中的研究報道較少。本研究中3個主要單元： 對羥基苯基(H)、愈創木脂(G)和丁香基(S)的基因表達量在接種6 h和12 h均有上調，而在24 h下調。王恒照等(2019)對促生菌B.pumilusLZP02接種后的水稻根部的研究發現，LZP02可提高苯丙素合成相關基因的表達。這一結果與本研究存在一定的差異，其主要原因是本研究為12 h與24 h的表達量比較，如果將24 h與未接種內生菌的試驗組比較，苯丙素合成基因的表達量仍為上升，同時Hong的取樣時間為促生菌接種48 h后，也與本研究存在一定不同。

豆香素是苯丙素合成途徑中另一類非常重要的植物次生代謝產物。香豆素是木質素合成的前體，并且香豆素本身就是一種植物抗毒素，可以防止病蟲的侵害，同時豆香素還對植物根部微生物組的裝配起著重要的作用(Stringlisetal., 2018)。本研究中香豆素合成基因表達量也出現了相應的變化，在12 h達到峰值，這表明油茶內生細菌可通過調節宿主信號傳導，從而提高宿主對病原菌的抗性。

4 結論

對內生細菌接種后不同時間點的油茶根部轉錄組進行測序分析，得到12 561個差異表達基因。差異表達基因在6 h最多，主要集中在植物-病原互作通路、植物激素信號轉導、苯丙素化合物等通路。其中htpG、鈣調蛋白基因(CAM2)、呼吸爆發氧化酶編碼基因(RBOH)、轉錄因子WRKY22、鈣調蛋白基因(CAM2)、木糖轉移酶基因(TCH4)、茉莉酸受體ZIM結構域包含蛋白編碼基因(TIFY10B)、Aux/IAA基因IAA1、過氧化物酶編碼基因等上調顯著。接種后6 h差異表達基因最多，為互作的關鍵時間。