廣東省與江蘇省松材線蟲種群遺傳結構差異分析*

丁曉磊 汪青桐 林司曦 趙瑞文 張 悅 葉建仁

(南方現代林業協同創新中心 南京林業大學林學院 南京 210037)

松材線蟲病是一種由松材線蟲(Bursaphelenchusxylophilus)引起松樹毀滅性死亡的森林病害。松材線蟲是一種生活史復雜且傳播迅速的入侵物種，自1982年在南京中山陵的黑松(Pinusthunbergii)上發現該病以來，疫情迅速蔓延(程瑚瑞等， 1983)。根據國家林業與草原局發布的2022年第6號公告， 我國大陸地區松材線蟲病已擴散至731個縣級行政區由南向北不斷擴散。20多年來松材線蟲病是我國發生最嚴重和造成損失最大的林業病害，絕大多數省區均屬于松材線蟲適生區(葉建仁， 2019)，而目前松材線蟲病在我國的傳播路徑尚不明確。有學者推測廣東省是松材線蟲在我國最初的定殖和擴散中心，江蘇省是新的擴散中心(謝丙炎等， 2009)。因此，不同發病地區間的種群結構和遺傳變化規律有待深入研究。

疫源追溯是明確和掌握松材線蟲病人為傳播路徑的重要手段。在開展疫源追溯的過程中，首先要了解松材線蟲在中國的種群分布及其遺傳結構變化規律。通過研究各種群間不同個體的基因組差異，可以對不同個體進行遺傳變異水平上的分類(黃文達等， 2010)。關于松材線蟲種群遺傳變異和地理來源的關系國內外已有許多研究，大多是基于DNA片段ISSR、AFLP和RFLP長度和酶切位點等的變異程度來進行相關的比較分析(Aikawaetal.， 2013; Guetal.， 2011; Jungetal.， 2010; 許峻榮等， 2014)。這些研究在闡明松材線蟲種群變異與地理分布及病害傳播進程相互關系中提供了重要參考。松林線蟲傳入我國已有40年，松材線蟲病主要依賴于人為因素傳播而呈現出跳躍式發展的特點，因此松材線蟲的種群分布特征越來越復雜(葉建仁等， 2012)，原有標記的靈敏性和分辨度已難以適應松材線蟲群體在不同地理區域間復雜的遺傳結構的變化，因此如何尋求高靈敏度和高分辨度的遺傳標記來研究當前松材線蟲在中國的種群結構與地理變化規律就顯得十分的重要。

SNP(Single nucleotide polymorphism)是基于第二代分子標記技術發展起來的第三代分子標記技術(Lander， 1996)，是由單核苷酸在基因組水平上的變異而產生的DNA序列多態性，包括堿基的轉換、顛換以及插入或缺失等(陳秋玲等， 2010)。在群體遺傳學中，SNP遺傳標記是發現物種和種群基因組多態性位點的方法(曾燕如等， 2003； Phillipetal.， 2004)。目前，SNP分子標記技術在人類基因組學(Baietal.， 2018)、農學(Woldegiorgisetal.， 2019)和動物學(Choietal.， 2016)等方面均廣泛應用，該技術在各領域的研究已經頗為成熟。本研究采用SNP分子標記方法，分析廣東與江蘇松材線蟲群體的遺傳多樣性，探究這2個代表性地區松材線蟲種群的遺傳變異和親緣關系，為松材線蟲病在中國的傳播與溯源提供參考。

1 材料與方法

1.1 蟲株收集及鑒定純化

采用貝爾曼漏斗法分離廣東省(GD)和江蘇省(JS)不同地區松材線蟲病疫木中的松林線蟲，所有疫木均由國家林業和草原局公告中的疫區所在地的林業部門現場采集，詳見表1。疫木分離獲得線蟲后參照謝輝的松材線蟲分類系統(謝輝， 2000)，對松材線蟲進行初步形態學鑒定，采用SCAR分子標記方法對松材線蟲進行分子生物學鑒定(陳鳳毛， 2005)。

