基于线粒体16S rRNA基因的舟山群岛碧蛾蜡蝉种群遗传多样性分析

李 奎, 周传江, 梁爱萍

(1. 河南师范大学水产学院, 河南新乡 453007; 2. 中国科学院动物研究所, 动物进化与系统学院重点实验室, 北京 100101; 3. 中国科学院大学生命科学学院, 北京 100049)

岛屿(island)是指散布在江河湖海中的、四面环水的、且在高潮时仍暴露于水面上的陆地。岛屿地理结构分明，且被海水阻隔，形成了彼此独立的微型生态位，是探究推动物种形成、生物多样性、适应辐射及跨海扩散等事件的天然实验室(Demosetal., 2016)。陆桥岛(continental islands)在历史上曾与大陆相连，海平面的抬升或地质构造事件使其与邻近的大陆隔离(Paulay, 1994)。岛上的生物曾是邻近大陆上生物区系的一部分，长期的隔离阻断了岛屿与邻近大陆种群间的基因流，使其形成独立的进化单元。这很可能造成岛屿上不同种群及岛屿与大陆种群间遗传多样性的差异。Frankham(1997)及Yamada和Maki(2012)采用等位酶和微卫星标记，对分布于日本本州岛东部沿岸及伊豆群岛上的水仙花Weigelaoraeensis遗传多样性分析表明，大陆种群的遗传多样性明显高于岛屿种群，且随着岛屿与大陆距离的增加，其遗传多样性逐渐减小。

舟山群岛位于我国浙江省东北部，北连上海，南接宁波，东邻公海，总面积约2.22万平方公里，是我国最大的群岛，由1 390个大小不同的岛屿组成(宋亚民，2001)。从地质构造上讲，舟山群岛隶属于闽浙隆褶带的东北端，是浙江天台山脉沿东北方向东海的延伸(叶青超和黄金森, 1958)。受第四纪冰期的影响，更新世末期东海海平面抬升，造成舟山群岛开始与大陆隔离；直到全新世时期(大约7-9 kya)，舟山群岛与大陆完全分隔(Wangetal., 2014)。因此，舟山群岛是探究岛屿隔离对种群多样性影响的良好模型。然而，有关舟山群岛上种群多样性的研究多局限于邻近大陆的岛屿上的个别物种(Huetal., 2006; Chenetal., 2008)。近期，对分布于整个舟山群岛上的黑斑侧褶蛙Pelophylaxnigromaculatus的遗传多样性评估表明，岛屿隔离致使岛屿种群遗传多样性相对较低(Wangetal., 2014)。目前，对舟山群岛上分布的昆虫的遗传多样性和种群历史动态的研究工作还不多(Lietal., 2020)，岛屿隔离效应对分布于该群岛上的昆虫的遗传多样性的影响还有待深入的探讨。

碧蛾蜡蝉Geishadistinctissima隶属于蜡蝉总科(Fulgoroidea)蛾蜡蝉科(Flatidae)碧蛾蜡蝉属Geisha，广泛分布于我国南方诸省(Medler, 1992; 王应伦等, 2005)。在样本采集过程中，本研究发现碧蛾蜡蝉广泛分布于舟山群岛上。尽管Li等(2020)对分布于该区域的碧蛾蜡蝉的种群遗传多样性和进化历史进行了探讨，但其所选的分子标记(COX1, COX2-COX3, CYTB)均为半翅目昆虫线粒体基因组中最为保守的区段(Wangetal., 2015)，无法充分解析岛屿隔离对碧蛾蜡蝉种群遗传多样性及进化历史的影响。线粒体16S rRNA不仅含有高度保守区，还包含中度保守区和高度可变区，被广泛应用于分子系统学及群体遗传多样性研究中(黄卫东等, 2019)。因此，本研究采用16S rRNA作为分子标记，以评估岛屿隔离效应对分布于舟山群岛及邻近大陆的碧蛾蜡蝉种群遗传多样性的影响。

