基于微卫星分子标记的浙江省不同滞育型迷卡斗蟋的遗传分化分析

王佰秋, 沈初泽, 武 雪, 谢慧聪, 宗靖淞, 余哲媛,白 芸, 李 恺,*, 何祝清,*

(1. 华东师范大学生命科学学院, 上海 200241; 2. 北京师范大学生命科学学院, 北京 100875)

迷卡斗蟋Velarifictorusmicado隶属于直翅目(Orthoptera)蟋蟀科(Gryllidae)斗蟋属Velarifictorus，从南至北广泛分布于东亚地区(Wangetal., 2020)，也有记录显示该物种已经成功入侵北美地区，并不断扩散(Bowles, 2018)。为适应不同纬度的气候条件，迷卡斗蟋演化出了两种截然不同的滞育型:卵滞育型和若虫滞育型。卵滞育型多分布于东亚的北部地区，以卵的形式过冬，春末孵化，秋季羽化成虫，交配产卵后死亡，完成整个世代；若虫滞育型多分布于东亚的南部地区，以若虫形式过冬，夏季羽化成虫，交配产卵后死亡，卵一般在夏末孵化，在冬季以大龄若虫形态过冬，翌年夏季羽化成虫，完成整个世代(Masaki and Walker, 1987; He and Takeda, 2013)。在中国的浙江省内发现有两种滞育型种群同域共存的现象(Wangetal., 2020)。线粒体基因分析结果显示，卵滞育种群与若虫滞育种群的线粒体基因已经发生了分化，相似滞育型的种群具有相似的COI基因，而不同滞育型的种群COI基因之间存在一定差异(Wangetal., 2020)。但线粒体基因是母系遗传，即在不存在突变的情况下后代的线粒体基因完全与母虫相同。虽然两种滞育型的迷卡斗蟋羽化成虫的时间不一致，但可能存在若虫滞育种群的成虫存活时间较长，卵滞育种群较早羽化成虫，两个种群的个体直接交配从而产生后代的情况。两个滞育型种群是否有交配，又关系到如何定义迷卡斗蟋这个物种。因为根据物种的定义，若两个种群长期不存在自然交配，且某一特征(包括生活史)存在显著差异，应将其划分为两个物种。因此，通过分子生物学技术对其进行更为深入的研究十分关键。

在蟋蟀总科中，已有不少基于微卫星技术对种群结构进行研究。例如Vanhala和Cottrell(2008)筛选出针蟋Nemobiussylvestris的25个多态性微卫星位点，用以分析森林破碎化对物种基因交流和遗传多样性的影响；Birge等(2007)筛选出针蟋属Allonemobius的10个多态性微卫星位点，用于研究其属内杂交情况；Bretman等(2008)筛选出可共同用于蟋蟀属Gryllusbimaculatus和G.campestris的微卫星位点；Tinghitella等(2011)对滨海油葫芦Teleogryllusoceanicus的7个微卫星位点进行分析，发现自其引入至夏威夷后，遗传多样性降低; Gupta等(2020)基于家蟋蟀的全基因组开发了9个微卫星位点，可用于种群结构分析。

综上所述，通过微卫星能够较好地研究属内以及种下的种群遗传结构。因此，本研究新开发获得11个微卫星位点，并据此对浙江省3个地区(天童山、大盘山和百山祖)迷卡斗蟋的卵滞育型和若虫滞育型种群进行测序分析，研究其遗传结构和基因交流情况，探究浙江省同域分布不同滞育型迷卡斗蟋种群之间的遗传分化的情况。

1 材料与方法

1.1 迷卡斗蟋样本的获取与基因组DNA的提取

