許秀玉 張 勇 甘先華 張華新 張應(yīng)中 王明懷 毛亦楊

（1 廣東省森林培育與保護(hù)利用重點(diǎn)實驗室/廣東省林業(yè)科學(xué)研究院，廣東 廣州 510520；2 中國林業(yè)科學(xué)研究院 國家林業(yè)局鹽堿地研究中心， 北京 100091； 3 中國林業(yè)科學(xué)研究院 熱帶林業(yè)研究所， 廣東 廣州 510520； 4 廣東省森林資源保育中心， 廣東 廣州 510173）

植物青枯病（Ralstonia solanacearum）抗性性狀受多基因控制，屬于數(shù)量性狀，采用常規(guī)育種方法選育抗病品系難度大，進(jìn)展慢[1]。木麻黃科（Casuarinaceae）植物自然分布區(qū)廣，不同樹種、種源及無性系間遺傳變異非常豐富，在青枯病抗性上也存在著極顯著分化[2]。近年來，隨著越來越多木本植物[3-6]完成全基因組測序及大量分子標(biāo)記的開發(fā)，數(shù)量性狀的研究，特別是利用關(guān)聯(lián)分析挖掘數(shù)量性狀相關(guān)功能基因已經(jīng)越來越成為國際植物基因組學(xué)研究的熱點(diǎn)。林木繁衍以自然雜交為主，群體遺傳多樣性水平高，連鎖不平衡水平低，特別適用于關(guān)聯(lián)分析[7]。目前，許多林木樹種開展了數(shù)量性狀與分子標(biāo)記的關(guān)聯(lián)分析，但大多集中在生長性狀與木材性狀[8]，利用關(guān)聯(lián)分析挖掘林木抗病基因位點(diǎn)的研究極少見報道，僅有火炬松（Pinus taeda）-抗鐮刀球菌（Fusarium circinatum）[9]、抗脂潰瘍?。‵usarium circinatum）[10]，石榴（Punica granatum）-抗白葉枯?。╔anthomonas axonopodis）[11]，巨尾桉（Eucalyptus grandis × E. urophylla）-抗銹病（Puccinia psidii）[12]等少數(shù)幾個樹種病害關(guān)聯(lián)分析。木麻黃不僅尚未開發(fā)足量的分子標(biāo)記，在分子輔助選育、遺傳圖譜構(gòu)建、抗病基因定位等分子生物學(xué)領(lǐng)域的研究嚴(yán)重滯后，木麻黃青枯病抗性相關(guān)分子標(biāo)記研究國內(nèi)外也極少見報道。本研究利用木麻黃223 個EST-SSR 標(biāo)記，對我國現(xiàn)有木麻黃種質(zhì)資源青枯病抗性進(jìn)行關(guān)聯(lián)分析，旨在挖掘與木麻黃青枯病抗性密切相關(guān)的標(biāo)記位點(diǎn)或具有特定功能的基因位點(diǎn)，為分子標(biāo)記輔助育種提供有益參考，加速抗病種質(zhì)資源的選育進(jìn)程。

1 材料與方法

1.1 試驗材料

供試53 份木麻黃種質(zhì)來源于廣東、福建、海南，是廣東、福建、海南各省根據(jù)不同育種目的選育并保存下來的無性系，由廣東省林業(yè)科學(xué)研究院培育。

供試材料青枯病抗性鑒定方法、病情指數(shù)、抗性分級等數(shù)據(jù)來源于許秀玉等[13]研究結(jié)果，為了統(tǒng)計分析方便，本文將病情指數(shù)換算成抗病指數(shù)，抗病指數(shù)=100-病情指數(shù)。供試木麻黃種質(zhì)材料信息見表1。由表1 可見，供試材料抗病指數(shù)變異范圍大，介于27.2～100.0 之間，基本呈連續(xù)正態(tài)分布，適合用于抗病關(guān)聯(lián)分析。

1.2 DNA 提取與標(biāo)記分型

采取木麻黃新鮮嫩枝，采用改良的 CTAB 法[14]提取基因組DNA。檢測濃度與純度后稀釋至100 ng/μL，作為PCR 擴(kuò)增模板。

