馬鈴薯品種（系）田間晚疫病抗性評價和全基因組遺傳多樣性分析

李曉川，王朝海，周平，馬維，吳瑞，宋治豪，梅艷

畢節市農業科學研究所，貴州畢節 551700

0 引言

【研究意義】晚疫病由致病疫霉菌（（Mont.）de Bary）引起，是馬鈴薯（L.）生產的第一大病害。貴州省所處的烏蒙山區由于氣候冷涼潮濕，易大規模爆發晚疫病。致病疫霉菌存在高度的變異性，簡單的單個抗病基因導入栽培種所帶來的晚疫病抗性，很容易被變異的致病疫霉菌克服。因此，研究馬鈴薯種質的遺傳多樣性及抗病種質基因組中可能影響抗性的遺傳區段，對馬鈴薯晚疫病抗性種質創新與利用具有重要意義。【前人研究進展】關于馬鈴薯晚疫病抗性育種的研究，多集中在克隆具有晚疫病致病疫霉菌小種專化抗性的R（Resistance）基因上，目前，已克隆獲得多個R基因，均是NBS-LRR類基因。但無論是具有專一的小種專化抗性的R基因（、、和），還是具有廣譜的小種專化抗性的R基因，將單個抗病基因導入栽培種后所具有的晚疫病抗性，在商業化種植后，均易被變異的致病疫霉菌克服。如，將廣譜抗性RB基因導入培育的栽培種Biogold，田間種植1年后，便被晚疫病菌克服。馬鈴薯晚疫病的田間抗性被認為具有相對的持久性，具有非小種專化性抗性，其形成的遺傳學基礎是由多個基因控制，認為多個R基因的堆疊可以增強馬鈴薯晚疫病田間抗性表現。【本研究切入點】西南烏蒙山區是中國馬鈴薯晚疫病發病最為嚴重的地區，而利用此環境條件，對馬鈴薯種質資源進行田間晚疫病抗性評價，以及隨后進行全基因組遺傳多樣性分析的研究相對較少。【擬解決的關鍵問題】本研究通過在晚疫病多發且發病嚴重地區對馬鈴薯種質進行田間晚疫病抗性鑒定，獲得抗病種質資源，運用現代高通量測序技術分析種質的遺傳多樣性及抗病種質基因組內可能影響抗性的遺傳區段，為馬鈴薯晚疫病抗性種質創新與利用提供理論依據。

1 材料與方法

1.1 馬鈴薯田間晚疫病抗性材料篩選及取樣

參加馬鈴薯田間晚疫病抗性篩選的材料，包括自家選育的中高代品系432份以及從各地搜集的馬鈴薯品種（系）458份（電子附表1）。于2018—2021年分別在貴州省畢節市威寧縣草海鎮鄭家營村、卯關村，威寧縣雪山鎮團崗村、安家營村、法地村以及七星關區朱昌鎮王家沖村內多點種植，每個品種（系）種植50株，種植過程不噴灑晚疫病防治藥劑，進行田間晚疫病抗性鑒定。抗性材料篩選：以在晚疫病爆發一個月后，植株依然能夠存活的馬鈴薯品種（系）為抗病材料（貴州畢節一般在3月中下旬種植馬鈴薯，4月中下旬出苗，6月中下旬晚疫病爆發，7月中下旬可達到當地主栽中晚熟馬鈴薯品種生長期90—100 d的要求）；以在發病10 d左右大部分植株死亡的馬鈴薯品種為感病材料；其余為半抗病材料。挑選待測序感病材料時，參考歐洲馬鈴薯栽培種數據庫（http://www.europotato.org/quick_search.php）關于晚疫病抗性及其他農藝性狀的數據，同時參考馬鈴薯系譜數據庫（http://www.plantbreeding.wur.nl/potatopedigree/index.html）的系譜學數據。取晚疫病抗感材料葉片，分別進行混合，將樣品經液氮速凍后保存備用。

1.2 dd-RAD簡化基因組測序

將樣品送派森諾生物科技公司構建 dd-RAD（double digest restriction site-associated DNA sequencing）簡化基因組測序文庫。經Illumina NovaSeq高通量平臺進行2×150 bp的雙端測序。采用FastQC（http://www.bioinformatics.babraham.ac.uk/projects/fastqc）對原始數據進行質量控制，使用fastp對數據進行過濾。采用 BWA-MEM 程序將過濾后的高質量數據比對馬鈴薯參考基因組（DM v6.1），采用GATK軟件包來進行SNP的檢測及過濾，采用ANNOVAR軟件對SNP位點進行注釋。