本研究所用的迷卡斗蟋共6个种群134个标本分别采集于浙江省天童山、大盘山和百山祖国家级自然保护区，采集时间地点及种群名称具体见表1，根据采集时间及若虫成虫形态判别其滞育型，即：9和10月份采集获得成虫及6月份采集获得若虫为卵滞育型；反之，9和10月份采集获得若虫及6月份采集获得成虫为若虫滞育型。所有材料采集后首先保存在无水乙醇中，后于实验室更换为新的无水乙醇并保存在-20℃冰箱中。利用AxyPrep Genomic DNA Miniprep Kit(AXYGEN)试剂盒对蟋蟀后腿肌肉组织提取基因组DNA，利用Thermo Scientific NanoDrop One超微量分光光度计测定总DNA浓度。

1.2 线粒体基因扩增和单倍型网络关系图构建

迷卡斗蟋线粒体COI基因引物序列参考Simon等(1994)，LCO: 5′-GGTCAACAAATCATAAAGATAT TGG-3′；HCO: 5′-TAAACTTTCAGGGTGACCAAAAA ATCA-3′，由上海铂尚生物公司合成，由上海生工公司对扩增产物进行测序，扩增长度658 bp。PCR反应体系(50 μL): 2×Taq PCR Mix 25 μL, ddH2O 19 μL, 正反向引物(10 μmol/L)各2 μL, 基因组DNA 2 μL。扩增程序: 94℃ 4 min; 94℃ 30 s, 45℃ 30 s, 72℃ 30 s, 共35个循环；最后72℃ 5 min。利用ATGC V.7.0.2(Genetyx Corporation, Tokyo, 日本)检测测序结果的峰图并拼接序列，对单峰结果的序列进行拼接编辑，对于具有双峰的序列的个体重新扩增测序。利用MEGA V.7.0.14(Kumaretal., 2016)中的Clustal W进行序列比对。利用DAMBE V.5.3.8(Xia and Xie, 2001)进行碱基替换饱和性检测，结果显示ISS

1.3 线粒体基因单倍型网络关系图构建

使用DnaSP V.5.10.01(Librado and Rozas, 2009)软件将COI基因的序列fasta文件按地理种群分组，生成单倍型文件.nex和.arp文件。在ARLEQUIN中导入.arp文件，计算各种群的单倍型频率。将单倍型文件统计为.txt文件。在PopART(Leigh and Bryant, 2015)中导入单倍型文件.nex和种群的单倍型频率.txt文件。利用中介(median-joining, MJ)网络方法构建单倍型网络关系图以检测不同组之间的系统发育关系。

1.4 微卫星分子标记开发及筛选

提取得到迷卡斗蟋的基因组DNA后，由上海派森诺公司进行测序，采用Illumina MiSeq系统收集数据并产生标准的数据文件格式，数据经过过滤后，共获得1.31 Gb数据量。根据引物重复单位和多态性的类型选择101对引物，其中包含16对二核甘酸重复模体的引物，36对三核甘酸重复模体的引物，44对四核甘酸重复模体的引物，2对五核甘酸重复模体的引物，3对六核甘酸重复模体的引物。选择4头迷卡斗蟋(包括2头卵滞育个体，2头若虫滞育个体)的基因组DNA进行微卫星位点的扩增，初步筛选出能够扩增出清晰条带并具有多态性的位点。

通过种群个体再次筛选初筛的引物以检测位点的遗传多样性参数。采用荧光PCR扩增微卫星基因，通过毛细管电泳检测扩增产物的条带大小。选择浙江天童山若虫滞育种群，对20头个体进行微卫星位点的再次筛选。由上海生工公司合成3种荧光标记的引物，分别为FAM, HEX和NED。引物信息见表2。采用96孔板扩增，不同荧光标记且大小相差50 bp以上的扩增产物混合送样以节约成本，PCR产物在ABI 3730XL测序仪(Applied Biosystems)上进行基因分型，在GeneMapper v4.0(Applied Biosystems)上进行分析。

1.5 遗传多样性分析

