基于线粒体基因和核基因的中华虎凤蝶指名亚种地理种群遗传分化分析

向 颖, 董婉未, 蒋国芳,,*, 洪 芳, 张佑祥, 张文武

(1. 吉首大学生物资源与环境科学学院, 湖南吉首 416000; 2. 泉州师范学院海洋与食品学院, 福建泉州 362000;3. 南京师范大学生命科学学院, 江苏省生物多样性与生物技术重点实验室, 南京 210023;4. 湖南乌云界国家级自然保护区管理局, 湖南常德 415700)

种群是生物进化的基本单位，种群遗传分化研究是生物遗传多样性与物种进化的核心问题之一(Frankhametal., 2004)。探讨珍稀保护生物的种群遗传多样性和分化、种群历史动态等，不仅可以帮助我们理解种群形成及其演化过程，外界因素对该种群遗传特性的影响以及种群的生态适应性(Slatkin, 1987)，还可以为扩大其种群数量及相关保护策略的制定提供理论依据。

随着分子生物学的发展和相关技术的实践应用，种群分化研究所应用的分子标记也越来越多，如线粒体、核糖体和微卫星等(Sei and Porter, 2007; Kodandaramaiahetal., 2012)。mtDNA已成为谱系地理学研究中最常用的分子标记(张代臻等, 2008)，而核基因则能够补充mtDNA不具备的信息。所以，将两类不同进化速率和遗传方式的分子标记(线粒体基因和核基因)相结合进行种群遗传结构分析，结果将更加全面可靠。近年来，已有学者利用这两类分子标记对蝴蝶种群遗传多样性(Bulletal., 2006)、种群遗传结构(Czajkowskaetal., 2020)、种群历史动态(Flanaganetal., 2004)和生物地理格局及其形成历史(Dincăetal., 2011)等进行了深入探讨。

中华虎凤蝶Luehdorfiachinensis隶属鳞翅目(Lepidoptera)凤蝶科(Papilionidae)虎凤蝶属Luehdorfia。该蝶仅分布于我国秦岭和长江中下游部分地区(董思雨等, 2014)，是中国特有的蝴蝶，属国家二级保护野生动物，极具生态保护与观赏价值。目前关于中华虎凤蝶的分布范围(何桂强等, 2011; 李密等, 2012)、生物学特性(童雪松和潜祖琪, 1992; 王文明等, 2010)、保护现状(袁德成等, 1998)和食性研究(姚洪渭等, 1999; 姚肖永等, 2008)等方面均有报道，但中华虎凤蝶种群遗传分化研究的相关报道缺乏。本研究首次采用4个分子遗传标记，即线粒体COI基因和3个非编码核基因(ITS2,Pgi-intron和Tpi-intron)片段，研究我国5个地区(江苏、浙江、安徽、江西和湖南)的中华虎凤蝶指名亚种L.chinensischinensis地理种群间的遗传分化，以期为中华虎凤蝶的种群保护提供重要科学依据。

1 材料与方法

1.1 供试虫源

本研究所用中华虎凤蝶指名亚种60头幼虫样本于2013-2014年间分别采集自江西、浙江、安徽、江苏和湖南5个地区。收集到的样本均浸泡在无水乙醇中保存备用。详细信息见表1。

表1 中华虎凤蝶指名亚种样本采集信息

1.2 PCR扩增及测序

使用Progema试剂盒分别提取中华虎凤蝶指名亚种各地理种群各头幼虫样本的总DNA，经1%的琼脂糖凝胶电泳检测，保存于-20℃冰箱中备用。本研究选取了线粒体基因COI以及非编码核基因ITS2,Pgi-intron和Tpi-intron作为目的基因。线粒体基因COI的引物是依据中华虎凤蝶线粒体全基因组序列(Liuetal., 2013)设计；而核基因引物是根据EPIC(exon-primed-intron-crossing)方法(Jennings and Etter, 2011)设计。引物均由南京金斯瑞生物科技有限公司合成，引物的详细信息见表2。PCR扩增反应体系(25 μL)：模板DNA 1 μL，上下游引物(10 μmol/L)各0.5 μL， 10×PCR Buffer 2.5 μL， MgCl2(20 μmol/L) 2 μL， dNTPs (10 μmol /L) 2 μL， Taq 酶(2 U/μL) 0.2 μL， ddH2O 16.3 μL。PCR扩增程序： 95℃预变性5 min； 95℃变性30 s， 55℃(核基因)或51℃(COI)退火30 s， 35个循环； 72℃延伸1 min(核基因)或30 s(COI)；循环结束后72℃延伸补齐10 min， 4℃保温。所有PCR产物送至南京金斯瑞生物科技有限公司进行测序。