1.3 遺傳多樣性分析

采用GCTA軟件對SNP數據（去除MAF＜0.05的SNPs）進行主成分分析。采用fastTree軟件中的Maximum likelihood算法構建系統發育樹，對系統發育樹分支的可靠性進行驗證（bootstrap，1 000 replications）。采用 Admixture軟件對分型獲得的SNP信息進行群體遺傳結構分析，設置K=2—20（即假設存在2—20個群體），計算cross-validation（cv）error值，模型選擇為混合模型，其余參數為軟件的默認設置。根據不同K的cv值選擇群體個數。采用stacks程序包中的populations命令進行群體遺傳多樣性參數的計算。

1.4 選擇性消除分析

采用 vcftools程序進行選擇性消除參數的計算。利用代表不同田間晚疫病抗性群體的SNP標記信息，以20 kb為窗口，5 kb為步長分別計算π值，并計算代表不同窗口和步長抗病種質和感病種質 π值的比值。同時，選取相同長度的窗口和步長計算不同晚疫病抗性群體間的Fst值。最終獲得不同窗口和步長的遺傳區段的π值比值和Fst值。其中，抗病和感病種質之間π值比值小，代表相對于感病種質，抗病種質的這段染色體區域核苷酸多樣性低，是受選擇的遺傳區段，選取π值比值最小的5%的遺傳區段，并選取這些遺傳區段中Fst值最大10%的遺傳區段進行后續分析，即選取占基因組0.5%的遺傳區段。通過基因本體論分析（gene ontology，GO）預測遺傳區段內的基因。基因的氨基酸序列利用 Clustal Omega進行多序列比對。進化樹利用MEGA6軟件的Neighbor-Joining參數通過計算1 000次生成樹圖。

1.5 馬鈴薯抗晚疫病種質全基因組關聯分析

以101個抗病種質和21個感病種質的群體SNP信息，利用GEMMA 0.98.1軟件的一般線性模型進行全基因組關聯分析，得到各SNP位點與性狀關聯程度的值，并計算-log()，-log()值越大，表示 SNP位點與性狀的關聯程度越強。利用 CMplot繪制-log()和期望-log()的QQ圖以及-log()的曼哈頓圖。選擇-log()＞6的SNP位點周圍50 kb的基因組序列，通過基因本體論分析預測遺傳區段內的基因。

2 結果

2.1 馬鈴薯田間晚疫病抗性鑒定

將收集的馬鈴薯品種（系）經2018—2021年在貴州省畢節市威寧縣和七星關區內6個地點種植，進行田間晚疫病抗性鑒定。有101個品種（系），在晚疫病爆發一個月后依然能夠存活，基本可以滿足當地主栽的中晚熟馬鈴薯品種生長期90—100 d的要求，確定為抗晚疫病品種。同時以在晚疫病爆發后10 d左右大部分植株死亡為標準，共鑒定得到感病品種175個，進一步從中挑選了薯條/薯片加工品種11個和優秀親本9個（兩者間有重疊）及擁有其他優良性狀的品種4個，共21個品種進行下一步分析（表1）。

表1 122份馬鈴薯種質資源和2d-RAD簡化基因組測序情況Table 1 The 122 potato germplasm and 2d-RAD simplified genome sequencing results

續表1 Continued table 1

續表1 Continued table 1

續表1 Continued table 1

2.2 馬鈴薯種質的簡化基因組測序

將122個（101個抗晚疫病馬鈴薯品種（系）以及21個感病品種）樣本進行dd-RAD簡化基因組測序。每個樣本各產生957 120 660—2 571 656 553 bp的測序數據。過濾掉低質量測序片段（Reads）和不合格序列后，每個樣本各產生6 750 190—18 433 648個高質量測序片段，其中，6 701 577—18 323 584個測序片段與參考基因組匹配（DM v6.1），占全部測序片段的99.28%—99.48%，覆蓋參考基因組15.32%—27.14%（表1）。共檢測到8 697 602個SNP，其中，在基因外顯子區域的有397 244個SNP，內含子區域的有970 743個SNP，基因上下游1 kb內有499 056個SNP，基因間區域有6 650 230個SNP（表2）。馬鈴薯12條染色體上，平均每1 Mb得到的SNP在2 828個以上（圖1）。

表2 122份馬鈴薯種質資源SNP基因分型情況Table 2 SNP identification in the 122 potato germplasm