在Genalex v6.5 (Peakall and Smouse, 2012)中计算种群各位点的等位基因数(number of alleles,Na)、有效等位基因数(number of effective alleles,Ne)、观测杂合度(observed heterozygosity,Ho)、期望杂合度(expected heterozygosity,He)和近交系数(inbreeding coefficient,Fis)；在Arlequin(Excoffier and Lischer, 2010)中进行哈德-温伯格平衡(Hardy-Weinberg equilibrium, HWE)和种群的连锁不平衡(linkage disequilibrium, LD)检测。 在Microchecker里检测是否存在无效等位基因，在Arlequin v3.5.2.2(Excoffier and Lischer, 2010)中计算遗传分化指数(Fst)、种群内近交系数(Fis)和基因流(gene flow,Nm)。在Cervus v3.0(Kalinowskietal., 2007)中计算多态信息含量(polymorphism information content,PIC)。

1.6 遗传分化分析

将筛选好的微卫星位点应用于6个种群共134个个体中。为检测种群间的分化程度并分析分歧来源，在Arlequin中计算遗传分化指数(Fst)并进行AMOVA分析。Fst值是利用成对遗传距离计算的，统计意义为10 000个置换。AMOVA分析时采用两种分析方法对种群进行分组并分析遗传分化来源，(1)按照地理位置划分：百山祖种群BSZ-E和BSZ-N为一组(group 1)，大盘山种群DP-E和DP-N为一组(group 2)，天童山种群TT-E和TT-N为一组(group 3)；(2)按照滞育型划分：BSZ-E, DP-E和TT-E为卵滞育组(group 4)，BSZ-N, DP-N和TT-N为若虫滞育组(group 5)。

1.7 遗传结构分析

利用Structure v2.3.4(Pritchardetal., 2000)软件进行种群结构分析，在线 http:∥taylor0.biology.ucla.edu/struct_harvest/确定最佳的种群组数，即最佳K值；利用Clumpp v1.1(Jakobsson and Rosenberg, 2007)进行Clusters重复抽样分析以减少结果的偶然性；利用Distruct v1.1图形化显示种群遗传结构。

2 结果

2.1 基于线粒体基因COI的单倍型网络关系

结果显示迷卡斗蟋134个个体共分为18个单倍型，6个种群显著分为两个世系，且同一滞育型的种群聚为一个支系，即卵滞育型种群BSZ-E, DP-E和TT-E聚为一支，若虫滞育型种群BSZ-N, DP-N和TT-N聚为另一支(图1)。

图1 浙江省3个地区不同滞育型迷卡斗蟋6个种群基于mtCOI基因的单倍型中介网络关系Fig. 1 Median-joining haplotype network based on mtCOIgene of six populations of Velarifictorus micado withdifferent diapause types from three areas in Zhejiang Province种群信息见表1。For population information, see Table 1. 下同The same below.

2.2 微卫星位点筛选

在101对微卫星引物中，共筛选出16对能够扩增出具有多态性的清晰条带，进一步检测这些微卫星位点的遗传多样性。结果显示16个微卫星位点中有4个位点(Vm11, Vm27, Vm92和Vm95)偏离哈德-温伯格平衡，Microchecker检测显示位点Vm11和Vm92存在无效等位基因，这很可能是造成其偏离的原因。位点Vm59的遗传多态性低，其观测杂合度为0.250，期望杂合度为0.219。其余11个微卫星位点遗传多态性较高，可用于后续分析(表2)。

2.3 遗传多样性

在6个种群共134个个体中，利用11对微卫星引物检测到每个位点的等位基因数(Na)为10.333(位点Vm37和Vm41)～17.833(位点Vm101)，平均13.803；观测杂合度(Ho)为0.327(位点Vm37)～0.898(位点Vm101)，平均0.743。 位点Vm54(Fis=-0.012)起表现为存在杂合子冗余，Fis值0.006(位点Vm101)～0.631(位点Vm37)表示杂合子缺失(表3)。 等位基因数(Na)的平均值范围为(6.000±2.966)(BSZ-E)～(18.000±2.449)(TT-E)个。 所有种群等位基因数(Na)平均为13.803。BSZ-E的观测杂合度(Ho)最低，为0.514±0.292，TT-N的Ho值最高，为0.850±0.100，平均0.743。Fis值为0.06990(TT-N)～0.18849(TT-E)，表示杂合子缺失(表1)。