1 材料与方法

1.1 供试材料

作者于2017年7-9月在舟山群岛境内由西南向东北方向共选择了11个较大的有人类居住的岛屿作为碧蛾蜡蝉采样点。其中舟山本岛和六横岛因面积较大，分别设置了2个采样点(舟山本岛：桥头施和展茅镇；六横岛：峧头和台门)。此外，本研究还分别在宁波市瑞岩寺国家森林公园(RYS)及象山县的贤痒镇(XY)设立了采样点。研究标本采用手抄网进行随机采样，详细的采样信息见图1和表1，与已发表的研究(Lietal., 2020)相比，本研究增加了5头样本。采集到的标本整理编号后置于无水乙醇中，并存放于-80℃冰箱中备用。

图1 舟山群岛地理位置及碧蛾蜡蝉样本采集点分布Fig. 1 Geographical location of the Zhoushan Archipelago and the sampling sites of Geisha distinctissima地理底图于2017年下载自国家基础地理信息系统http:∥www.ngcc.cn/ngcc/;采集样点种群代码和全称同表1，括号内的数字表示样本数。The geographical base map was downloaded in 2017 from the National Fundamental Geographic Information System (http:∥www.ngcc.cn/ngcc). The population codes and full names of sampling sites are the same as in Table 1, and the value in brackets is the number of individuals collected in the site.

1.2 舟山群岛碧蛾蜡蝉种群16S rRNA基因扩增

剪取碧蛾蜡蝉3对足，采用动物组织与血液DNA提取试剂盒(天根，北京)，根据标准步骤提取该物种的总基因组DNA。利用Primer Premier 6.0参照NCBI上已有的碧蛾蜡蝉线粒体基因组(GenBank登录号: NC_012617.1)序列设计16S rRNA基因引物(正向引物: 5′-TTAATCCAACAT CGAGGTC-3′；反向引物: 5′-TACCTTTTGTGTCA GGGTT-3′)，并由上海生工生物工程有限公司合成。反应体系(25 mL): 样本DNA 0.8 μL, 正反向引物(10 mmol/L)各0.8 mL, LA Taq酶混合液12.5 μL, 双蒸水10.1 μL。扩增反应条件: 95℃预变性5 min; 95℃变性30 s, 55℃退火30 s, 72℃延伸1 min，共35个循环；最后，72℃终延伸5 min，4℃保存。扩增产物经分离纯化后送至擎科生物公司(北京)进行双向测序。

1.3 生物信息学分析

将获取的碧蛾蜡蝉种群16S rRNA基因序列采用Lasergene.v7.1软件包中的SeqMan程序进行比对拼接，同时核查所获16S rRNA基因序列是否有错配、双峰及杂峰等现象，以确保序列的准确性。采用MEGA7(Kumaretal., 2016)中的MUSCLE进行多序列比对。比对完成的所有个体的rRNA基因序列按照BankIt(www.ncbi.nlm.nih.gov/WebSub/?tool=genbank)使用说明在线上传至NCBI中。

表1 舟山群岛碧蛾蜡蝉种群遗传多样性信息Table 1 Genetic diversity of Geisha distinctissima populations in the Zhoushan Archipelago

1.4 碧蛾蜡蝉种群遗传进化分析