圖1 2d-RAD簡化基因組測序分型得到的SNP密度熱圖Fig. 1 SNP density heat map from 2d-RAD simplified genome sequencing

2.3 馬鈴薯種質的遺傳多樣性分析

利用檢測到的SNP信息對122份馬鈴薯種質進行種群分析。首先，假定這些種質資源由2—20個群體（population）組成（即K=2—20），利用Admixture軟件分析種群結構，當K=6時，cross-validation（cv）error值最小，但當K=3—9時，cv值差別不大（圖2-A）。說明這些種質資源可能是由3—9個群體組成，但最可能由6個群體（population）組成（圖2-B）。同時，將群體的SNP信息簡化為3個維度的信息，進行主成分分析（PCA），將結果根據種群結構分析得到的 6個群體，繪制三維散點圖（圖2-C），顯示雖然各群體成員有聚集在一起的趨勢，但界限并不明顯。進一步繪制遺傳進化樹（圖2-D），顯示這些種質資源可以分成6個群體（population），群體內又有親緣關系更近的資源聚集成亞群（sub-population），如，在群體Ⅱ中來源于云薯的品種（系）聚集在亞群3中，幾個著名的加工品種（如Atlantic、Shepody等）聚集在群體Ⅳ中，親本來源于國際馬鈴薯中心的幾個品系聚集在群體Ⅵ內的亞群2和亞群3中。同時，發現種群結構的劃分與資源的晚疫病抗性并不相關。

圖2 122份馬鈴薯種質資源的種群分析Fig. 2 Population structure analysis of the 122 potato germplasm

進一步利用種群結構信息進行種質遺傳多樣性分析。在6個群體內平均核苷酸多樣性值（π），范圍為0.2055—0.2572（表3），說明群體內樣本間的遺傳差異較大，范圍為0.156909—0.187336（表4），說明群體間遺傳分化程度也較大。6個群體內的期望雜合度（e）范圍為0.187—0.2297，觀測雜合度（o）范圍為0.0829—0.1186（表3），觀測雜合度均小于期望雜合度。同時6個群體內的近交系數（Fis）范圍為0.2412—0.3554（表3），說明6個群體內存在近交的現象。通過進行遺傳多樣性參數的計算，發現盡管101份馬鈴薯抗晚疫病種質和21份感病種質存在較豐富單核苷酸遺傳多樣性，但在選育這些種質資源的過程中存在著親緣關系較近親本相互雜交的情況。

表3 群體內的遺傳多樣性參數Table 3 Genetic diversity parameters within the population

表4 群體間群體分化指數（Fst）的矩陣Table 4 The matrix of fixation index (Fst) among populations

2.4 馬鈴薯種質選擇性消除分析

以分離基因組內可能影響馬鈴薯晚疫病抗病性狀的遺傳區段進行選擇性消除分析。選取抗病和感病種質之間 π值比值最小的 5%的遺傳區段，并進一步選取這些遺傳區段中Fst值最大的10%的遺傳區段進行后續分析，共獲得993個遺傳區段，合并相鄰的遺傳區段后，共有745個遺傳區段，其中，最大的遺傳區段長75 kb（圖3-A）。這些遺傳區段中共包含507個基因，通過GO基因功能預測，膜的組成部分有最多數目的基因；與防御反應相關的有14個基因；許多基因預測擁有分子結合功能，如：ATP結合、ADP結合、DNA結合、氨基酸結合、鋅離子結合、鐵離子結合等（圖3-B）。在這些基因中，還發現了 4個擁有晚疫病抗病（R）基因結構域特性的NBS-LRR類基因（圖3-C），與已知的晚疫病抗病基因進行序列比對，發現與、的序列相近；、和與/、和的序列相近。

圖3 選擇性清除分析Fig. 3 Selective sweep analysis

2.5 馬鈴薯種質全基因組關聯分析

以101個抗病種質和21個感病種質的群體SNP信息，進行全基因組關聯分析，得到各SNP位點與性狀關聯程度的-log()。在QQ圖右側，部分SNP位點的觀測-log()與期望-log()發生偏離（圖4-B），說明與總體SNP位點不同，這部分SNP的觀測-log()不符合正態分布，表明這部分SNP位點附近的染色體區域可能受到了選擇。在曼哈頓圖中，與晚疫病抗性性狀關聯最強的9個SNP位點（-log()＞6）分布于馬鈴薯的不同染色體上（圖4-A）。在這些SNP位點周圍50 kb的基因組范圍內，共有69個基因，其中，41個基因可以依據細胞成分、生物過程及分子功能通過GO預測得到基因功能。在41個基因中，12個基因編碼的蛋白預測為質體的組成成分，15個基因預測參與到應激反應，12個基因預測參與清除過氧化物自由基的過程，4個基因預測參與半胱氨酸合成，3個基因預測參與跨膜運輸（圖4-C）。

