基于COⅠ基因的我国悬铃木方翅网蝽(半翅目： 网蝽科) 种群遗传多样性和遗传结构分析*

张元臣 卢绍辉 龚东风 马兴莉

(1.安阳工学院生物与食品工程学院 安阳 455000； 2.河南省太行山林业有害生物野外科学观测研究站 林州 456550；3.河南省林业科学研究院 郑州 450002； 4.北卡罗来纳州立大学农业与生命科学学院有害生物综合治理中心 罗利27695)

物种遗传变异与其适应环境的能力和进化程度密切相关，也是物种在新环境下能够定殖和迁移的关键因素，而生物入侵能够引起种群遗传变异，这为研究入侵生物的生态和进化问题提供了理想的模式(Amourouxetal., 2013; Lau, 2015)。分子标记是研究物种遗传变异非常有效的工具，能够为物种的起源和分化提供直接证据，可揭示入侵物种的扩散路径和结果，并推测未来的种群发生情况(李菁等, 2010; 宾淑英等, 2014; 杨顺义等, 2018)。目前已经开发出了多种分子标记技术，而线粒体基因具有进化迅速、垂直母系遗传、重组困难等特点，使其广泛应用于昆虫种群的遗传变异和遗传结构的研究中(王克勤等， 2018; Cameron, 2014; Menetal., 2017)。很多研究已经利用昆虫线粒体细胞色素氧化酶亚基Ⅰ(cytochrome C oxidase subunit Ⅰ，COⅠ)基因对种群遗传分化、入侵路径和入侵来源等进行了深入研究(王蒙等, 2014; 王戎疆等， 2018; Seraphimetal., 2016; Ariasetal., 2018; Chinetal., 2018; Aslametal., 2019)。

悬铃木方翅网蝽(Corythuchaciliata)(半翅目Hemiptera: 网蝽科Tingidae)(以下简称网蝽)具有快速适应新栖息地的特点，在新的生境里经常大面积暴发成灾(卢绍辉等, 2018; Maceljskietal., 1974)； 目前已经成为悬铃木属(Platanus)树木上的首要害虫(Bures, 1997)，防治非常困难(Cravedi, 2000)。2002年在长沙首次发现该网蝽(Streito, 2006)， 2006年开始在武汉市区大面积暴发(李传仁等, 2007)，此后，该网蝽在长江中下游的多个区域相续被发现，严重为害悬铃木属植物，给城市绿化带来了非常大的危害(Juetal., 2009)。自该网蝽入侵我国以来，昆虫研究者对其生态生物学特性、防治方法、重要基因的功能等进行了深入的研究(王福莲等, 2008; 郝德君等, 2012; 付宁宁等, 2017; 李峰奇等, 2018; Juetal., 2011)。Yang 等(2014)发掘了9个该网蝽微卫星标记位点，并对这9个微卫星标记位点在不同地理种群中的有效性进行了验证。进一步利用这些微卫星标记位点和3个线粒体基因的部分核苷酸序列，对收集于2009年到2014年间我国的21个网蝽地理种群和国外的2个网蝽地理种群进行了遗传结构研究，结果发现这些种群间具有高度的遗传分化，并推测该网蝽从东海岸入侵我国，之后向全国范围进行扩散(Yangetal., 2017)。然而，随着时间的推移，悬铃木方翅网蝽在我国入侵范围急剧扩大，因此有必要对该网蝽遗传结构进行分析和扩散路径的重建。基于此，本研究采用线粒体COⅠ标记基因对采集于2015年我国12个不同地理种群网蝽的遗传多样性与遗传结构进行分析，旨在探讨网蝽在中国不同地理种群间的遗传变异情况、分化程度及基因流水平，据此分析遗传分化产生的原因及在我国扩散的中心区域，并进行了入侵趋势的分析，从而为有效防治网蝽提供基础理论依据。

1 材料与方法

1.1 试验材料的采集