表2 引物信息

1.3 数据分析

利用DNASTAR (Burland, 2000)和MEGA 7.0 (Kumaretal., 2016)对扩增基因进行校正比对，分析序列碱基组成及序列变异情况。使用SequenceMatrix软件将4个分子标记串联后进行以下分析：使用DnaSP 5.10.1 (Librado and Rozas, 2009)计算所有采集到的5个种群的多态性位点数(S)、核苷酸多样性(π)和平均核苷酸差异数(K)。使用软件Arlequin v3.5.1.3 (Excoffier and Lischer, 2010)进行种群结构和遗传变异分析，计算各地理种群的Tajima’sD值(Tajima, 1989)、Fu’sFs值(Fu, 1997)、基因流(Nm)及群体间的分化指数(Fst)。群体间Fst系数和基因流(Nm)分别是衡量种群间遗传分化程度和基因交流程度大小的指数(Hudsonetal., 1992; Rousset, 1997)，Fst>0.25和Nm<1，则种群遗传分化大，基因流小。

为检验种群间的遗传分化与种群间的地理距离是否具有相关性，在MEGA 7.0中利用Kimura 2-Parameter (K2P)模型分别计算地理种群的遗传距离。根据经纬度计算各采样地点的地理距离及其自然对数，利用PC-ORD 5软件(Grandin, 2006)对遗传距离及地理距离的自然对数进行Mantel检验，分析地理种群间遗传距离和地理距离的相关性。

2 结果

2.1 碱基组成及序列变异

对5个中华虎凤蝶指名亚种地理种群共60头供试样本分别进行了COI和ITS2,Pgi-intron和Tpi-intron4个基因片段扩增测序，每个基因片段各获得了60条序列。各基因片段的GenBank登录号见表1，序列特征和变异信息见表3。COI基因的变异位点和简约信息位点最少，A, G, C和T平均含量分别为30.62%, 13.72%, 14.74%和40.92%，符合昆虫线粒体基因序列组成(Liu and Beckenbach, 1992)。使用3对核基因序列扩增引物对中华虎凤蝶指名亚种样本DNA进行扩增，获得的基因序列长度分别为760 bp (ITS2), 1 122 bp (Pgi-intron)和757 bp (Tpi-intron)，Tpi-intron变异位点和简约信息位点最多，其次是Pgi-intron和ITS2。4个基因序列串联后的基因序列(以下简称“合并序列”)全长为3 277 bp，A+T%(66.6%)>G+C%(33.4%)，和COI基因一样，有较强的AT偏好性。

表3 中华虎凤蝶指名亚种地理种群线粒体基因和核基因序列的碱基组成及变异

2.2 遗传多样性

基于合并序列对中华虎凤蝶指名亚种5个地理种群进行遗传多样性分析，结果见表4。总群体核苷酸多样性(π)为0.00686，平均核苷酸差异数(K)为18.289。各地理种群内核苷酸多样性(π)为0.00121～0.01138，多态性位点数(S)为7～92，其中宝华山种群(BHS)的π和S均为最小，而桃红岭种群(THL)的均为最大。

表4 基于线粒体基因和核基因合并序列的中华虎凤蝶指名亚种地理种群的遗传多样性和中性检验

2.3 遗传分化

基于合并序列的种群间分化指数和基因流如表5所示，总群体Fst为0.4910，Nm为0.2600，说明地理种群间的遗传分化程度较大。宝华山种群(BHS)和乌云界种群(WYJ)Fst最大(0.7819)，Nm最小(0.1394)，两地理种群间分化程度高，基因交流不频繁；而宝华山种群(BHS)和滁州种群(CZ)Fst最小(0.0812)，Nm最大(5.6582)，说明它们基因交流程度高，遗传分化较小。对中华虎凤蝶指名亚种地理种群进行AMOVA分析(表6)，结果表明中华虎凤蝶指明亚种的变异主要是来自于种群内部(52.79%)，群体间的遗传变异略小。

表5 基于线粒体基因和核基因合并序列的中华虎凤蝶指名亚种地理种群的配对种群间分化指数(Fst)(下三角)和基因流(Nm)(上三角)