對檢測后發現含松材線蟲的24棵疫木中分離挑取松材線蟲雌雄成蟲各15條左右，接入長好灰葡萄孢(Botrytiscinerea)的PDA平板中，28℃培養7天左右至松材線蟲取食完灰葡萄孢。使用貝爾曼漏斗法分離培養基中線蟲，進行形態學鑒定和分子生物學鑒定，確定為松材線蟲后使用體積分數0.05%硫酸鏈霉素浸泡8～10min，再使用無菌水洗滌3遍，放入4℃冰箱保存備用?；铙w線蟲蟲株均保存于南京林業大學森林病理實驗室松材線蟲蟲株資源庫。

表1 24株蟲株樣本來源

1.2 供試線蟲DNA提取

參照陳鳳毛等(2005)方法提取松材線蟲DNA。吸取200μL線蟲液置于1.5 mL Eppendorf 管中，加入200μL線蟲裂解液并混勻； 在65℃水浴中處理1 h，待冷卻至室溫后，加入等體積的酚抽提液，4℃，離心機13 000 r·min-1離心10min，取上清液； 用等體積酚∶氯仿∶異戊醇為25∶ 24∶1(V/V)混合液抽提，13 000 r·min-1離心10min，取上清液;用等體積氯仿： 異戊醇為24∶1(V/V)，13 000 r·min-1離心10min，取上清液； 在上清液中加入1/10體積的3mmol·L-1乙酸鈉(pH為4.6)和2倍體積的無水乙醇(-20℃)，于-80℃冷凍20min;13 000 r·min-1離心10min，沉淀用70%的乙醇(-20℃預冷)洗滌2次，65℃烘干10min； 加入100μL TE 重新懸浮沉淀，取其中5μL用無菌水稀釋200倍，用紫外光分光光度計測定DNA 質量濃度，再通過質量分數1%瓊脂糖凝膠電泳檢測DNA 質量，-20℃條件下保存，備用。松材線蟲蟲株DNA均在收集培養完成后進行提取并保存于南京林業大學森林病理實驗室松材線蟲蟲株DNA資源庫。

1.3 供試線蟲全基因組測序

將24個松材線蟲DNA進行Illumina高通量基因組測序(武漢未來組生物公司)，測序使用平臺為HiSeq 4000。基因組重測序方法為非鏈特異性、150 bp雙末端測序，測序深度>40×，每個樣本數據量大小為>8 G。

1.4 數據處理

1.4.1 下機數據處理及SNPs位點信息統計 1) 測序質量控制： 通過FastQC軟件，查看數據的基本信息、堿基的質量值、序列中duplication程度等基本信息。

2) 基因組比對分析： BWA (-M -t 12)、SAMtools (sort+index)和Picard (REMOVE_DUPLICATES=true)軟件將測序結果與松材線蟲全基因組進行比對，并去除可能存在的污染序列。

3) SNPs位點鑒定： 用Freebayes軟件，參數設置為“-u -C 5 -e 50 --standard-filters --min-coverage 10”，其他參數均為軟件默認值。

4) 樣本基因型分析與統計： 主要由VCFtools插件vcf-stat完成初步統計，之后利用Perl個性化腳本進行深度信息挖掘，并使用Prism進行統計作圖。

1.4.2 松材線蟲種群遺傳分化分析 對松材線蟲測序蟲株的SNP位點信息進行生物軟件分析，應用R語言包中的SNPRealte(https:∥bioconductor.org/packages/release/bioc/html/SNPRelate.html)軟件對等位基因頻率較低、連鎖不平衡(linkage disequilibrium)較高的SNPs位點進行過濾，同時繪制PCA主成分分析圖； 使用VCF-kit(https:∥vcf-kit.readthedocs.io/en/latest/)進行種群分化樹構建。

差異位點分析和過濾使用PLINK2(http:∥www.cog-genomics.org/plink2/)軟件進行挖掘，基因功能富集分析由agriGO網站在線完成(http:∥systemsbiology.cau.edu.cn/agriGOv2/)。

1.5 焦磷酸測序法驗證特異性SNP位點

應用PyroMark Assay Design Software 2.0軟件設計焦磷酸測序SNP分型引物。下游引物在5′端進行生物素標記。引物由南京金斯瑞生物科技公司合成(表2)。