选择网蝽为害的悬铃木，随机取下一片叶片，用毛笔将叶片上的网蝽成虫挑到含有100%酒精的2 mL离心管中，带回室内用显微镜进行确认，之后将确定的网蝽成虫样品放到-20 ℃冰箱中保存。每个点随机选取间距在200 m以上的悬铃木20棵，每一棵树收集1头成虫，每个收集点共20头网蝽成虫被收集。网蝽成虫收集时间在2015年7月，收集地分别在长沙、武汉、上海、北京、石家庄、济南、太原、杭州、成都、重庆、贵阳和郑州共12个城市。

1.2 基因组DNA的提取

将100%酒精保存的网蝽成虫从离心管中挑出来，放到10 cm培养皿中，用ddH2O冲洗2次。将网蝽放到吸水纸上将水吸干，按照DNA提取试剂盒(DP304-03)(天根生化科技公司)提取DNA。用1.2%～1.5%的琼脂糖电泳和NanoDrop(赛默飞世尔科技公司)检测基因组DNA的质量和浓度后，将其稀释同样的浓度(30 μg·μL-1)，放置于-20 ℃条件下保存备用。

1.3 COⅠ基因扩增、测序及拼接

根据Yang等(2017)的方法对COⅠ基因进行扩增，将PCR产物送往上海生工进行双向测序。将测序获得的核苷酸序列去除掉引物，剩余核苷酸序列放到NCBI数据库中进行BLASTn比对，根据比对结果验证网蝽COⅠ基因核苷酸序列的正确性，随后以网蝽COⅠ核苷酸序列为基础进行后续的数据分析。

1.4 数据分析

参考呼晓庆等(2019)的方法对不同地理种群网蝽线粒体COⅠ序列进行核苷酸组分、保守位点、多态性位点、单倍型组成、单倍型数量、单倍型分布(Hap)、单倍型多样性(haplotype diversity,Hd)、核苷酸多样性(nucleotide diversity,π)、核苷酸平均差异数(average number of nucleotide differences,k)、不同地理种群间的固定系数FST，并依据(1-FST)/2FST公式计算种群间基因流动值(Nm)。根据Qin等(2016)的方法构建了9个网蝽单倍型网状进化图，并进行网蝽不同种群间的遗传距离与地理距离的相关性分析。为了方便对网蝽在主要分布区的遗传多样性进行分析，本文根据不同地理种群间地理位置将主要分布区分为华中组群(长沙、武汉)、华东组群(上海、杭州)、华北组群(北京、石家庄、太原、郑州和济南)、西南组群(成都、重庆和贵阳)。对4个组群进行分子方差分析(analysis of molecular variance,AMOVA)，并计算不同组群网蝽间的FST值和Nm值，同时对不同组群间进行种群历史动态分析得到扩张前后网蝽有效种群的大小(region mutation parameter,θ)、网蝽种群的扩张时间(T)、Tajima’sD值、Fu’sFS值及期望错配分布值与观察错配分布值之间的平方和(SSD)。通过MIGRATE V.3.2.1软件(Beerli, 2006)计算网蝽组群有效迁移率(migration rate,M)； 利用MEGA v7.0.26软件根据本研究中12个网蝽种群和NCBI上其他国内外网蝽种群的COⅠ基因核苷酸序列，基于Kimura 2-Parameter模型计算距离矩阵并根据邻接法建立系统发育树。

2 结果与分析

2.1 网蝽COⅠ基因核苷酸的碱基构成及序列

本研究共获得了240条761 bp的网蝽COⅠ基因序列，A、T、G和C核苷酸含量分别为35.0%、32.6%、17.5%和14.9%。A+T核苷酸含量(67.6%)明显高于G+C核苷酸含量(32.4%)，核苷酸组成类型符合昆虫线粒体核苷酸的特征。240个COⅠ序列含有13个多态性位点(polymorphic sites)，为总碱基数的1.71%。发现的多态性位点全部是核苷酸间的相互替换引起的，含有4个自裔位点和9个简约信息位点。

2.2 网蝽地理种群的单倍型和核苷酸多样性

