广东省与江苏省松材线虫种群遗传结构差异分析*

丁晓磊 汪青桐 林司曦 赵瑞文 张 悦 叶建仁

(南方现代林业协同创新中心 南京林业大学林学院 南京 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℃冰箱保存备用。活体线虫虫株均保存于南京林业大学森林病理实验室松材线虫虫株资源库。

表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。总体来看，广东省虫株的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标记。为探究全国松材线虫的传播路线，后续研究应扩大取样的地理范围，进一步分析各个地理来源的种群遗传结构，为检疫和溯源工作的实施提供数据基础。