依据Kimura 2-Parameter双参数模型，基于合并序列计算中华虎凤蝶指名亚种地理种群间的遗传距离(表7)。各地理种群的遗传距离在0.0019～0.0105之间，遗传距离最大的是杭州种群(HZ)和桃红岭种群(THL)(0.0105)，最小的为滁州种群(CZ)和宝华山种群(BHS) (0.0019)。就地理距离来讲，相距最近的是滁州种群(CZ)和宝华山种群(BHS)，相距最远的为乌云界种群(WYJ)和杭州种群(HZ)。对地理种群间的地理距离的对数值和遗传距离进行Mantel相关性检验，结果显示R=0.8519,P=0.018<0.05，表明中华虎凤蝶指名亚种各地理种群间的地理距离和遗传距离呈正相关关系。

表6 基于线粒体基因和核基因合并序列的中华虎凤蝶指名亚种地理种群的分子方差分析(AMOVA)

表7 基于线粒体基因和核基因合并序列的中华虎凤蝶指名亚种地理种群间的成对遗传距离(下三角)与地理距离(km)自然对数值(上三角)

2.4 种群历史动态

基于合并序列对中华虎凤蝶指名亚种地理种群中性检验和错配分布进行分析(图1和表4)，总种群的Tajima’sD和Fu’sFs分别为2.52028和-0.32628，但不显著，而且总种群的错配分布图呈现多峰分布。对各种群来讲，Tajima’sD和Fu’sFs均为正值，这表明近期各种群都没有出现扩张。

3 讨论

核苷酸多样性(π)是评价物种遗传多样性的重要参数(Grant and Bowen, 1998)，基因联合片段序列π>0.005，表明中华虎凤蝶指名亚种总种群遗传多样性较高。合并基因序列的AMOVA分析结果表明遗传分化主要来自种群内部，呼晓庆和杨兆富(2019)对中国草地螟Loxostegesticticalis地理种群的遗传分析中也得出了类似结果。

图1 基于线粒体基因和核基因合并序列的中华虎凤蝶指名亚种地理种群的错配分布

有研究表明，对植食性昆虫而言，地理隔离和寄主植物分化是使不同种群产生遗传分化重要原因(Forbesetal., 2017; 朱诗禹等, 2017; Goodmanetal., 2019)。中华虎凤蝶指名亚种总种群Fst>0.4,Nm<1，表明种群间的遗传分化大，基因交流少。Mantel检验结果表明遗传距离与地理距离呈现正相关性，而许多其他研究也得出类似结论(薛进等, 2007; 谢艳兰等, 2019)。值得注意的是WYJ种群和BHS, CZ及HZ种群的遗传分化非常大(Fst>0.65,Nm<1)，我们认为这不仅仅和地理隔离有关，寄主植物的分化也产生了重要影响；BHS, CZ及HZ种群属于低海拔类群，寄主植物为杜蘅Asarumforbesii(袁德成等, 1998; 王文明等, 2010)，而WYJ种群是高海拔种群，寄主植物为细辛A.sieboldii(谭济才等, 2010)，中华虎凤蝶不同地理种群对不同寄主植物产生了适应性进化。

Tajima’sD和Fu’sFs检验是可反映出种群变化动态历史的中性检验方法。Fu’sFs和Tajima’sD检验值<0时，则表明种群有扩张及瓶颈效应，反之则种群稳定(Harpendingetal., 1998; Ramos-Onsins and Rozas, 2002)。中华虎凤蝶指名亚种各种群Tajima’sD值均大于0，Fu’sFs值不显著，由此推测，中华虎凤蝶指名亚种在历史上并没有发生种群扩张。错配分布是以可视化方式体现群体历史动态，当种群形成历史较短或发生扩张时，其错配分布曲线呈现单峰；而当种群为存在时间较长且稳定的种群时，错配分布曲线则呈现多峰曲线分布(Excoffier and Lischer, 2010; 苏军虎等, 2017)。中华虎凤蝶指名亚种总种群错配分布图呈现多峰分布，其种群在历史上较稳定。综上，中华虎凤蝶指名亚种各地理种群较稳定，近期没有经历明显的种群扩张。

本研究采用4个分子标记对中华虎凤蝶指名亚种5个地理种群进行研究，结果显示各地理种群遗传多样性较高，遗传变异较丰富，物种的环境适应能力较强，较容易在新的环境中生存，表明近年来中华虎凤蝶的保护工作是有一定成效的。但值得注意的是，中华虎凤蝶种群分布格局具有密度低、种群间隔大的特征(袁德成等, 1998; 董思雨等, 2014)。因此，在后续的保护工作中不能掉以轻心，应继续保护好其栖息地和寄主植物资源，严禁采挖；对于种群数量小的野外种群也可以输入人工繁育个体保证其有效种群大小。

致谢本研究所用部分蝴蝶幼虫由江西中医学院贾凤海教授、滁州学院褚立新教授和浙江杭州楼信权先生提供，在此深表谢意！