反應體系： 25μL Etaq pre-Mix，DNA模板與正反向引物各2μL， 19μL ddH2O，總體系50μL。反應程序： 94℃預變性5min； 94℃ 變性 30 s，55℃ 退火 30 s，72℃ 延伸30 s，共30個循環； 最后72℃延伸7min。PCR產物用1%瓊脂糖凝膠電泳進行質量檢測。

表2 引物及其序列

按照QIAGEN PyroMark Q96 ID實時定量焦磷酸序列分析儀說明書對篩選出的SNP位點進行上機驗證。

2 結果與分析

2.1 蟲株分離純化

對疫木線蟲進行分離鑒定，最后純化獲得24株松材線蟲，廣東省和江蘇省各12株，具體蟲株樣品信息如表1。其中11株松材線蟲都分離自馬尾松，僅有1株來自江蘇省連云港市的松材線蟲分離自赤松。所有分離后的蟲株均經由南京林業大學自主研發的松材線蟲自動檢測系統鑒定，確定其為松材線蟲。

2.2 松材線蟲SNPs位點統計

通過生物信息學分析軟件，對24株松材線蟲蟲株進行SNPs位點統計(表3)。在所有蟲株中共發現15 100 281個SNP位點，其中廣東省蟲株具有12 990 503個SNP位點，江蘇省蟲株則具有2 109 778個SNP位點，廣東省蟲株的SNP位點總量遠大于江蘇省蟲株。從表3中可以看出，SNPs數量最多的是GD30，最少的是JS09，除GD09和GD12外，廣東省蟲株SNPs數量均多于江蘇省蟲株； 純合子數量最多的是GD02，最少的是JS09； 缺失SNP數量最多的為GD04，最少的是GD30； 特有SNP數量最多的是GD13，最少的是JS01?？傮w來看，廣東省蟲株的SNP數量和純合子數量明顯比江蘇省蟲株多，且差異性極顯著(P<0.01)(圖1)。綜上所述，除GD09和GD12外，廣東省蟲株的遺傳多樣性比江蘇省蟲株高，而GD09和GD12與江蘇省蟲株遺傳多樣性差異小。

表3 24個松材線蟲樣本中的SNPs位點統計

****：差異極顯著Significant difference (P<0.0001)

圖2 SNP突變基因型統計

圖3 聚類分析PCA圖及系統發育樹

2.3 突變基因型統計

所有蟲株種均有G→T、C→G、G→C、T→A、A→G、T→G、G→A、T→C、A→C、C→A、A→T、C→T，這12種基因型，但不同蟲株中具有體基因型出現頻率有差別。廣東蟲株中，除了GD09和GD12蟲株外，出現頻率最高的基因型為A→G、C→T、G→A、T→C這4種類型，而所有江蘇省蟲株出現頻率最高的基因型為A→G、C→G、G→C、T→C，與GD09和GD12兩蟲株相同(圖2)。廣東省蟲株除GD09和GD12兩蟲株外，突變基因型數量多于江蘇省蟲株，且主要突變基因型種類不同，表現明顯的遺傳差異性。

2.4 樣本親緣關系分析

對所有的蟲株進行主成分分析(PCA)，并繪制系統進化樹。JS蟲株全部聚為一類，且JS蟲株互相之間遺傳距離小，蟲株之間未發生明顯的遺傳分化。GD蟲株之間遺傳距離較大，蟲株之間存在一定的分化，且與地理來源關系密切： GD26和GD27均來自梅州市，這2株蟲株遺傳關系近； GD04、GD07和GD13均來自廣東省，它們遺傳關系相近； GD03和GD17均來自廣東省，之間遺傳距離小。JS蟲株和GD大部分蟲株之間遺傳分化明顯，但是GD09、GD12和JS蟲株遺傳距離近，與JS蟲株無明顯的遺傳差異(圖3)。

2.5 差異SNP位點挖掘及焦磷酸測序驗證