在240条网蝽COⅠ基因核苷酸中发现了9个单倍型，分别有4个共享单倍型和5个独享单倍型(表1)，表明网蝽地理种群间已经出现了遗传分化。所有单倍型中分布最多的是Hap1和Hap3单倍型，都被9个种群共享，而Hap4和Hap7单倍型分别被3个和2个种群共享，其他单倍型为独享单倍型。12个地理种群单倍型在1～6个之间，其中GY种群拥有最高的单倍型数量(6个)，而SH和CS种群拥有最低的单倍型数量(1个)(表1)。单倍型网络图表明，Hap1单倍型居于网络图的中心，与其相关联的单倍型数量最多(图1)。Hap3单倍型共享个体数量也较多，然而与其相关联的单倍型数量比较少，而且Hap3单倍型是从Hap7单倍型分化出来的，Hap7单倍型直接从Hap1单倍型分化而来。基于此，本研究认为Hap1为祖先单倍型。

表1 网蝽COⅠ基因单倍型的分布①Tab.1 Distribution information of haplotype from COⅠ gene in C. ciliata

图1 网蝽单倍型网络进化Fig. 1 Haplotype Median-Joining network using COⅠ gene from C. ciliata

网蝽总种群单倍型多样度Hd值为0.584，核苷酸多样度π值为0.003，核苷酸平均差异数k值为2.276(表2)。不同地理种群间Hd值在0.000～0.758之间；π值在0.000～0.005之间；k值在0.000～3.853之间(表2)。在分析的12个地理种群中，GY种群拥有最高的单倍型遗传多样性(Hd值为0.758)，而CS和SH种群拥有最低的单倍型多样性(Hd值均为0)(表2)。总体上，12个网蝽地理种群中遗传多样性差异比较明显，CS、SH种群无多态性，12个种群遗传多样性从中国南部——贵阳到华中地区——郑州等地区逐渐减少，而从华中到华北——北京等地区遗传多样性逐渐增大。

表2 12个网蝽地理种群COⅠ 基因单倍型多样度及核苷酸多样性①Tab.2 Haplotype diversity and nucleotide diversity of COⅠ gene in twelve geographic populations of C. ciliata

2.3 系统发育

种群系统发育树结果表明，法国、斯洛文尼亚与德国种群聚集在一起，表明这3个区域网蝽种群的亲缘关系非常近(图2)。在国内种群中，CS、SH、WH种群与这3个欧洲种群亲缘关系也非常近(图2)，如果我国网蝽来源于欧洲种群，CS、WH和SH都有可能是我国网蝽的首次入侵地。

图2 根据COⅠ基因核苷酸序列建立的网蝽种群系统发育树Fig. 2 Phylogenetic tree of the nucleotide sequence of COⅠ gene from C. ciliata populations

其他国内外网蝽种群并未呈现出地理区域的相关性，甚至在同一区域内聚集有多个遗传关系比较远的多个种群，例如斯洛文尼亚1种群与韩国种群聚在一支，而4种群与法国种群聚在一支，5种群却与我国保定1种群聚在一支。作为网蝽发源地的美国种群与国内保定1种群和扬州种群聚集在一支，说明它们的亲缘关系比较近(图2)。

2.4 分子变异

组群间AMOVA分析发现，其在组群间、组群内种群间和种群内个体间的固定系数，即FCT、FSC和FST的P值都小于0.05，说明上述分组群进行分析是合理的。组群间、组群内种群间和种群内个体间变异分别占总变异的26.03%、21.78%和52.19%(表3)。这样的结果说明，本研究中网蝽种群间变异主要是由种群内个体之间的变异引起的。

表3 12个网蝽地理种群的分子方差分析①Tab.3 AMOVA analysis of twelve geographic populations of C. ciliata

2.5 网蝽种群间的遗传关系