圖4 全基因組關聯分析Fig. 4 Genome-wide association analysis

3 討論

3.1 dd-RAD簡化基因組測序是馬鈴薯基因分型的有效手段

構建dd-RAD簡化基因組文庫首先通過對參考基因組進行分析，選取合適的限制性內切酶對基因組DNA進行酶切，因此，在樣本間基因組DNA酶切位置基本相同，并且可以在基因組內獲得位置間隔相對均勻的待測序片段。相比較于一些利用物理手段進行基因組 DNA片段化的文庫構建方法，dd-RAD是一種可預測的文庫構建方法。文庫構建完成后，通過2×150 bp的雙端測序，得到150 bp的測序片段，相對于其他利用限制性內切酶構建測序文庫的方法，如：2b-RAD只得到36 bp的測序片段，dd-RAD得到的測序片段長度更長，因此，將測序片段比對到參考基因組上時，比對結果更加可靠，分型結果也更加準確。本研究利用 dd-RAD簡化基因組測序，對122份馬鈴薯種質資源進行基因分型，測序片段覆蓋了參考基因組 15.32%—27.14%，共檢測到8 697 602個SNP標記。在馬鈴薯12條染色體上，平均每1 Mb含有2 828個SNP以上（圖1）。馬鈴薯的連鎖不平衡（LD）衰減值最小為600 kb，平均在2—4 Mb，表明本研究的SNP密度足以滿足要求。

3.2 群體結構分析揭示馬鈴薯種質資源之間的遺傳關系

群體結構分析能夠揭示馬鈴薯種質資源之間的遺傳關系。研究表明，通過群體結構分析，可以將普通栽培種、安第斯栽培種及野生種區分開。并且，部分來源相同、育種地相近的種質，可以在群體結構分析中聚集在一起。如，DUAN等研究表明大部分中薯系列品種聚集在一個類群中。WANG等研究表明不同時間引入中國的國外品種、來源于國際馬鈴薯中心的種質資源、以及中國不同地區選育的馬鈴薯品種分別聚集在不同的亞群。與這些研究類似，本研究來源于云薯的種質聚集在群體Ⅱ的亞群3中，親本來源于國際馬鈴薯中心的幾個品系聚集在群體Ⅵ的亞群 2和亞群3中（圖2-D)。

群體結構分析可以在一定程度上區分薯條、薯片加工品種與其他普通栽培種，本研究幾個著名的加工品種（如大西洋（Atlantic）、夏坡蒂（Shepody）等）聚集在群體Ⅳ中（圖2-D）。但種質種群結構的劃分與種質的晚疫病抗性并不相關（圖2-D），這可能是由于晚疫病抗性基因在群體里屬于稀有位點，它們的選擇與否與群體結構之間沒有必然聯系。

同時，為了緩解馬鈴薯遺傳基礎狹窄的問題，可以根據群體結構，分析種質資源之間的親緣關系，并在雜交育種中選擇親本親緣關系較遠的種質。例如，本研究所涉及的幾個易感晚疫病品種Agria、Desiree、Innovator、Felsina、Nicola和Superior，它們除了擁有優良的農藝性狀外，更重要的是它們具有較高的一般配合力，是國際范圍內近期作為直接親本育成品種數最多的幾個品種，可以在親本選配中選擇群體Ⅵ內的種質雜交。同時，本研究中所涉及的種質，尤其是感病種質較少，可以利用在不同研究中相同種質所處的種群結構位置，分析其他研究中的種質與本研究種質間的親緣關系。如通過相同品種（Favorite、Desiree、中薯17號）分析與189份不同晚疫病抗性的普通栽培種、安第斯栽培種及野生種之間的親緣關系。通過相同的云薯系列品種，分析與292份不同時間引入中國的國外品種、來源于國際馬鈴薯中心的種質資源、以及中國不同地區選育的馬鈴薯品種之間的關系。通過相同品種（Atlantic、Favorita、Katahdin）分析與209份中國育成的品種之間的親緣關系。