表2 浙江省天童山迷卡斗蟋若虫滞育种群的16个微卫星位点信息及遗传多样性Table 2 Information of sixteen microsatellite loci in nymphal diapause populations of Velarifictorus micadofrom Tiantongshan, Zhejiang Province and their genetic diversities

续表2 Table 2 continued

表3 浙江省3个地区不同滞育型迷卡斗蟋6个种群基于11个微卫星位点的遗传多样性信息Table 3 Genetic diversity information of six populationsof Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on11 microsatellite loci

2.4 遗传分化和基因流

迷卡斗蟋种群间的Fst值范围为0.01340(DP-N与TT-N)～0.21087(BSZ-E与BSZ-N)，其中BSZ-E与其他种群Fst差异大(P<0.001)。BSZ-E与BSZ-N分化程度最高，Fst值为0.21087，与其他种群的之间的Fst值为0.15783～0.19513。BSZ-N与其他种群的Fst值差异较大，Fst值为0.03242～0.05796(P<0.001)(表4)，所有种群之间的Nm值均大于1，为1.87111～36.81399(表4)，表明基因交流频繁。

表4 浙江省3个地区不同滞育型迷卡斗蟋6个种群间基于11个微卫星位点的基因流Nm值(上三角)和遗传分化指数Fst值(下三角)Table 4 Gene flow Nm values (above the diagonal) and pairwise genetic differentiation coefficient Fst values(below the diagonal) between six populations of Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on 11 microsatellite loci

基于地理位置和滞育型两种分组方法的AMOVA分析结果显示，总体遗传变异均来源于个体之间。基于地理位置分为3组(group 1-3)，组间差异占比-1.87%，组间种群内差异占比9.35%，组内种群间差异占比12.69%，个体间差异占比79.83%。基于滞育型类型分为卵滞育组和若虫滞育组两组(group 4和5)，组间差异占比1.27%，组间种群内差异占比7.02%，组内种群间差异占比12.60%，个体间差异占比79.12%(表5)。因此种群间差异主要来源于个体之间，与地理隔离或者滞育型关系无显著相关性。

表5 浙江省3个地区不同滞育型迷卡斗蟋6个种群基于11个微卫星位点的AMOVA分析Table 5 AMOVA of six populations of Velarifictorus micado with different diapause typesfrom three areas in Zhejiang Province based on 11 microsatellite loci

2.5 遗传结构

使用Structure对浙江6个种群进行遗传结构分析，基于基因型频率分析得最佳分组K值为3，当K=3时，DeltaK值最大，即6个种群分为3组，主要划分为BSZ-E一组，BSZ-N中有一半个体与DP-E, DP-N, TT-E和TT-N为一组，其余个体为一组(图2: A)。当K=2时，主要划分为BSZ-E一组，其余个体为一组(图2: B)。

图2 基于Structure v2.3.4分析当K=3(A)和K=2(B)时浙江省3个地区不同滞育型迷卡斗蟋6个种群的遗传结构Fig. 2 Genetic structure of six populations of Velarifictorus micado from three areas in Zhejiang Provinceanalyzed by Structure v2.3.4 with K=3 (A) and K=2 (B)

3 讨论与结论

3.1 迷卡斗蟋遗传多样性

等位基因数(Na)、观测杂合度(He)、期望杂合度(Ho)和多态信息含量(PIC)是反映种群遗传多样性的重要指标(Maudetretal., 2002)。本研究共成功筛选出迷卡斗蟋11个具有多态性的微卫星位点，每个位点的平均等位基因数(Na)为13.803，平均观测杂合度(He)为0.743，PIC值范围为0.826～0.928，筛选的微卫星具有较高的多态性(表3)。利用这11对微卫星位点分析了浙江省天童山、大盘山和百山祖3个地区两种滞育型共6个种群的遗传多样性。 每个种群的等位基因数为(6.000±2.966)～(18.000±2.449)个，其中BSZ-E的最低(表1)。所有种群的观测杂合度(Ho)平均为0.743(表3)，其中BSZ-E的观测杂合度最低，为0.514±0.292，由此可见BSZ-E的遗传多样性不高，可能是由于采集时仅局限于一个区域；TT-N的观测杂合度(He)最高，为0.850±0.100，与此对应的TT-N的自交系数(Fis)最低(表1)。