种群间遗传分化程度及基因流动分别采用FST及Nm来衡量，FST和Nm值代表的意义参考卢绍辉(2020)。4个组群间，华中与华东组群间网蝽具有较低的遗传分化，而华中与西南、华北组群间网蝽具有非常高的遗传分化，而华东与西南、华北组群间网蝽也具有较高程度的遗传分化(表4)。在66个种群对间FST值在-0.05～0.89中，其中有43个种群对间有显著差异(表5)。有10个种群对间的FST值小于0.05(表5)，说明这些种群对间遗传分化程度很低，甚至没有发生遗传分化，而其他33个种群对间的FST值大于0.05(表5)，说明这些种群对间遗传分化程度比较大。不同地理种群对间Nm值范围在-85.25～47.58间，有10个种群对间的Nm值大于4(表5)，暗示这些种群对间的基因交流经常发生； 17个种群对间Nm值在1～4间(表5)，表明这些种群对间有中等程度的基因交流； 其他种群对间的Nm值均小于1(表5)，表明种群对间几乎不发生基因交流。

表4 网蝽组群对的FST值(对角线下)和Nm值(对角线上)①Tab.4 Pairwise FST(below the diagonal) and Nm(above the diagonal) among group populations of C. ciliata

表5 12个网蝽地理种群对间的FST值(对角线下)和Nm值(对角线上)Tab.5 Pairwise FST(below the diagonal) and Nm(above the diagonal) among twelve geographical populations of C. ciliata

2.6 网蝽遗传距离与地理距离的相关性

对不同网蝽种群的遗传距离与地理距离间进行相关性分析，发现二者之间相关性系数为0.15，但不存在显著差异(图3)。说明网蝽不同地理种群间的遗传分化不是由地理隔离引起的，没有距离隔离的现象。网蝽种群位置的经度和纬度与单倍型多样性的相关系数分别为-0.127 和0.103，都没有表现出显著的相关性。

图3 网蝽种群遗传距离和地理距离间的相关性Fig. 3 Relevance between genetic and geographic distance from twelve populations of C. ciliata

2.7 不同区域种群历史动态和有效迁移率

对4个区域种群的历史动态进行了分析。Tajima’sD和Fu’sFS分别为负值且P<0.05时，显著偏离了中性突变，认为在历史上发生过种群扩张事件，否则认为种群在历史上没有扩张事件；PSSD>0.05，表明不能拒绝群体扩张的假说，即符合原来群体发生了扩张，而PSSD<0.05说明种群偏离突然扩张假设模型，即原来群体不能发生扩张(Tajima, 1989)。华东和华中组群内只有2种单倍型，因此这2个组群不能绘制错配图。在主要分布区，θ0和θ1值没有显著的改变，Tajima’sD及Fu’sFS值均不显著(表6)，说明中国区域网蝽种群没有发生扩张事件。同时，错配分布图有非常显著的双峰(图4C)以及PSSD值小于0.05(表6)，也证实中国区域网蝽种群没有发生扩张事件。

表6 网蝽种群历史分析①Tab.6 Demographic history parameters of C. ciliata

图4 中国主要分布区网蝽的错配分布Fig. 4 Mismatch distributions of the main distribution range of C. ciliata in ChinaA： 西南地区 Southwest China； B： 华北地区 North China； C： 中国主要分布区 Main distribution areas in China.

华北、西南和华东组群网蝽扩张前有效种群θ0和扩张后有效种群θ1间数值变化不大(表6)，表明这3个组群网蝽数量增加不明显，且Tajima’sD和Fu’sFS值均没有显著性(表6)，呈现双峰或者是下降型的错配分布图(图4A和B)，表明这3个地区网蝽未经历过明显的种群扩张(卢绍辉， 2020)。华中区域网蝽θ1值比较大，表明该地区网蝽种群增加较大(表 6)。同时，Tajima’sD为负值且Fu’sFS为显著的负值，表明该区域网蝽种群经历过扩张事件(表 6)。同时，PSSD值大于0.05(表 6)，也证实了该组群网蝽经历过种群扩张事件。

4个网蝽区域组群间的有效迁移率(M)都比较高，在421.9～648.0间(表7)。华中组群向华东、华北、西南的M值分别为580.9、610.4和598.4，而华东、华北、西南组群向华中组群的M值分别为449.8、421.9和455.7，表明华中组群迁出要高于迁入(表7)，这与网蝽种群历史分析中得到的“华中地区种群经历过种群扩张的事件”结论一致。