3.3 馬鈴薯田間抗晚疫病種質資源擁有較豐富的單核苷酸多態性但在其選育過程中存在近交現象

馬鈴薯作為營養器官繁殖的作物，其品種選育經歷的有性雜交過程十分有限。系譜學研究顯示，有記錄的馬鈴薯雜交育種始于一個智利普通栽培品種Rough Purple Chili，它的后代與歐洲品種雜交后，表現突出，從此馬鈴薯雜交育種開始得到重視并廣泛應用。美國馬鈴薯協會列出北美177個應用于生產的品種，其中，156個品種含有Rough Purple Chili的血緣。在全世界育成的有記錄的4 397個品種中，15個直接親本育成的品種數占總品種數的 14.9%；523個有記錄的中國品種中，8個直接親本育成的品種數占總品種數的27.34%。因此，馬鈴薯存在著種質資源的遺傳基礎狹窄的問題。在本研究中，通過利用122份馬鈴薯種質的群體 SNP信息，可以將這些種質資源分成6個群體（圖2）。6個群體內的觀測雜合度均小于期望雜合度，并且 6個群體內的近交系數（Fis）為 0.2412—0.3554（表3），說明選育這些種質的過程中存在近交的現象，也反映了這些種質資源的遺傳基礎狹窄的問題。同時，馬鈴薯作為基因組高度雜合的多倍體，本身擁有較豐富的單核苷酸多態性。在6個群體內，反映群體內樣本間平均遺傳多樣性差異程度的平均核苷酸多樣性值（π）為0.2055—0.2572（表3），反映群體間的分化程度群體分化指數（Fst）為 0.156909—0.187336（表4），說明田間抗晚疫病種質資源擁有較豐富的單核苷酸多態性。

3.4 選擇性清除分析和關聯分析有助于分離影響農藝性狀染色體區段

選擇性清除分析利用在受選擇的遺傳區段中，遺傳多樣性降低并且遺傳分化升高，來分離影響農藝性狀的染色體區段。本研究以20 kb為窗口5 kb為步長，在基因組相同位置，分別計算不同晚疫病抗性馬鈴薯種質間平均核苷酸多樣性（π）值比值和群體分化指數（Fst）值。相較于其他類似的研究，本研究選擇了相對較小的窗口計算π值和Fst值，一方面是因為簡化基因組測序只對部分基因組進行了基因分型，過大的窗口可能平均化已測序和未測序遺傳區段的單核苷酸多態性信息，并導致結果不能反映該區段的真實情況；另一方面本研究種群結構的劃分與其晚疫病抗性并不相關，選擇較小的窗口更能精細化區分遺傳區段的受選擇情況。以 π值比值極小及Fst值極大，選擇可能影響晚疫病抗性的745個遺傳區段進一步分析，這些遺傳區段中共包含507個基因，其中，有14個基因可能參與防御反應，更為重要的是有4個類基因，1個與的序列相近，另外3個與擁有廣譜晚疫病抗性、和的序列相近（圖3-C），其中，/克隆自馬鈴薯野生種，是克隆的第一個具有廣譜晚疫病抗性的基因，和克隆自馬鈴薯野生種，是的同源基因。說明選擇性清除分析有助于分離影響晚疫病性狀的遺傳區段。

關聯分離利用連鎖不平衡，使用自然群體挖掘與性狀相關聯的基因，在馬鈴薯中也有廣泛的應用。例如，通過分析抗黃瘟病（）基因的等位基因座，在雙倍體及四倍體馬鈴薯群體中，對抗黃瘟病性狀與SNP標記之間進行了關聯分析，克隆出了一個抗黃瘟病的基因。同樣也有研究利用全基因組關聯分析對馬鈴薯產量、淀粉含量、抗旱等性狀進行了研究。在本研究中，得到了 9個可能與晚疫病性狀相關聯的SNP位點，在它們周圍的50 kb的基因組范圍內，有15個基因預測參與應激反應，有12個基因預測參與清除過氧化物自由基，這些基因可能參與了晚疫病侵染后的細胞的修復過程（圖4-C）。

4 結論

dd-RAD簡化基因組測序可以在馬鈴薯基因組中分型到數量較多分布相對均勻的SNP標記。馬鈴薯田間抗晚疫病種質資源擁有較豐富的單核苷酸多態性但在其選育過程中存在近交的現象。選擇性清除和關聯分析有助于分離影響晚疫病抗病性狀的遺傳區段。