使用DnaSP v5(Librado and Rozas, 2009)对不同岛屿上碧蛾蜡蝉种群的16S rRNA基因单倍型个数(haplotype number,Nh)、单倍型多态性(haplotype diversity,Hd)以及核苷酸多态性(nucleotide diversity,π)进行评估，以衡量种群间的遗传多样性差异；结合先前的研究结果(Wangetal., 2014)，采用Spearman检验分别对种群单倍型多态性及核苷酸多态性与岛屿-大陆间距离和岛屿被隔离时间的相关性进行统计分析，以探究种群遗传多样性差异的潜在原因。为探究种群的遗传结构，本研究运用MrBayes v3.2.6(Ronquistetal., 2012)构建种群间的系统发育树。碱基替代模型采用软件jModelTest v2.1.7 根据Bayesian information criterion(BIC)方法进行评估(Darribaetal., 2012)。系统发育分析采用两个独立的重复进程，每个进程运行4条MCMC链(3条热链和1条冷链)，每条链运行10 000 000代，每1 000代抽样一次，并选取前25%的样本作为老化样本舍去。当两个重复进程间的分歧频率的平均标准差(standard deviation,SD)小于0.01时，即认为MCMC链已经收敛。本研究还采用Network v4.6构建单倍型网络图(Bandeltetal., 2000)，以进一步证实系统发育结果的可靠性；并采用Arlequin v3.5进行分子方差分析(AMOVA)，评估不同谱系间的遗传分化程度(Excoffier and Lischer, 2010)。为探究舟山群岛碧蛾蜡蝉的种群历史动态，本研究还运用Arlequin v3.5进行中性检验(neutrality test)和错配分布分析(mismatch distribution analysis)。此外，本研究采用BEAST v1.8.4绘制了贝叶斯天空线(Bayesian skyline plots, BSP)，以评估种群大小随时间的变化。BSP分析以jModelTest v2.1.7估算的核苷酸替代模型为最适模型，采用松散的分子钟(uncorrelated lognormal relaxed clock)。碱基突变速率采用昆虫线粒体 COX1的进化速率(1.15% Myr-1)(Brower, 1994)，MCMC链运行50 000 000代，每5 000代抽样一次，前25%的老化样本舍去。采用Tracer v1.6确保所有参数的有效种群大小(ESS)大于200，并绘制贝叶斯天空线。

2 结果

2.1 碧蛾蜡蝉种群遗传多样性

本研究扩增获得247头碧蛾蜡蝉个体长度为740 bp的16S rRNA基因片段(GenBank登录号: MN548464-MN548710)，共检测到35个单倍型，其中共享单倍型11个(表2)。共享率最高的单倍型为H_1, 5, 7, 8, 28，均被15个以上的个体所共享。其中H_5被来自14个不同种群中的87个个体所共享，是涵盖种群数最多的一个单倍型。这表明上述几个共享单倍型可能为祖先单倍型。整体而言，该地区的碧蛾蜡蝉种群单倍型多态性(Hd)为0.817，核苷酸多态性(π)为0.008(表1)。对不同种群而言，桥头施种群(QTS)的单倍型多态性最高(Hd=0.889)，但泗礁岛种群(SJ)的核苷酸多态性最高(π=0.009)。贤氧镇种群(XY)的单倍型多态性和核苷酸多态性均为0；此外，六横峧头种群(LHJT)的单倍型多态性和核苷酸多态性在除贤氧镇外的各种群中最低。尽管不同岛屿种群的单倍型多态性和核苷酸多态性不尽相同，但统计结果显示，二者与岛屿-大陆间距离和岛屿被隔离的时间均无显著的相关性，相关系数(R2)非常小，且不显著(P>0.05)(图2)，表明二者在不同种群间的差异与岛屿-大陆间距离和岛屿被隔离的时间的相关性不显著。

表2 基于16S rRNA基因序列的舟山群岛碧蛾蜡蝉种群单倍型统计Table 2 Statistics of haplotypes of Geisha distinctissima populations in the Zhoushan Archipelagobased on 16S rRNA gene sequence

个体编码由种群代码和个体编号组成。Individual code consists of population code and individual no.

图2 基于16S rRNA基因序列的舟山群岛碧蛾蜡蝉种群遗传多样性与岛屿特征间的相关性Fig. 2 Correlation between the genetic diversity of Geisha distinctissima populations and island featuresof the Zhoushan Archipelago based on 16S rRNA gene sequenceA: 遗传多样性与岛屿到大陆距离的相关性Correlation between the genetic diversity and island distance from the mainland; B: 遗传多样性与岛屿隔离时间的相关性Correlation between the genetic diversity and the time of island isolation. 菱形为单倍型多态性，三角形为核苷酸多态性。The diamonds represent haplotype diversity, and the triangles represent nucleotide diversity.

