基于 28S rDNA D2区序列的陕南凹缘菱纹叶蝉分子系统学与遗传多样性研究

杨金宏，孔卫青

(安康学院，陕西省蚕桑重点实验室, 陕西安康 725000)

凹缘菱纹叶蝉(HishimonussellatusUhler)俗称绿头菱纹叶蝉，隶属于昆虫纲(Insecta)半翅目(Hemiptera)叶蝉科(Cicadellidae)殃叶蝉亚科(Euscelinae)[1]。凹缘菱纹叶蝉是桑树黄花型萎缩病、枣疯病等的重要媒介昆虫，在中国分布广泛[2]。成虫和若虫均可造成危害。它们以口针刺吸病树汁液，亦将病原吸入体内，后者移至唾液腺并繁殖；再次取食时病原菌传入新寄主，从而导致健康植株发病。桑树染病后通常会出现花变叶、矮化、黄化、衰退、簇生、丛枝、小叶等症状，严重影响蚕桑产业的健康发展。凹缘菱纹叶蝉每年4月下旬至10月下旬发生危害，1 a达3～4代，繁殖速度快，世代周期短，世代重叠，危害严重，且与养蚕周期重合，难以防治。

秦巴山区有众多的小盆地和山间谷地相连接，安康市汉滨、石泉、汉阴、平利、紫阳等县(区)和汉中市西乡、勉县等均位于该地理区域，形成陕南蚕桑产业主产区。各地自然地理条件差异较大，具有很高的环境异质性。桑园种植多分布于坡地和山地，可能阻碍凹缘菱纹叶蝉扩散，形成地理隔离，进而导致不同种群间受到不同的选择压力而出现种群分化。随着近年来国家退耕还林、土地流转政策以及蚕桑产业生产模式的变更，较多的连片桑园和强村大户出现，这为凹缘菱纹叶蝉扩散和传播提供了条件。为了解陕南桑园凹缘菱纹叶蝉种群的遗传信息，分析其分布特点，本研究利用 28S rDNA基因D2区对其种群遗传结构和遗传分化进行分析。

28S rDNA存在于编码细胞质核糖体大亚基中，在生物体内含量较高，在进化过程中比较保守，碱基替换率低，是研究生物系统发育较好的分子标记[3-4]。同时在 28S rDNA基因的保守序列中，含有12个高变区D1～D12，可解释从属到科级水平上的系统发生关系以及属种间的亲缘关系[5-8]，是非常有用的遗传进化分析工具。Dietrich等[9]用 28S rDNA D2～D10区序列分析角蝉总科(Membracoidea)与叶蝉科的亲缘关系。戴仁怀等[10]利用 28S rDNA基因 D2区、 16S rDNA序列和形态特征对中国角顶叶蝉亚科(Deltocephalinae)的系统发育进行研究。目前，未见利用 28S rDNA基因对凹缘菱纹叶蝉进行分子系统学和遗传多样性研究的报道。本试验尝试这项研究，其结果将促进对陕南蚕桑产区凹缘菱纹叶蝉的种群遗传结构及其分布概况的了解，并为凹缘菱纹叶蝉相近种的分类和快速分子诊断提供依据。

1 材料与方法

1.1 样本采集

2017年9月至10月，用黄板诱虫法采集秦巴山区陕南各蚕桑区凹缘菱纹叶蝉成虫。采集的新鲜标本浸泡于φ=75% 乙醇中，于12 h内置于-20 ℃的冰箱保存，备用。样品采集地点及相关信息见表1。

表1 陕南凹缘菱纹叶蝉种群样本采集地及其数量Table 1 Locality and sample size of Hishimones sellatus populations in South Shaanxi

1.2 基因组DNA提取

1.3 PCR扩增及测序

目标片段扩增用28S D2～D3引物[9]完成，引物序列为：28SD2F：5′-AGTCGKGTTGCTT-GAKAGTGCAG-3′；28SD2R：5′-TTCAATT-TCATTKCGCCTT-3′。反应总体系50 μL，包括基因组DNA 2.0 μL，10×PCR扩增缓冲液5.0 μL，Mg2+(25 mmol/μL)5.0 μL，dNTPs(10 mmol/μL)1.5 μL，上下游引物(10 pmol/μL)各1.0 μL，高保真TaqDNA聚合酶(5 U/μL)1.0 μL，加双纯水补足至50.0 μL。PCR扩增程序为[8]：94 ℃预变性3 min，36×(94 ℃变性1 min，52.2 ℃退火 1 min，72 ℃延伸50 s)，72 ℃终延伸10 min。取 3.0 μL PCR扩增产物进行10 g/L的琼脂糖凝胶电泳，检测成功后直接送生工生物工程(上海)股份有限公司进行纯化测序。

