基于SNP分子标记的华东地区松材线虫种群遗传分化研究

汪青桐，丁晓磊，叶建仁，史秀峰

(1.南京林业大学林学院，南方现代林业协同创新中心，江苏 南京 210037;2.上海市林业总站，上海 200072)

松材线虫病又称松树萎蔫病，是由病原松材线虫(Bursaphelenchusxylophilus)引起的松树毁灭性病害。松材线虫原产自北美洲，对美国本土松种并不造成严重危害[1-2]，但是对中国、日本和韩国等亚洲国家及葡萄牙和西班牙等欧洲国家的松林则造成巨大危害[3]。目前中国是受松材线虫病危害最严重的国家，自1982年在南京中山陵的黑松(Pinusthunbergii)上首次发现松材线虫病以来，至2021年已有19个省(自治区、直辖市)的728个县级行政区域成为松材线虫病疫区[4]。松材线虫作为一种入侵生物，在到达新环境后，由于遗传漂变以及环境的变化，会产生一定适应性变异，导致入侵种群在经历一个时期后可能会在遗传上出现相应的变化[5]。DNA分子标记技术是研究生物种群遗传变化的重要手段，它是以个体间核苷酸的差异作为标记基础，在DNA水平上对遗传物质变化的直接反应，能精准揭示生物种间或种内的遗传差异[6]。

SNP(single nucleotide polymorphism)分子标记技术是第3代DNA分子标记技术，其原理是染色体基因组中单个碱基的插入(insertion)、缺失(deletion)、转换(transition)、颠换(transversion)引起染色体基因组单个核苷酸突变，从而导致DNA序列的多态性[7]。SNP检测方法主要有酶切扩增多态性序列(CAPS)、单链构象多态性(SSCP)、直接测序和基因芯片等，直接测序是最直接可靠检测SNP的方法[8],现已在人类基因组学、医学、植物基因组学[9-12]等学科广泛应用。除人类外，对水稻(Oryzasativa)、高粱(Sorghumbicolor)、家犬(Canislupusfamiliars)、狼(Canislupus)等物种也已经建立了成熟的SNP数据库[13-15]。目前SNP标记技术在线虫上的应用主要集中在线虫抗性基因SNP位点的筛选和开发方面，研究较多的包括甜菜胞囊线虫(HeteroderaSchachtii)、大豆胞囊线虫(Heteroderaglycines)等[16-18]。2013年Figueiredo等[19]利用SNP标记分析了5个国家的松材线虫虫株的遗传结构，这是SNP标记技术在松材线虫这个物种上的首次应用。2015年完成了日本松材线虫全基因组的SNP分析，在JP基因组中发现了2 772 939个SNPs[20]。SNP分子标记技术在松材线虫种群遗传学方面的应用已经日趋成熟，成为研究松材线虫群体基因组学常用的方法之一。

本研究基于高通量测序方法，通过SNP分子标记技术，分析华东地区不同地理种群松材线虫虫株的SNPs信息，并对其进行聚类分析，研究其SNP遗传多样性、遗传分化和基因流动，为建立我国松材线虫病的疫源追溯体系提供必要的基础信息。

1.1 虫株收集及鉴定纯化

采集华东地区安徽、福建、江苏、江西、山东和浙江6个省份松材线虫病疫木样本，采用贝尔曼漏斗法分离松材线虫。分离后获得的线虫，参照谢辉[21]的松材线虫分类系统，在显微镜下观察，对松材线虫初步进行形态学鉴定后，再应用SCAR标记[22]方法对各虫株进行松材线虫进行分子鉴定。

鉴定为松材线虫后，挑取分离虫液中的松材线虫雌雄成虫各15条左右，接入长好灰葡萄孢(Botrytiscinerea)的potato-dextrose-agar(PDA)平板中，放置于28 ℃培养箱，至松材线虫取食完毕后，使用贝尔曼漏斗法分离线虫，重复进行1次形态学鉴定，确定为松材线虫后，使用0.05%质量分数硫酸链霉素清洗，再用无菌水洗涤3遍，放入冰箱4 ℃条件下保存备用。活体虫株均保存于南京林业大学森林病理实验室松材线虫虫株资源库。

经纯化培养，共获得67株松材线虫样本，来自华东地区60个县级行政区，其中：安徽(AH)11株、福建(FJ)11株、江苏(JS)18株、江西(JX)6株、山东(SD)9株以及浙江(ZJ)12株。67个松材线虫虫株的样本信息及采样地松材线虫病首次发生时间(入侵年份)见表1。