3.2 迷卡斗蟋遗传分化和基因交流

Wright(1978)指出Fst<0.05表示种群间分化程度较低。BSZ-E与其他种群的遗传分化系数(Fst)的范围为0.15783～0.21087，其次为BSZ-N与其他种群的，这表明百山祖的种群与其他地区的迷卡斗蟋种群遗传距离大(表4)。Wright(1978)认为当Nm>1时，种群间存在基因交流，减缓了种群间的分化。基因流结果表示，各种群之间存在着频繁的基因交流，虽然其中BSZ-E与其他种群的基因交流程度较低，为1.87111～2.66805，但是其他种群之间的基因交流十分频繁，其中DP-N与TT-N的基因交流最频繁，基因流值为36.81399(表4)。可见地理位置更近且属于同种滞育型之间的迷卡斗蟋基因交流最大。

3.3 迷卡斗蟋遗传结构

对迷卡斗蟋的遗传结构进行分析，当分为两组时(K=2)，BSZ-E独立为一组，其他种群聚为一组。当进一步分为3组时(K=3)，BSZ-N中的一半与其他种群聚为一组，剩下的一半单独为一组(图2)，可见，百山祖地区与其他地区的迷卡斗蟋之间也许由于地理距离纬度的原因基因交流减少。而AMOVA方差分析结果显示，基于地理位置分为3组时，个体间差异占比最高(表5)，说明组与组之间、种群与种群之间并没有显著的遗传分化，迷卡斗蟋间的遗传差异主要来源于种群内部个体之间的差异，而不是地理上的差异或者不同的生活周期。综上研究结果表明，不同滞育型之间的种群之间仍存在着频繁的基因交流。

3.4 迷卡斗蟋生活史分化

基于线粒体基因的种群地理学分析显示，迷卡斗蟋的mtCOI基因共18个单倍型，分化为两个单倍型组(图1)，分别对应于卵滞育种群和若虫滞育种群，这与Wang等(2020)的研究结果一致。然而利用微卫星数据的种群遗传学分析显示，具有两种不同滞育型的种群间并未呈现显著的遗传分化，同域分布的卵滞育和若虫滞育迷卡斗蟋存在广泛的基因交流(图2)，这说明在同一地区两种不同滞育型的种群的成虫阶段可能仍具有相当程度的重合。有一种可能性为若虫滞育种群在夏季成虫后，一直存活至秋季，与秋季刚成虫的卵滞育种群成虫交配。综合前人线粒体基因分析的结果，我们推测调控迷卡斗蟋生活史分化的主效基因可能位于线粒体基因组上，或很大程度伴随母系遗传。类似的现象在部分其他的昆虫类群中也有报道，如大猿叶虫甲Colaphellusbowringi不同地理种群杂交后代的滞育率主要取决于其母本(Laietal., 2008)，红头毛蛾Zygaenatrifolii父本对后代滞育率完全没有贡献(Wipking and Kurtz, 2000)。所以不同滞育型的两者“杂交”后，卵由母虫在其发生季节产出，只有与母虫保持一致的滞育型，才能继续存活，因此mtDNA清晰地分化为两支。假设“杂交”后代保持雄虫的滞育型且能够继续存活并产生后代，在mtDNA结果中应该能发现交错的现象，而非如今清晰地划分为两组。综上所述，两个种群存在频繁的基因交流，但“杂交”后代只保持了母虫的滞育型且继续繁衍。这既解释了为何同域不同滞育型种群之间存在频繁的基因交流，而mtCOI基因却能够清晰地分为两支。将来，要理清不同滞育型迷卡斗蟋的遗传分化的调控机制，还要对同域分布种群进行进一步的GWAS分析。