選取雜交、30、41、45、W6、G1、平5 和X1等8 個親緣關(guān)系較遠(yuǎn)、抗病指數(shù)大于95 的無性系構(gòu)建抗病基因池；選取2、34、37、59、82、G88、K18 和C7 等8 個親緣關(guān)系較遠(yuǎn)、抗病指數(shù)小于50的無性系構(gòu)建感病基因池。利用223 個木麻黃ESTSSR 標(biāo)記[15]對抗、感池進(jìn)行標(biāo)記分型與篩選。初步篩選標(biāo)記方法[16]：比較兩個基因池的標(biāo)記-等位基因頻率，相差80%以上的標(biāo)記入選，篩選出用于關(guān)聯(lián)分析的標(biāo)記。

反應(yīng)體系、分型過程及參試223 個標(biāo)記具體信息參見Xu 等[15]研究結(jié)果。

1.3 數(shù)據(jù)分析

1.3.1 群體結(jié)構(gòu)及關(guān)聯(lián)分析 標(biāo)記分型與篩選采用GenAlEx 6.4.1 檢測[17]，篩選出在抗感基因池中能差異表達(dá)且不偏離哈-溫平衡的標(biāo)記用于關(guān)聯(lián)分析。利用Structure 2.3.4[18]軟件分析木麻黃供試材料群體結(jié)構(gòu)，設(shè)定亞群數(shù)（K）的取值范圍1～16，將MCMC（Markov Chain Monte Carlo）的Length of Burn-in Period（不作數(shù)迭代）設(shè)為 100 000 次，每個K 值重復(fù)運(yùn)行10 次，其余參數(shù)采用默認(rèn)設(shè)置。根據(jù) lnP(D)計算ΔK，以ΔK 最大原則選取最適合的亞群數(shù)（K），并計算各參試材料的基因型屬于第K 亞群的概率（Q 值）。利用SPAGe-Di1-5a 軟件[19]獲取供試材料個體間親緣關(guān)系系數(shù)K 矩陣?；赒+K 模型，選擇TASSEL 3.0 軟件[19]的MLM 模型，確定P＜0.05 時的關(guān)聯(lián)位點(diǎn)，計算標(biāo)記對表型變異的解釋率（R2）。

表1 供試無性系及抗病指數(shù)Tab.1 The clones and the disease index

1.3.2 等位變異表型效應(yīng)分析 采用無效等位基因（null allele）法[20]計算各增效/減效等位變異表型效應(yīng)值（Ai），確定關(guān)聯(lián)標(biāo)記的增效/減效等位變異，將某一標(biāo)記增效/減效等位變異的平均效應(yīng)作為該標(biāo)記的增效/減效應(yīng)值（AAE），計算表型效應(yīng)比（AN）。Ai=∑xij/ ni-∑Nk/ nk；AAE= ∑ac/ nc；AN=AAE / ( ∑Nk/ nk) × 100%。（xij：第j 個具有i 等位片段個體的表型值；ni：具有i 等位片段的個體數(shù)；Nk：第k 個具有無效等位基因個體的表型值；nk：具有無效等位基因的個體數(shù)；ac：關(guān)聯(lián)位點(diǎn)內(nèi)第c 個增效/減效等位變異表型效應(yīng)值；nc：關(guān)聯(lián)位點(diǎn)內(nèi)增效/減效等位變 異數(shù)。）

2 結(jié)果與分析

2.1 木麻黃EST-SSR 標(biāo)記分型與篩選

檢測發(fā)現(xiàn)，223 個木麻黃EST-SSR 標(biāo)記中，共有49 個標(biāo)記在抗、感基因池中等位基因頻率差異達(dá)80%以上，且CASeSSR244 與CASeSSR064兩個標(biāo)記偏離哈-溫平衡（表2），因此選擇其余47 個標(biāo)記作為關(guān)聯(lián)分析標(biāo)記，利用這47 個標(biāo)記對參試材料進(jìn)行群體結(jié)構(gòu)分析與親緣關(guān)系系數(shù)計算。

2.2 群體結(jié)構(gòu)分析

根據(jù)模型后驗概率lnP(D)計算ΔK，通過ΔK 來確定亞群數(shù)（K），以ΔK 值最大時的K 值即為最佳亞群數(shù)。由圖1A 可知，當(dāng)K=2，ΔK 值最大，即參試群體明顯存在2 個亞群，表明將53份木麻黃種質(zhì)資源分為2 個組群時，具有最大的后驗概率ln P(D)，組群間具有最大的遺傳差異，組群內(nèi)植株具有相對一致的遺傳背景。Q 值計算結(jié)果發(fā)現(xiàn)：所有參試材料在各自所屬組群中的Q值全部大于0.5，這表明所有參試材料都可歸入2個組群中某一個，群體結(jié)構(gòu)簡單，遺傳組分相對單一，有利于減少關(guān)聯(lián)分析的假陽性，提高關(guān)聯(lián)分析效果。