表1 松材线虫67个虫株样本信息

1.2 供试线虫DNA提取、保存及线虫全基因组测序

参照陈凤毛[23]的方法提取松材线虫DNA。提取的DNA保存1份于南京林业大学森林病理实验室松材线虫虫株DNA资源库。

将合格的松材线虫DNA送到武汉未来组生物公司，进行Illumina高通量基因组测序，测序平台为HiSeq 4000。基因组重测序方法为非链特异性、150 bp双末端测序。

1.3 数据处理

1.3.1 数据及SNP位点信息统计

1)质量控制：通过FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)软件，检查数据的基本信息、碱基的质量值、每条reads序列的质量值、每条序列的ATCG组成以及长度分布、序列中duplication程度等基本信息。

2)去除测序接头污染序列：用Cutadapt软件去除测序接头(adapters)。

3)基因组比对：通过BWA和SAMtools(http://samtools.sourceforge.net/samtools.shtml)软件将测序结果与松材线虫全基因组进行比对，获得用于后续分析的基因组位置信息。

4)PCR重复序列去除：用Picard(http://broadinstitute.github.io/picard/)软件去除由PCR扩增而形成的重复片段。

5)SNP位点鉴定：使用Freebayes，参数设置为“-u-C 5-e 50--standard-filters--min-coverage 10”，其他参数均为软件默认参数。

