基于SNP分子標記的華東地區(qū)松材線蟲種群遺傳分化研究

汪青桐，丁曉磊，葉建仁，史秀峰

(1.南京林業(yè)大學(xué)林學(xué)院，南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心，江蘇 南京 210037;2.上海市林業(yè)總站，上海 200072)

松材線蟲病又稱松樹萎蔫病，是由病原松材線蟲(Bursaphelenchusxylophilus)引起的松樹毀滅性病害。松材線蟲原產(chǎn)自北美洲，對美國本土松種并不造成嚴重危害[1-2]，但是對中國、日本和韓國等亞洲國家及葡萄牙和西班牙等歐洲國家的松林則造成巨大危害[3]。目前中國是受松材線蟲病危害最嚴重的國家，自1982年在南京中山陵的黑松(Pinusthunbergii)上首次發(fā)現(xiàn)松材線蟲病以來，至2021年已有19個省(自治區(qū)、直轄市)的728個縣級行政區(qū)域成為松材線蟲病疫區(qū)[4]。松材線蟲作為一種入侵生物，在到達新環(huán)境后，由于遺傳漂變以及環(huán)境的變化，會產(chǎn)生一定適應(yīng)性變異，導(dǎo)致入侵種群在經(jīng)歷一個時期后可能會在遺傳上出現(xiàn)相應(yīng)的變化[5]。DNA分子標記技術(shù)是研究生物種群遺傳變化的重要手段，它是以個體間核苷酸的差異作為標記基礎(chǔ)，在DNA水平上對遺傳物質(zhì)變化的直接反應(yīng)，能精準揭示生物種間或種內(nèi)的遺傳差異[6]。

SNP(single nucleotide polymorphism)分子標記技術(shù)是第3代DNA分子標記技術(shù)，其原理是染色體基因組中單個堿基的插入(insertion)、缺失(deletion)、轉(zhuǎn)換(transition)、顛換(transversion)引起染色體基因組單個核苷酸突變，從而導(dǎo)致DNA序列的多態(tài)性[7]。SNP檢測方法主要有酶切擴增多態(tài)性序列(CAPS)、單鏈構(gòu)象多態(tài)性(SSCP)、直接測序和基因芯片等，直接測序是最直接可靠檢測SNP的方法[8],現(xiàn)已在人類基因組學(xué)、醫(yī)學(xué)、植物基因組學(xué)[9-12]等學(xué)科廣泛應(yīng)用。除人類外，對水稻(Oryzasativa)、高粱(Sorghumbicolor)、家犬(Canislupusfamiliars)、狼(Canislupus)等物種也已經(jīng)建立了成熟的SNP數(shù)據(jù)庫[13-15]。目前SNP標記技術(shù)在線蟲上的應(yīng)用主要集中在線蟲抗性基因SNP位點的篩選和開發(fā)方面，研究較多的包括甜菜胞囊線蟲(HeteroderaSchachtii)、大豆胞囊線蟲(Heteroderaglycines)等[16-18]。2013年Figueiredo等[19]利用SNP標記分析了5個國家的松材線蟲蟲株的遺傳結(jié)構(gòu)，這是SNP標記技術(shù)在松材線蟲這個物種上的首次應(yīng)用。2015年完成了日本松材線蟲全基因組的SNP分析，在JP基因組中發(fā)現(xiàn)了2 772 939個SNPs[20]。SNP分子標記技術(shù)在松材線蟲種群遺傳學(xué)方面的應(yīng)用已經(jīng)日趨成熟，成為研究松材線蟲群體基因組學(xué)常用的方法之一。

本研究基于高通量測序方法，通過SNP分子標記技術(shù)，分析華東地區(qū)不同地理種群松材線蟲蟲株的SNPs信息，并對其進行聚類分析，研究其SNP遺傳多樣性、遺傳分化和基因流動，為建立我國松材線蟲病的疫源追溯體系提供必要的基礎(chǔ)信息。

1 材料與方法

1.1 蟲株收集及鑒定純化

采集華東地區(qū)安徽、福建、江蘇、江西、山東和浙江6個省份松材線蟲病疫木樣本，采用貝爾曼漏斗法分離松材線蟲。分離后獲得的線蟲，參照謝輝[21]的松材線蟲分類系統(tǒng)，在顯微鏡下觀察，對松材線蟲初步進行形態(tài)學(xué)鑒定后，再應(yīng)用SCAR標記[22]方法對各蟲株進行松材線蟲進行分子鑒定。