表2 49 個標(biāo)記哈-溫平衡卡方檢驗結(jié)果Tab.2 Results of 49 EST-SSRs Chi-square tests for Hardy-Weinberg equilibrium

構(gòu)建的群體結(jié)構(gòu)圖如圖1B，亞群1 包含34個無性系，亞群2 包含19 個無性系，來源不同的參試材料在不同的亞群中均有分布，這表明亞群的劃分與地理來源之間沒有必然聯(lián)系，這與前期木麻黃種質(zhì)資源親緣關(guān)系分析結(jié)果相一致[21]。

2.3 關(guān)聯(lián)分析

由表3 關(guān)聯(lián)分析結(jié)果可見：有10 個標(biāo)記與青枯病抗性顯著相關(guān)（P ＜ 0.05），相關(guān)標(biāo)記對表型變異的解釋率較高，介于12.6%～60.1%之間。其中有5 個標(biāo)記具有增效效應(yīng)，標(biāo)記CASeSSR359、CASeSSR052 與CASeSSR056 增效效應(yīng)值最大，分別為19.1、18.1、11.7，其相對的表型變異解釋率分別為19.1%、18.1%、43.7%，推斷這些標(biāo)記與某個貢獻(xiàn)率較大的抗病位點(diǎn)連鎖的可能性較大，可利用這些標(biāo)記對木麻黃種質(zhì)資源的青枯病抗性潛力做出評價。除了CASeSSR359、CASeSSR052標(biāo)記外的其余8 個標(biāo)記均具有減效效應(yīng)，CASeSSR022、CASeSSR242、CASeSSR176 減效效應(yīng)值最大，分別為32.2，22.0，21.4，相應(yīng)的表型變異解釋率分別為35.9%，20.7%，12.6%，初步推斷這些標(biāo)記與感病位點(diǎn)連鎖的可能性較大，可利用這些標(biāo)記淘汰易感植物材料。

2.4 位點(diǎn)和等位片段的表型效應(yīng)

由表4 可見，具有增效表型效應(yīng)等位變異的 關(guān) 聯(lián) 位 點(diǎn) 有CASeSSR052、CASeSSR056、CASeSSR241、CASeSSR263 及CASeSSR359，其表型效應(yīng)范圍分別為：2.9～37.1、6.1～17.3、2.0～5.3、2.1～13.8、 11.5～27.3，其優(yōu)異增效等位變異分別為：CASeSSR052-255 bp（PE=37.1）、CASeSSR052-237 bp（PE=37.1）、CASeSSR359-370 bp（PE=24.0）、CASeSSR359-386bp（PE=27.3）。具有減效表型效應(yīng)的關(guān)聯(lián)位點(diǎn)有CASeSSR022、CASeSSR056、CASeSSR176、CASeSSR241、CASeSSR242、CASeSSR253、CASeSSR258、CASeSSR263，其中CASeSSR022-196 bp、 CASeSSR022-190 bp、CASeSSR022-188 bp、CASeSSR253-202 bp、CASeSSR242-227 bp、CASeSSR253-186 bp、CASeSSR176-204 bp 的等位變異減效表型效應(yīng)均高于20.0。育種中，可優(yōu)先利用表型效應(yīng)值大的優(yōu)異等位變異，利用優(yōu)異增效等位變異進(jìn)行優(yōu)質(zhì)親本及抗病植物材料的篩選，為后續(xù)標(biāo)記輔助育種及雜交育種提供有益參考。

圖1 木麻黃53 個無性系群體結(jié)構(gòu)Fig.1 Population structure of Casuarinaceae

表3 關(guān)聯(lián)標(biāo)記的平均增效效應(yīng)、減效效應(yīng)Tab.3 EST-SSR associated with R. solanacearum resistance could be detected with MLM model in Casuarinaceae and average positive (negative) allele effect

表4 木麻黃無性系青枯病抗病關(guān)聯(lián)位點(diǎn)及其等位變異Tab.4 Associated loci and phenotypic effect for R. solanacearum resistance

3 結(jié)論與討論