1.4 序列多样性分析

以Clustal X 1.83软件[11]基于邻接法(Neighbor-joining)对所得序列进行比对，去除两端不准确序列后，最终输出fasta格式的对齐序列。利用Dnasp 61101软件[12]统计各序列间的碱基组成、变异位点，以及基于对齐序列的种群核苷酸多样性Pi[13]、单倍型多样性Hd[Hd=n(1-∑Pi2)/(n-1)，其中n为样本数]和种群遗传分化系数等，利用Tajima’sD方法对序列进行中性检验，并根据数值服从β分布的特点，利用随机抽样模拟的种群中性进化为零假设，通过与观测值的比较进行显著性检验[14]。

1.5 种群间遗传差异分析

利用Arlequin 3.5.22软件[15]，基于Locus-by-Locus模型分析种群变异的分子方差AMOVA(Analysis of molecular variance)以及种群间的固定系数Fst(F-statistics)，以成对单倍型差异比较法(Pairwise difference)计算种群间遗传距离。运行Dnasp 61101软件[12]计算种群间的基因流系数Nst(Population subdivision at the nucleotide level)和种群间的地理距离(D)，Nst= 1 /(1 + 2Nm)；D= 111.199{(a1-a2)2+ (b1-b2)2cos2[(a1+a2) / 2]}1/2，其中Nm为有效迁移个体数，a和b分别代表经度和纬度，并进一步计算D的自然对数值以及与Nst之间的决定系数(R2)。

采用Bandelt等[16]的Median-joining模型绘制单倍型间的网络进化图，以叶蝉科角顶叶蝉亚科的Dwightlaacutipennis作外群，MEGA 7.0.26软件[17]运行ML极大似然法构建各个单倍型的系统进化树，采用的碱基替换模型为GTR+G，各分支置信度以1 000次自展重复进行检验。

2 结果与分析

2.1 凹缘菱纹叶蝉 28S rDNA D2区序列的特征及多样性

经检测，提取的132头凹缘菱纹叶蝉的总DNA全部合格； 28S rDNA基因D2区片段的PCR扩增和纯化测序，均获得良好的测序结果。测序序列比对并去除2端不准确序列后，得到包括缺失位点的有效片段695 bp(Alignment score=138 469)，序列比对结果显著性检验P值为 0.001～0.041。综合所有的多样性位点共定义8个单倍型(图1)，其中存在多态性位点63个，突变67个；其中4个位点存在3种变体，59个位点有2种变体。分析各单倍型序列数和地理分布，H_1有98条序列(表2)，占序列总数的 74.2%，且都存在于所有7个地理种群中，为优势共享单倍型；H_2和H_3仅在1个种群中出现，为私有单倍型。各种群内单倍型的分布从1种～6种不等，其中3号种群最多6种，且出现29个个体共享1个单倍型，具有较高的保守性。核苷酸位点多样性Pi和单倍型多样性指数Hd分别从核苷酸的突变位点数和差异基因型两个层面展示种群的遗传多样性。7个种群的Pi分布为0～ 0.019 03，Pi最高的是3号种群，Hd为0～ 0.581 82，7号种群最高(表2)。Tajima’sD中性检验是基于等位基因频率和单倍型数目变化差异的一种种内自然选择检验方法，对每一个种群和总群体进行检验(表2)，发现除4号种群大于0外，其他种群均为负数，其中5号种群为 -2.222 40，且显著性检验P<0.01，说明该种群受到强烈的负选择影响[14]。

图1 28S rDNA基因D2区序列单倍型及变异位点信息Fig.1 Haplotypes and variation sites of 28S rDNA D2 region sequence

地理种群Geographic population核苷酸多样性(Pi)Nucleotide diversity单倍型多样性(Hd)Haplotype diversity单倍型数量及其分布Haplotype number and distribution中性检验Tajima’s DP值P value10.015 240.466 67H_1(15), H_7(2), H_8(4), -0.520 080.320200H_1(7)30.019 030.553 91H_1(29), H_4(3), H_5(2), H_6(3), H_7(3), H_8(4)-0.213 670.42940.013 680.485 29H_1(12), H_3(3), H_8(2)0.260 430.67350.006 510.166 67H_1(11), H_5(1)-2.222 40*0.00160.008 420.278 95H_1(17), H_4(1), H_6(2)-0.798 470.27970.013 260.581 82H_1(7), H_2(2), H_5(2)-0.191 360.445合计 Total0.013 780.440 40H_1 ～ H_8-0.742 410.450

注：“( )”中的数字代表叶蝉数目。“*”表示差异显著。

Note: Numbers in “( )” represent the number of leafhoppers.“*” indicate significant difference.