鑒定為松材線蟲后，挑取分離蟲液中的松材線蟲雌雄成蟲各15條左右，接入長好灰葡萄孢(Botrytiscinerea)的potato-dextrose-agar(PDA)平板中，放置于28 ℃培養(yǎng)箱，至松材線蟲取食完畢后，使用貝爾曼漏斗法分離線蟲，重復(fù)進行1次形態(tài)學(xué)鑒定，確定為松材線蟲后，使用0.05%質(zhì)量分數(shù)硫酸鏈霉素清洗，再用無菌水洗滌3遍，放入冰箱4 ℃條件下保存?zhèn)溆谩；铙w蟲株均保存于南京林業(yè)大學(xué)森林病理實驗室松材線蟲蟲株資源庫。

經(jīng)純化培養(yǎng)，共獲得67株松材線蟲樣本，來自華東地區(qū)60個縣級行政區(qū)，其中：安徽(AH)11株、福建(FJ)11株、江蘇(JS)18株、江西(JX)6株、山東(SD)9株以及浙江(ZJ)12株。67個松材線蟲蟲株的樣本信息及采樣地松材線蟲病首次發(fā)生時間(入侵年份)見表1。

表1 松材線蟲67個蟲株樣本信息

1.2 供試線蟲DNA提取、保存及線蟲全基因組測序

參照陳鳳毛[23]的方法提取松材線蟲DNA。提取的DNA保存1份于南京林業(yè)大學(xué)森林病理實驗室松材線蟲蟲株DNA資源庫。

將合格的松材線蟲DNA送到武漢未來組生物公司，進行Illumina高通量基因組測序，測序平臺為HiSeq 4000。基因組重測序方法為非鏈特異性、150 bp雙末端測序。

1.3 數(shù)據(jù)處理

1.3.1 數(shù)據(jù)及SNP位點信息統(tǒng)計

1)質(zhì)量控制：通過FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)軟件，檢查數(shù)據(jù)的基本信息、堿基的質(zhì)量值、每條reads序列的質(zhì)量值、每條序列的ATCG組成以及長度分布、序列中duplication程度等基本信息。

2)去除測序接頭污染序列：用Cutadapt軟件去除測序接頭(adapters)。

3)基因組比對：通過BWA和SAMtools(http://samtools.sourceforge.net/samtools.shtml)軟件將測序結(jié)果與松材線蟲全基因組進行比對，獲得用于后續(xù)分析的基因組位置信息。

4)PCR重復(fù)序列去除：用Picard(http://broadinstitute.github.io/picard/)軟件去除由PCR擴增而形成的重復(fù)片段。

5)SNP位點鑒定：使用Freebayes，參數(shù)設(shè)置為“-u-C 5-e 50--standard-filters--min-coverage 10”，其他參數(shù)均為軟件默認參數(shù)。