2.2 碧蛾蜡蝉种群遗传结构

基于jModelTest估算的舟山群岛碧蛾蜡蝉16S rRNA基因的最适替代模型为F81+I，根据此模型的系统发育分析表明，尽管个别末端枝的支持度较高，但该地区的碧蛾蜡蝉种群并没有形成明显的进化支系(图3)，提示该物种在舟山群岛上可能没有形成明显的遗传结构。单倍型网络图表明，该地区种群以单倍型H_1, 5, 7, 8为核心，形成一个网状网络图(图4)，进一步证实了上述的推断。分子方差分析发现，种群内部的变异占总变异的68.26%(表3)，是种群间变异(31.74%)的2倍多，表明该物种在岛屿之间及岛屿与大陆之间并没有呈现出明显的遗传分化。

图3 舟山群岛碧蛾蜡蝉种群16S rRNA基因序列单倍型的贝叶斯系统发生树Fig. 3 Bayesian tree of 16S rRNA gene sequence haplotype of Geisha distinctissima populations in the Zhoushan Archipelago系统树节点处的数值为后验概率。The values on the note of the tree are the posterior probability.

图4 舟山群岛碧蛾蜡蝉种群16S rRNA基因序列单倍型网络图Fig. 4 Haplotype network of Geisha distinctissima populations in the Zhoushan Archipelago based on 16S rRNA gene sequence

表3 基于16S rRNA基因序列的舟山群岛碧蛾蜡蝉种群分子变异方差分析(AMOVA)Table 3 Analysis of molecular variance (AMOVA) of Geisha distinctissima populations in the Zhoushan Archipelagobased on 16S rRNA gene sequence

***显著水平小于0.0001。***The significant level is less than 0.0001.FST: 遗传分化指数Genetic differentiation index.

2.3 碧蛾蜡蝉种群历史动态

由于舟山群岛碧蛾蜡蝉种群没有呈现出明显的遗传结构，本研究将该地区所有种群合并成一个大种群，并进行种群历史动态分析。中性检验分析发现，Tajima’sD(-1.912,P=0.003)和Fu’sFs(-27.180,P<0.001)均呈显著的负值，表明碧蛾蜡蝉种群在舟山群岛及邻近的大陆上可能发生了近期的种群扩张事件。错配分布分析呈典型的单峰(图5)，且Rg. index值(Rg. index: 0.049,P=0.33)较小，SSD值(SSD=0.002,P=0.12)不显著，进一步证实该种群在舟山群岛及邻近大陆的近期扩张事件。 BSP分析表明，舟山群岛的碧蛾蜡蝉种群的扩张事件大约发生在0.5-3.5 kya，在该时期内其有效种群的大小发生了明显的增长现象(图6)。

图5 基于16S rRNA基因序列的舟山群岛碧蛾蜡蝉种群的错配分布分析Fig. 5 Mismatch distribution of Geisha distinctissima populations in the Zhoushan Archipelago based on 16S rRNA gene sequenceSSD: 平方差之和Sum of square deviations. Rg. index: 崎岖参数Raggedness index.

图6 基于16S rRNA基因序列的舟山群岛碧蛾蜡蝉种群的贝叶斯天空线Fig. 6 Bayesian skyline plots of Geisha distinctissima populations in the Zhoushan Archipelago based on 16S rRNA gene sequence实线为有效种群大小的中值，虚线为95% 的最高后验密度。Solid lines indicate the median values of the effective population size, and the dotted lines denote the 95% highest posterior density.

3 讨论