2.2 种群遗传差异分析

Excel软件计算基因流系数Nst和种群地理距离D的自然对数间的决定系数(R2)，可知R2=0.019 7(P<0.001)，接近于0，说明该基因在本研究范围内与各个种群之间的基因流动和地理距离不相关，可能存在网状基因交流。

2.3 单倍型进化分析

Median-joining模型是在maximum-parsimony模型[20]和Kruskal’s算法[21]基础上发展的，据此绘制的进化图呈网状分布。本研究的凹缘菱纹叶蝉 28S rDNA D2区8种单倍型间的进化网络图呈线状分布(图2)，在所有的8个单倍型中，只有H_4和H_1是种内网络进化图中的中间单倍型[22]，其他单倍型均处于进化网络边缘，通过6个缺失的中间单倍型(mv1～mv6)连接，这与各单倍型在各个地理种群的分布结果一致(表2)，说明凹缘菱纹叶蝉的地理聚集性不强，交叉散布于整个进化网络图中。

灰色圆圈代表单倍型 Grey circle represents the haplotype; 灰色圆圈面积代表单倍型相对频率 The size of grey circle represents the relative frequency; 实心黑点代表中间单倍型 Solid black spots represent intermediate haplotype; 文本框内数字代表单倍型间的突变位点数 Boxed numbers represent the mutation sites between haplotypes

图2 基于 28S rDNA基因的D2区的8个单倍型的Median-joining 网络关系图
Fig.2 Median-joining network of eight haplotypes based on 28S rDNA Gene D2 region sequences

进一步以叶蝉科角顶叶蝉亚科的Dwightlaacutipennis的 28S rDNA D2区(JX845493.1)作外群，构建单倍型的ML系统进化树，结果其 log likelihood值为-1 474.09。自展支持值bootstrap法是对初始DNA数据进行放回抽样生成伪样本，进而估计节点可靠度的一种非参数方法，本研究节点的自展支持值均大于90，且与网络进化图的聚集特征一致，进化关系较为可靠。根据外群的节点位置，进化树中单倍型H_1、H_2、H_7、H_8形成一个分支，H_3～H_6形成一个分支，其中H_7的枝长最长，碱基替换率最高，说明其与其他单倍型的实际进化距离最远(图3)。

枝长代表碱基替换率 Branch lengths measured in the number of substitutions per site

图3 基于 28S rDNA D2区的单倍型的ML系统发育树
Fig.3 ML phylogenetic tree of haplotypes based on 28S rDNA D2 region

3 讨 论

种群为适应不同地域的自然环境而独自演化，可能出现地理种群的遗传分化和形态变异[23]，基于 28S rDNA基因D2区对陕南桑园凹缘菱纹叶蝉各种群和总群体进行Tajima’sD中性检验和显著性分析，结果表明该种群和总群体近期没有受到正选择的影响，5号种群受到负选择的影响，种群中存在较多的低频率等位基因，而其他种群等位基因频率处于稳定状态，遵循中性进化模型[14]。陕南桑园凹缘菱纹叶蝉 28S r DNA基因D2区的核苷酸多样性指数处于较低水平，变异主要来自种群内且差异不显著，表明该群体的大部分差异是随机分布，且存在广泛的基因交流。

选择恰当的外群即亲缘关系合适的外源物种，是进化树分析的关键，外群选择不当，容易造成进化树的拓扑结构发生变化[24]。叶蝉科被分为16个亚科，亚科内部的亲缘关系知之甚少，且不同亚科间可能存在许多并系[25]。因此本研究首先以蝽科昆虫为外群，分析整个叶蝉科的DNA序列，分析获得角顶叶蝉亚科是分析凹缘菱纹叶蝉的合适外群。以此构建的进化树各节点的自展支持值均大于90%，且与网络图一致，说明该外群的选择是合适的。

基因流是种群间基因通过个体或其配子迁移到其他种群而导致等位基因频率发生变化的现象，是影响种群间和种群内遗传变异的主要因素[26]。本研究种群内的变异远大于种群间，说明种群间的基因流增加了种群间的相似度和种群内的变异度，也是陕南桑园凹缘菱纹叶蝉没有形成地理格局的原因。基因流和地理距离之间的相关性分析显示，该基因区域各个种群之间的基因流动和地理距离不相关，出现这种状况的原因，一方面可能与凹缘菱纹叶蝉具有一定的飞行能力和迁飞的习性有关[27]，成虫和若虫能够适应各种桑园生态环境，能够渐次或连续发生一些近距离的扩散，导致不同地理种群间发生遗传物质的交流。另一方面可能是人类生产活动比如桑苗等宿主植物远距离运输等的作用。