表7 中国网蝽组群间的迁移率①Tab.7 Immigration rate among group populations of C. ciliata in China

3 讨论

本研究对收集的12个地理种群的240个悬铃木方翅网蝽个体进行COⅠ基因测序，获得了240条761 bp的核苷酸序列，根据这些序列进一步分析了不同地理种群间的核苷酸序列差异及遗传分化。本研究测定的COⅠ基因核苷酸序列间未发现有碱基缺失和插入的情况。腺嘌呤(A)和胸腺嘧啶(T)的含量(67.6%)明显高于鸟嘌呤(G)和胞嘧啶(C)的含量(32.4%)，表现出A+T碱基偏倚性，符合昆虫线粒体构成的典型特征(Jermiinetal., 1994)。

多种因素影响生物向新栖息地扩散转移，通常这些因素也都可以导致生物种群发生遗传分化(陈圣宾等， 2010)。种群分化程度与其栖息地的距离存在着正相关，距离越远分化程度越大，这种现象称为距离隔离(Kimuraetal., 1964； Slatkin, 1993)。然而，本文发现中国分布区内的网蝽不同地理种群间并不存在距离隔离，说明中国区域网蝽发生的遗传分化不是由距离隔离产生的。对其他昆虫研究也发现了类似的结果，如草地螟(Loxostegesticticalis)(呼晓庆等， 2019)。研究发现，一些山脉和高原，比如云贵高原等，能明显影响物种的遗传结构(Yeetal., 2014; Xieetal., 2017)。GY、CQ和CD种群位于我国云贵高原区域，该区具有复杂的地形和独特的气候特点，被认为是生物高度富集的重要区域(Barthlottetal., 2007)。西南区域的GY、CQ和CD网蝽种群，具有很高的单倍体多样性，可能是因为更适应这里的环境所造成的，因为生物在适宜的栖息地往往有高的遗传多样性(Petitetal., 2003； 卢绍辉, 2020)。华中区域的WH、CS种群和华东区域的SH、HZ种群具有非常少的单倍型个数，并且2种群间未发现遗传分化。网蝽可以通过人类活动传播(鞠瑞亭等， 2011)，长江中下游地区是我国经济高度发达的区域之一，人员和货物流动频繁(姚瑞华等， 2014)，因此华中和华东组群间没有遗传分化很有可能是因为人类活动导致经常发生基因交流造成的。另外，华中和华东区域处于我国长江中下游，湖泊河流密集。每年的夏季为网蝽的发生高峰期(李峰奇等， 2018； 卢绍辉, 2020)，也是长江中下游地区洪水频繁发生的时期(许有鹏等， 2005; 刘蕾等， 2018)，网蝽有可能以悬铃木叶片为载体顺着河流进行扩散，加快了网蝽种群间基因交流。因此，地理障碍、人为传播和气候条件等都能影响种群之间的基因交流，从而导致种群间出现遗传分化。