6)樣本基因型分析與統(tǒng)計：主要由VCFtools(http://vcftools.sourceforge.net/)完成初步統(tǒng)計，之后利用個性化腳本進行深度信息挖掘。

1.3.2 松材線蟲種群遺傳分化分析

對松材線蟲測序蟲株的SNP位點信息進行生物軟件分析，應(yīng)用R語言包中的SNPRealte軟件對等位基因頻率較低、連鎖不平衡(linkage disequilibrium)較高的SNP位點進行過濾，同時繪制PCA主成分分析圖。

使用Treemix軟件構(gòu)建系統(tǒng)發(fā)育樹[24]。由于Treemix是基于等位基因頻率作為遺傳距離計算依據(jù)，所以預(yù)先使用Plink[25]對所有SNP位點進行等位基因頻率計算，之后使用參數(shù)為“-noss-global”進行最大似然樹構(gòu)建，“-m”添加可能的遷移事件，并在最大似然樹上進行標記。

2 結(jié)果與分析

2.1 華東區(qū)地區(qū)松材線蟲蟲株的全基因組測序

應(yīng)用Illumina HiSeq 4000平臺對所采集的華東地區(qū)67個松材線蟲蟲株進行重測序，獲得的測序數(shù)據(jù)量539.95 G，蟲株獲得的reads數(shù)為39 032 238～92 873 782，測序數(shù)據(jù)Q30為87.0%～94.4%，與參考基因組(AH1)的比對率在94.5%左右，覆蓋深度大于40×，具體各蟲株測序結(jié)果見表2。

表2 67株松材線蟲蟲株測序質(zhì)量

2.2 華東地區(qū)松材線蟲SNP基因型統(tǒng)計

利用生物信息學(xué)軟件統(tǒng)計SNP基因型發(fā)現(xiàn)，67個松材線蟲蟲株中SNP基因型種類共有12種，分別為：A→C、A→G、A→T、C→A、C→G、C→T、G→A、G→C、G→T、T→A、T→C和T→G，所有蟲株均有以上基因型。所有樣本A→G、C→G、C→T、G→A、G→C、T→C這6種基因型出現(xiàn)的頻率明顯高于其他6種基因型(圖1)。

2.3 華東地區(qū)松材線蟲SNP位點分析

對67個松材線蟲蟲株的全基因組中的SNP位點信息進行統(tǒng)計分析，結(jié)果見圖2。在67個蟲株的基因組中，共獲得6 531 684個SNP位點，平均每個蟲株97 488個SNP，不同蟲株的SNP位點數(shù)量差距很大，SNP總量最多的是ZJ01(827 950)，來自浙江省寧波市鎮(zhèn)海區(qū)，最少的是JS09(29 015)，來自江蘇省淮安市盱眙縣；SNP純合子數(shù)量最多的也是ZJ01(710 030)，最少的也是JS09(12 089)；缺失SNP數(shù)量最多為JS09(1 416 861)，最少為SD06(339 209)，來自山東省煙臺市芝罘區(qū)；特有SNP數(shù)量最多的是ZJ01(18 385)，最少的是SD01(1)，來自山東省威海市榮成市。由圖2可以看出，華東地區(qū)6個省份內(nèi)不同蟲株間的SNP總量、純合子數(shù)量、缺失SNP數(shù)量以及特有SNP數(shù)量均存在明顯的差異。但是顯著性差異分析表明，SNP總量、純合子數(shù)量、缺失SNP數(shù)量以及特有SNP數(shù)量在各省份之間均無顯著性差異(df=5，F(xiàn)=0.5191，P>0.05)。說明雖然不同蟲株間的SNP位點數(shù)量和種類差異很大，它們間的差異與省際行政區(qū)劃沒有呈現(xiàn)特別明顯的相關(guān)性。

將華東地區(qū)松材線蟲病入侵年份劃分為4個階段：1982—1990年、1991—2000年、2001—2010年、2011—2016年，其對應(yīng)SNP類型及數(shù)量分布見圖3。通過顯著性差異分析發(fā)現(xiàn)，4個階段之間在SNP總量、純合子數(shù)量、缺失SNP數(shù)量以及特有SNP數(shù)量上，均無差異性顯著(df=3，F(xiàn)=0.521 9，P>0.05)，說明華東地區(qū)松材線蟲蟲株的SNP多樣性變化與入侵時間也沒有呈現(xiàn)明顯規(guī)律性。

2.4 華東地區(qū)松材線蟲的主成分分析

對來自6個省份的67株蟲株的SNP位點信息和基因型特征進行主成分分析(PCA)，結(jié)果(圖4)顯示67個蟲株被分為明顯的3個類群：類群1最大，包括了所有浙江蟲株和江西蟲株，以及江蘇、安徽、山東和福建的大部分蟲株，共有46個蟲株；類群2涉及江蘇、安徽、山東和福建的部分蟲株，包括14個蟲株，即JS02、JS05、JS06、JS08、JS12、JS13、JS14、JS18、AH02、AH03、AH05、SD01、SD03和FJ01；類群3僅涉及安徽、福建和山東的7個蟲株，即AH14、FJ08、FJ11、FJ12、SD06、SD07和SD08。結(jié)合華東地區(qū)6省3個類群蟲株的地理來源發(fā)現(xiàn)，華東地區(qū)6省的松材線蟲蟲株大多歸于類群1，3個類群與地理區(qū)域間存在著一定的聯(lián)系，浙江和江西的蟲株都屬于類群1，類群2和類群3中只有其他4個省的少量蟲株。

2.5 華東地區(qū)松材線蟲的基因交流情況

根據(jù)PCA結(jié)果，將3個類群中相同省份的蟲株劃分為1個群體(如1-AH包含類群1中的所有安徽蟲株)，使用Treemix繪制最大似然樹，標記基因遷移路徑(圖5)。由于Treemix是一種混群檢測種群間基因流動的手段，類群2中山東省蟲株以及類群3中安徽省蟲株樣本量過少，所以未加入分析，其余松材線蟲蟲株劃分為10個群體。從繪制的最大似然樹中看出，類群2和類群3中省份間遺傳分化小，類群1在各省份存在一定的遺傳分化現(xiàn)象，其中浙江和江西遺傳距離較近。通過Treemix檢測，得到了兩條可能的基因遷移路徑分別為2-JS向1-AH、3-SD向2-AH遷移，表明1-AH和2-AH區(qū)域存在被其他區(qū)域其他類群蟲株再次入侵的情況。

3 討 論

我國松材線蟲病主要存在兩個擴散中心，即廣東省和江蘇省[26]，呈現(xiàn)出兩個顯著的聚集分布區(qū)，一個位于廣東境內(nèi)，另一個包括江蘇、安徽大部和浙江西北部[27]。蘇皖聚集分布區(qū)域的形成是以南京為中心，由早期疫點向外內(nèi)源性擴散所致，蘇皖聚集分布區(qū)各疫點松材線蟲具有較高的同源性，很可能來自同一個疫源[28]。本研究發(fā)現(xiàn)，華東地區(qū)蟲株的主要SNP基因型均相同，并且類群1包含了華東地區(qū)6個省份大部分的蟲株，結(jié)合聚類結(jié)果和各地病害發(fā)生時間推測：大部分華東地區(qū)的松材線蟲種群間具有相同的傳播來源，并且以蘇皖為中心向外擴散。黃金思等[29]曾對廣東省松材線蟲蟲株進行SNP標記分析，將廣東省松材線蟲蟲株分為a、b、c等3個類群，且認為a和b、c類群具有不同的傳播來源。本研究結(jié)果與廣東省松材線蟲蟲株的SNP統(tǒng)計結(jié)果對比發(fā)現(xiàn)，華東地區(qū)蟲株的SNP數(shù)量明顯低于廣東省的b、c類群蟲株，主要基因型種類也有明顯差異，而與廣東省a類群蟲株相比，在SNP數(shù)量和主要基因型種類上均無明顯差異，由此推測華東地區(qū)的松材線蟲與廣東省部分蟲株親緣關(guān)系相近，蟲株間可能存在相同的傳播來源。

大部分物種在入侵期間通常會由于環(huán)境脅迫，發(fā)生自然選擇，產(chǎn)生奠基者效應(yīng)和遺傳漂變，導(dǎo)致物種遺傳多樣性的減少和喪失[30-31]。一些報道證實松材線蟲也會由于奠基者效應(yīng)和遺傳漂變，發(fā)生遺傳多樣性的喪失，導(dǎo)致種群間出現(xiàn)遺傳分化，且與地理來源有明顯的相關(guān)性[32-34]。Cheng等[35]使用AFLP標記分析我國不同地區(qū)松材線蟲遺傳多樣性發(fā)現(xiàn)，廣東省蟲株多樣性最高，安徽、江蘇、貴州和湖北地區(qū)蟲株遺傳多樣性要高于浙江沿海地區(qū)，存在明顯的地理來源差異。本研究中，華東地區(qū)松材線蟲蟲株的SNP多樣性和3個類群蟲株間的遺傳分化，沒有表現(xiàn)出按區(qū)域劃分的現(xiàn)象。這主要是因為松材線蟲依賴于人為因素傳播呈現(xiàn)出跳躍式發(fā)展的特點[28]，而華東地區(qū)各省份間松材線蟲可能存在多條傳播路徑。Cheng等[35]也分析了我國不同入侵時期的松材線蟲種群的遺傳多樣性，發(fā)現(xiàn)其遺傳多樣性在不同入侵階段沒有明顯的變化，推測原因為未發(fā)生明顯的奠基者效應(yīng)和遺傳漂變。本研究結(jié)果顯示，華東地區(qū)不同入侵階段的松材線蟲遺傳多樣性也沒有顯著差異。上述結(jié)果可能存在的原因：一方面，華東地區(qū)的氣候、寄主植物、媒介昆蟲等都較為相似且適宜松材線蟲生存和繁殖，大多蟲株到達新環(huán)境后沒有受到或者受到較小的環(huán)境脅迫，能夠快速定殖，所以沒有出現(xiàn)遺傳漂變導(dǎo)致遺傳多樣性降低的現(xiàn)象；另一方面，類群間基因流結(jié)果表明，華東地區(qū)的松材線蟲種群間可能存在多次入侵的現(xiàn)象，不同地理來源的種群之間基因交流頻繁，彌補了遺傳漂變帶來的影響，并且未發(fā)生奠基者效應(yīng)。針對華東6省松材線蟲種群間可能存在多次入侵的現(xiàn)象，應(yīng)該進一步加大檢疫力度，控制疫木的流通。