6)样本基因型分析与统计：主要由VCFtools(http://vcftools.sourceforge.net/)完成初步统计，之后利用个性化脚本进行深度信息挖掘。

1.3.2 松材线虫种群遗传分化分析

对松材线虫测序虫株的SNP位点信息进行生物软件分析，应用R语言包中的SNPRealte软件对等位基因频率较低、连锁不平衡(linkage disequilibrium)较高的SNP位点进行过滤，同时绘制PCA主成分分析图。

使用Treemix软件构建系统发育树[24]。由于Treemix是基于等位基因频率作为遗传距离计算依据，所以预先使用Plink[25]对所有SNP位点进行等位基因频率计算，之后使用参数为“-noss-global”进行最大似然树构建，“-m”添加可能的迁移事件，并在最大似然树上进行标记。

2.1 华东区地区松材线虫虫株的全基因组测序

应用Illumina HiSeq 4000平台对所采集的华东地区67个松材线虫虫株进行重测序，获得的测序数据量539.95 G，虫株获得的reads数为39 032 238～92 873 782，测序数据Q30为87.0%～94.4%，与参考基因组(AH1)的比对率在94.5%左右，覆盖深度大于40×，具体各虫株测序结果见表2。

表2 67株松材线虫虫株测序质量

2.2 华东地区松材线虫SNP基因型统计

利用生物信息学软件统计SNP基因型发现，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种基因型出现的频率明显高于其他6种基因型(图1)。

2.3 华东地区松材线虫SNP位点分析

对67个松材线虫虫株的全基因组中的SNP位点信息进行统计分析，结果见图2。在67个虫株的基因组中，共获得6 531 684个SNP位点，平均每个虫株97 488个SNP，不同虫株的SNP位点数量差距很大，SNP总量最多的是ZJ01(827 950)，来自浙江省宁波市镇海区，最少的是JS09(29 015)，来自江苏省淮安市盱眙县；
SNP纯合子数量最多的也是ZJ01(710 030)，最少的也是JS09(12 089)；
缺失SNP数量最多为JS09(1 416 861)，最少为SD06(339 209)，来自山东省烟台市芝罘区；
特有SNP数量最多的是ZJ01(18 385)，最少的是SD01(1)，来自山东省威海市荣成市。由图2可以看出，华东地区6个省份内不同虫株间的SNP总量、纯合子数量、缺失SNP数量以及特有SNP数量均存在明显的差异。但是显著性差异分析表明，SNP总量、纯合子数量、缺失SNP数量以及特有SNP数量在各省份之间均无显著性差异(df=5，F=0.5191，P>0.05)。说明虽然不同虫株间的SNP位点数量和种类差异很大，它们间的差异与省际行政区划没有呈现特别明显的相关性。

将华东地区松材线虫病入侵年份划分为4个阶段：1982—1990年、1991—2000年、2001—2010年、2011—2016年，其对应SNP类型及数量分布见图3。通过显著性差异分析发现，4个阶段之间在SNP总量、纯合子数量、缺失SNP数量以及特有SNP数量上，均无差异性显著(df=3，F=0.521 9，P>0.05)，说明华东地区松材线虫虫株的SNP多样性变化与入侵时间也没有呈现明显规律性。

2.4 华东地区松材线虫的主成分分析

对来自6个省份的67株虫株的SNP位点信息和基因型特征进行主成分分析(PCA)，结果(图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。结合华东地区6省3个类群虫株的地理来源发现，华东地区6省的松材线虫虫株大多归于类群1，3个类群与地理区域间存在着一定的联系，浙江和江西的虫株都属于类群1，类群2和类群3中只有其他4个省的少量虫株。

2.5 华东地区松材线虫的基因交流情况

根据PCA结果，将3个类群中相同省份的虫株划分为1个群体(如1-AH包含类群1中的所有安徽虫株)，使用Treemix绘制最大似然树，标记基因迁移路径(图5)。由于Treemix是一种混群检测种群间基因流动的手段，类群2中山东省虫株以及类群3中安徽省虫株样本量过少，所以未加入分析，其余松材线虫虫株划分为10个群体。从绘制的最大似然树中看出，类群2和类群3中省份间遗传分化小，类群1在各省份存在一定的遗传分化现象，其中浙江和江西遗传距离较近。通过Treemix检测，得到了两条可能的基因迁移路径分别为2-JS向1-AH、3-SD向2-AH迁移，表明1-AH和2-AH区域存在被其他区域其他类群虫株再次入侵的情况。

我国松材线虫病主要存在两个扩散中心，即广东省和江苏省[26]，呈现出两个显著的聚集分布区，一个位于广东境内，另一个包括江苏、安徽大部和浙江西北部[27]。苏皖聚集分布区域的形成是以南京为中心，由早期疫点向外内源性扩散所致，苏皖聚集分布区各疫点松材线虫具有较高的同源性，很可能来自同一个疫源[28]。本研究发现，华东地区虫株的主要SNP基因型均相同，并且类群1包含了华东地区6个省份大部分的虫株，结合聚类结果和各地病害发生时间推测：大部分华东地区的松材线虫种群间具有相同的传播来源，并且以苏皖为中心向外扩散。黄金思等[29]曾对广东省松材线虫虫株进行SNP标记分析，将广东省松材线虫虫株分为a、b、c等3个类群，且认为a和b、c类群具有不同的传播来源。本研究结果与广东省松材线虫虫株的SNP统计结果对比发现，华东地区虫株的SNP数量明显低于广东省的b、c类群虫株，主要基因型种类也有明显差异，而与广东省a类群虫株相比，在SNP数量和主要基因型种类上均无明显差异，由此推测华东地区的松材线虫与广东省部分虫株亲缘关系相近，虫株间可能存在相同的传播来源。

大部分物种在入侵期间通常会由于环境胁迫，发生自然选择，产生奠基者效应和遗传漂变，导致物种遗传多样性的减少和丧失[30-31]。一些报道证实松材线虫也会由于奠基者效应和遗传漂变，发生遗传多样性的丧失，导致种群间出现遗传分化，且与地理来源有明显的相关性[32-34]。Cheng等[35]使用AFLP标记分析我国不同地区松材线虫遗传多样性发现，广东省虫株多样性最高，安徽、江苏、贵州和湖北地区虫株遗传多样性要高于浙江沿海地区，存在明显的地理来源差异。本研究中，华东地区松材线虫虫株的SNP多样性和3个类群虫株间的遗传分化，没有表现出按区域划分的现象。这主要是因为松材线虫依赖于人为因素传播呈现出跳跃式发展的特点[28]，而华东地区各省份间松材线虫可能存在多条传播路径。Cheng等[35]也分析了我国不同入侵时期的松材线虫种群的遗传多样性，发现其遗传多样性在不同入侵阶段没有明显的变化，推测原因为未发生明显的奠基者效应和遗传漂变。本研究结果显示，华东地区不同入侵阶段的松材线虫遗传多样性也没有显著差异。上述结果可能存在的原因：一方面，华东地区的气候、寄主植物、媒介昆虫等都较为相似且适宜松材线虫生存和繁殖，大多虫株到达新环境后没有受到或者受到较小的环境胁迫，能够快速定殖，所以没有出现遗传漂变导致遗传多样性降低的现象；
另一方面，类群间基因流结果表明，华东地区的松材线虫种群间可能存在多次入侵的现象，不同地理来源的种群之间基因交流频繁，弥补了遗传漂变带来的影响，并且未发生奠基者效应。针对华东6省松材线虫种群间可能存在多次入侵的现象，应该进一步加大检疫力度，控制疫木的流通。