經過分析和過濾后共篩選到廣東省與江蘇省蟲株的差異SNP位點378個，差異位點主要富集功能與信傳導(GO:0023052)、應激反應(GO:0050896)、生物調控(GO： 0050789，GO： 0065007)等有關(圖4a)。選取廣東省與江蘇省蟲株的1個差異SNP位點(Contig007:3426052，A→G)，對測序分析結果進行試驗驗證。依據焦磷酸測序軟件生成待分析序列： TA/GTGCTAATAATTGAGGTGGGAGTAAG進行測序驗證。由焦磷酸測序結果可知，GD02、GD03、GD04、GD07、GD10、GD13、GD17、GD26、GD27和GD30在該位點上為G，而GD09、GD12和所有JS蟲株在該位點均為A(圖4b)。24株松材線蟲樣本經焦磷酸測序后的結果與待分析序列吻合，SNP位點的焦磷酸測序結果與重測序結果完全一致，且與聚類分析結果一致，驗證了相關位點作為鑒別不同地區松材線蟲的可行性。

3 討論

SNP分子標記技術在松材線蟲種群遺傳學方面的研究日趨成熟，該技術已用來研究松材線蟲種群遺傳差異性。Figueiredo等(2013)利用SNP分子標記技術分析來自5個國家的不同地理種群的松材線蟲，發現葡萄牙松材線蟲種群與中國和韓國的松材線蟲種群親緣關系更接近。黃金思等(2019)運用SNP標記分析中國廣東省松材線蟲蟲株，發現其具有不同的傳播來源，但尚未與其他地區的線蟲進行比較。

圖4 差異位點功能注釋及焦磷酸測序驗證

潘宏陽等(2009)的研究表明,松材線蟲病在中國有2個重要集中聚集區域，一個在江蘇和安徽，一個在廣東。筆者對廣東省蟲株和江蘇省蟲株綜合分析，除確定廣東省蟲株存在不同的傳播來源以外，還發現廣東省有部分蟲株和江蘇省蟲株存在相同的傳播來源。結合我國松材線蟲病發生時間， 1982年首次發現在江蘇南京，廣東省1988年首次發現，推測江蘇省和廣東省為我國2個松材線蟲病最早的傳播中心，且江蘇省的松材線蟲曾經傳播至廣東省。2個省份的松材線早蟲株聚類分析表明，江蘇省蟲株比較集中，各蟲株間的遺傳距離較近，而廣東省各蟲株間遺傳距離較遠。因此，江蘇省的松材線蟲蟲株之間可能存在交叉侵染的現象，基因流動比較頻繁，導致江蘇省蟲株遺傳差異性小。而廣東省蟲株中有2個與江蘇省蟲株具有相同的傳播來源，剩下的蟲株與其他所有樣本差異較大，說明該地區不僅有本土傳播的可能還存在國外傳入的外來種群。除此之外，江蘇和廣東省的年平均溫度差異較大，且在差異SNP富集中也發現了環境適應性相關功能。因此，松材線蟲的種群差異可能與外界環境因素如溫度等有密切關系。

歷史時期物種的遷移和擴散是決定種群遺傳結構的重要因素(黃族豪等， 2008)。一般認為種群受奠基者效應的影響，隨著種群的擴散，種群遺傳多樣性會逐漸的減少(Neietal.， 1975； Stone等， 2002)。特別是長距離的傳播擴散或者是環境差異大時，種群面臨的自然選擇壓力會更大，導致只有部分基因變異可能保存下來，因此新種群的遺傳多樣性會更低(黃族豪等， 2008)。

4 結論

比較江蘇省和廣東省松材線蟲種群的遺傳多樣性表明，大部分廣東省蟲株的遺傳多樣性要明顯的比江蘇省蟲株高。從江蘇省和廣東省遺傳多樣性的差異推測，廣東省部分蟲株擁有比江蘇省蟲株更特異的基因圖譜，可能更早發生松材線蟲病并可能是外源祖先的后代，而部分廣東省蟲株與江蘇省蟲株有相同的傳播來源。此外，經過驗證的相關差異位點可作為鑒別不同地區松材線蟲的新的候選SNP標記。為探究全國松材線蟲的傳播路線，后續研究應擴大取樣的地理范圍，進一步分析各個地理來源的種群遺傳結構，為檢疫和溯源工作的實施提供數據基礎。