木麻黃分子生物學(xué)研究基礎(chǔ)薄弱，早期利用農(nóng)桿菌開展了木麻黃基因工程研究[22]，利用FISSR、ISSR、SSR、SCAR 等標(biāo)記開展了木麻黃遺傳多樣性[23]、群體結(jié)構(gòu)[23-24]、群體親緣關(guān)系[21]、木麻黃物種鑒定[25]等少量分子育種研究，2018 年開 發(fā) 了223 個 木 麻 黃EST-SSR 標(biāo) 記[15]，2019 年首次公布了短枝木麻黃（Casuarina equisetifolia）的基因組數(shù)據(jù)[3]，本研究在此基礎(chǔ)上首次開展了木麻黃青枯病抗病性狀與分子標(biāo)記的關(guān)聯(lián)分析，共發(fā)掘出CASeSSR241、CASeSSR263、CASeSSR056 等10 個與木麻黃青枯病抗病性狀顯著關(guān)聯(lián)的EST-SSR 位點(diǎn)及CASeSSR052-255 bp、CASeSSR052-237 bp、CASeSSR359-370 bp、CASeSSR359-386 bp 等優(yōu)異增效等位變異，為今后預(yù)測控制性狀的主效QTL、新基因發(fā)掘、標(biāo)記輔助育種、雜交育種等研究提供參考。其次，CASeSSR253-202 bp、CASeSSR022-190 bp、CASeSSR022-196 bp 等多個等位片段對木麻黃青枯病抗性具有較高的減效效應(yīng)，利用這些等位變異可以淘汰木麻黃感病材料。

本研究從223 個木麻黃EST-SSR 標(biāo)記中篩選得到10 個與青枯病抗性相關(guān)的標(biāo)記，關(guān)聯(lián)標(biāo)記數(shù)量偏少，這可能由于與木麻黃青枯病抗性相關(guān)的基因原本就不多，也可能由于所用的EST-SSR 標(biāo)記數(shù)量不夠多，今后可擴(kuò)大標(biāo)記數(shù)量與種類以獲得更多的關(guān)聯(lián)位點(diǎn)。本研究篩選出的10 個關(guān)聯(lián)標(biāo)記表型變異解釋率介于12.6%～60.1%，高于大多數(shù)植物抗病蟲害關(guān)聯(lián)標(biāo)記，如辣椒（Capsicum annuum）疫?。≒hytophthora capsici Leonian）抗性關(guān)聯(lián)分析中表型變異解釋率2.98%～16.05%[26]，煙草（Nicotiana tabacum）青枯病抗性關(guān)聯(lián)分析中表型變異解釋率5%～14%[19]，陸地棉（Gossypium hirsutum）黃萎病（Verticillium dahliae）抗性關(guān)聯(lián)分析表型變異解釋率3%～13%[27]等，這可能由于本次試驗群體數(shù)量少造成的，也可能由于這些標(biāo)記位于基因編碼區(qū)或基因表達(dá)調(diào)控區(qū)，是青枯病抗性性狀重要位點(diǎn)。

木麻黃為外來樹種，本研究參試群體為我國海南、福建、廣東各省根據(jù)抗風(fēng)、速生、抗旱等不同育種目的選育出來的，為我國現(xiàn)存的木麻黃重要種質(zhì)資源，因此群體數(shù)量少。利用這個群體進(jìn)行關(guān)聯(lián)分析難免存在著LD 偏差及表型取樣誤差[26]，為了提高關(guān)聯(lián)分析的準(zhǔn)確性，今后要進(jìn)一步加強(qiáng)木麻黃育種群體的收集、積累、選育及雜交品系的創(chuàng)制，不斷擴(kuò)大用于木麻黃關(guān)聯(lián)分析的群體數(shù)量，對獲得的標(biāo)記位點(diǎn)與優(yōu)異等位變異進(jìn)一步驗證。其次，本研究中木麻黃青枯病抗性表型數(shù)據(jù)為室內(nèi)接種試驗所得，而林木對青枯病抗性的強(qiáng)弱除了受基因型影響外，還受自身生長勢、侵染方式及環(huán)境條件等因素的影響，今后還應(yīng)開展多地點(diǎn)田間試驗，獲得不同地點(diǎn)、不同生境條件下，不同木麻黃資源對青枯菌的抗病性能，研究環(huán)境與基因型的互作，以得到更多更準(zhǔn)確的抗性相關(guān)位點(diǎn)與優(yōu)異等位變異。