本研究基于16S rRNA基因序列对舟山群岛及邻近大陆的碧蛾蜡蝉种群遗传多样性进行了评估，以揭示岛屿隔离对其遗传多样性差异的影响。尽管碧蛾蜡蝉不同种群间的单倍型多态性和核苷酸多态性不尽相同，系统发育重建和单倍型网络图均没有解析出该地区碧蛾蜡蝉种群明显的遗传结构；且AMOVA分析表明，该地区碧蛾蜡蝉种群间没有明显的遗传分化。这些结果表明，各岛屿之间及岛屿与大陆种群间的遗传多样性没有呈现出明显的差异。本研究结果与近期基于保守的线粒体基因片段(COX1,COX2-COX3,CYTB)和核基因对该区域碧蛾蜡蝉和筛豆龟蝽Megacoptacribraria的比较谱系地理学的研究结果(Lietal., 2020)一致。然而，对整个舟山群岛上黑斑侧褶蛙的遗传多样性评估表明，冰后期海平面上升造成的岛屿隔离可阻断大陆与岛屿种群间的基因流，进而造成种群间遗传多样性的差异和遗传分化(Wangetal., 2014)。舟山群岛碧蛾蜡蝉种群遗传多样性的均质化表明，岛屿隔离效应似乎没有促使该物种遗传多样性的差异。单倍型多态性和核苷酸多态性与岛屿-大陆间距离和岛屿被隔离的时间的不相关进一步证实了上述猜想。这种种群间没有发生明显的遗传分化及遗传多样性的均质化表明，种群间可能存在明显的基因流。这一推测被多基因联合的比较谱系地理学研究得到了证实(Lietal., 2020)。人类在舟山群岛上活动的历史超过5 000年，且随着我国航海技术和造船术的不断发展与成熟，大陆与岛屿间的人类活动日趋频繁(赵利平, 2007; 方牧, 2008; 郑丽波等, 2016)，这可能是该地区碧蛾蜡蝉种群间遗传多样性均质化的原因。种群动态分析表明，舟山群岛碧蛾蜡蝉种群在大约0.5-3.5 kya发生了种群扩张。本研究估算的碧蛾蜡蝉种群扩张时间晚于多基因联合分析的种群扩张时间(Lietal., 2020)。且先前的研究认为，舟山群岛碧蛾蜡蝉种群在近5 000年内种群大小维持稳定。尽管本研究和先前的研究都使用了昆虫COX1的替代速率(0.0115 /site/Myr)作为分子钟(Lietal., 2020)，但不同的分子标记替代率不同(李可群, 2015)。因此，基于不同的分子标记估算的种群历史动态很可能不同。尽管如此，无论基于16S rRNA还是多基因的联合分析，均证实舟山群岛碧蛾蜡蝉发生了种群扩张。

本研究与先前的研究均证实，舟山群岛碧蛾蜡蝉种群间没有呈现出明显的遗传多样性的差异，且都经历了近期的种群的扩张事件。然而，不同的是，早前的多基因联合分析表明，该地区碧蛾蜡蝉种群具有明显的遗传分化，形成了3个遗传谱系，且在大约5-9 kya发生了一次种群扩张(Lietal., 2020)。本研究所获取的分子标记较单一，且片段长度较短，仅获取了35个单倍型，可能不足以很好地解析舟山群岛碧蛾蜡蝉种群的遗传分化。尽管如此，考虑到先前的研究所选分子标记的保守性及样本数的不一致性，本研究并没有将已有的分子标记与本研究获取的分子标记合并使用。分子钟是相对的，且不同分子标记进化速率不同(李可群, 2015)，可能会造成基于不同的分子标记估算的种群历史动态存在差异。由于本研究与已有的研究均存在大陆样本偏少，且研究结果均主要源于保守的线粒体基因片段，因此没有充分解析分布于舟山群岛及邻近大陆上的碧蛾蜡蝉种群的进化历史。后续的研究需广泛增加大陆样本，采用更有效的分子标记(如增加核分子标记或采用基因组数据)，以进一步深入探讨碧蛾蜡蝉在舟山群岛与大陆间的遗传分化和种群历史动态。