网蝽地理种群之间的单倍型多样数量随着种群栖息地经度的减小而增加，且位于华东和华中组群具有非常低的单倍型多样性。研究发现，瓶颈和自然选择作用通常能影响入侵生物的遗传多样性，往往初始种群有更高的多态性，相反有比较低的多态性(Lanzavecchiaetal.， 2008; 卢绍辉, 2020)。笔者的研究结果表明网蝽在西南地区具有最高的遗传多态性，按此推理，西南种群应该为原始种群。然而，对4个组群的历史动态和错配图进行分析发现，西南、华北和华东组群并没有发生过种群向外扩张事件，仅仅华中组群发生了种群向外扩张事件。华中组群迁出率远高于迁入率，也暗示华中组群经历过种群迁出事件。结合网蝽2002年在华中地区首次报道的事实(Streito, 2006)，笔者推测长沙和武汉可能是虫源地之一。同时系统发育分析表明，WH、SH和CS种群聚集在一支，具有非常高的亲缘关系，因此上海也可能是虫源地。Yang等(2017)也发现网蝽首先入侵中国的东海岸，之后向全国扩散，结合这项研究结果，本研究中WH、SH和CS种群仅有SH种群居于东海岸，那么更有可能上海为网蝽的虫源地，并向四周扩散。根据系统发育分析结果，SH种群与欧洲种群的亲缘关系最近，推测我国网蝽首先是从欧洲入侵的。上海是我国对外开放的最大窗口，与欧洲贸易频繁发生，因此很有可能通过人类活动进行长距离传播。之后以上海为虫源地，通过人类活动的长距离传播和飞行等短距离传播在全国范围内进行扩散(李峰奇等, 2018)。然而，网蝽在我国并没有最先发现于上海(肖娱玉等， 2006; Streito, 2006)。入侵昆虫在新发地初期常因种群数量有限而往往不容易被发现, 如白蜡窄吉丁(Agrilusplanipennis)(Siegertetal., 2014)、喜菘蝽(Bagradahilaris)(Palumboetal., 2016)和斑翅果蝇(Drosophilasuzukii)(Bieńkowskietal., 2020)。网蝽也很有可能在上海首先发生入侵后，由于种群数量有限而未被发现，导致了这种现象的出现。

SH、WH和CS种群对间没有发生种群的分化(FST=0)，说明它们之间的亲缘关系非常近。依据地理距离推测网蝽迁移路线是上海→武汉→长沙。研究认为网蝽能够借助自然风向较远距离进行传播(李峰奇等, 2018; 卢绍辉, 2020)。夏季在副热带高压的作用下，我国沿海地区风向以东南风为主(李慧等， 2013)。因此网蝽可能依靠东南风向西北方向迁移，扩散至武汉、长沙等地区。另外，网蝽广泛分布于城市悬铃木上，交通工具从悬铃木树下通过时也可能将网蝽运走，从而产生了远距离的扩散。CS与GY、CQ、CD 种群FST值分别为0.36、0.42和0.79，说明西南地区网蝽种群可能是CS种群以贵阳为中转站迁移进入的，形成了长沙→贵阳→重庆→成都的迁移路线。然而云贵高原地理障碍能力非常强，外部种群很难进入(Wangetal., 2020)。再加上长沙与贵阳距离遥远，网蝽不太可能通过飞行进入贵阳，更可能是通过人为传播进入贵阳地区。研究发现，人类活动，比如船、高速公路、苗木运输等，都能远距离传播网蝽(Halbertetal., 1983; 虞国跃等, 2014; 鞠瑞亭等, 2011)。这些因素都可能导致网蝽扩散到我国西南地区。与其他种群相比，WH与ZZ种群间分化非常低(FST=0.10)，具有非常高的亲缘性，说明ZZ种群很有可能从武汉传入的。ZZ与SJZ、JN、TY种群间的FST值分别为0.69、0.57、0.70，说明ZZ种群与这3个地区的种群产生了显著的遗传分化。假设SJZ、JN和TY种群是从郑州地区迁入形成的，之后扩散到北京，那么BJ与ZZ种群应该也有显著的遗传分化，然而研究结果表明，ZZ与BJ种群的FST值为0.15，说明2个种群间遗传分化程度非常低。近年来国内外频繁的贸易往来，很可能在这些地区发生了网蝽多次入侵事件，从而导致这种现象的出现。多次入侵事件的发生，导致很难推测网蝽从郑州向北扩散的路线。因此，网蝽向西的主要迁移路线为上海→武汉→长沙→贵阳→重庆→成都； 向北的主要迁移路线为上海→武汉→郑州，而郑州向北迁移途径需要进一步研究。

4 结论

本文对12个悬铃木方翅网蝽种群遗传结构进行了研究，在不同组群间和不同地理种群间具有显著分化的遗传结构，而网蝽种群个体间的变异是种群发生遗传分化的主要因素，同时地理障碍、人类活动和气候条件均有可能影响网蝽种群间的遗传分化程度。我国主要分布区内网蝽可能以上海为中心向西、向北往全国进行扩散。