日本鳗鲡离散SNP标记筛选及群体遗传多样性分析❋

于 磊， 刘炳舰， 刘进贤

(1. 中国科学院海洋研究所海洋生态与环境科学重点实验室，山东 青岛 266071； 2. 中国科学院大学，北京 100049；3. 青岛海洋科学与技术国家实验室海洋生态与环境科学功能实验室，山东 青岛 266071)

生物与其所生存的环境密不可分。在本地环境的选择压力下，本地群体会进化出与其生存环境相适应的表型特征，更好地适应本地环境，这就是“本地适应性”(Local adaptation)[1]。阐明本地适应性发生的程度和动态，能够加深对于生物多样性的认识，同时也可以为生物资源的保护提供参考[2]。

早期本地适应性研究多采用候选基因法。研究者从已知生物学功能的编码蛋白质的基因入手，通过群体遗传学分析检测本地适应性的信号[3-6]。该方法能够成功的关键在于目标候选基因的选择，然而准确地选择合适的基因在实际操作中难度极大。除此之外，候选基因法也很难应用到非模式生物的研究中。由于技术方法的束缚，关于生物种群本地适应性的研究进展缓慢。近年来，随着二代测序技术(Next-generation sequencing technology)的迅猛发展，基于全基因组扫描的方法被应用到生物种群本地适应性的研究中。该方法可以开发大量单核苷酸多态性标记(Single nucleotide polymorphism, SNP)，大大促进了该领域的发展。SNP标记主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。单核苷酸多态性标记具有在基因组中数量多，分布密度高等特点，而且在基因分型过程中不需要根据片段大小将DNA分带，可以大规模地高度自动化地完成，相比微卫星标记工作效率更高，更适于大样本量检测分析，SNP标记的使用是群体遗传学研究的一大突破[7]。研究者可以使用整个基因组范围内大量的SNP标记进行群体遗传学分析，对种群遗传学参数的解析力度骤增(例如，基因流和有效群体大小)。除此之外，基于全基因组扫描的方法在非模式生物的研究中体现出基于候选基因的方法不可比拟的优势。目前，基于二代测序技术检测生物种群本地适应性的研究开发了大量离散遗传位点(Outlier markers)，已经在多种海洋生物中检测到了适应性分化的信号，例如：大西洋鳕鱼(Gadusmorhua)[8-9]、大西洋鲱鱼(Clupeaharengus)[10-11]和太平洋七鳃鳗(Entosphenustridentatus)[12]等。然而，对于存在中性遗传分化的物种，在离散遗传位点上检测到的生物种群的遗传结构是中性力量和自然选择力量共同作用的结果，很难准确解读适应性分化的模式[13]。

日本鳗鲡(Anguillajaponica)隶属于鳗鲡目(Anguilliformes)鳗鲡科(Anguillidae)鳗鲡属(Anguilla)，是一种长距离降海洄游性鱼类[14]。成鱼广泛分布于东亚河口和淡水河流，在日本、中国，以及朝鲜半岛都有分布[15]。日本鳗鲡的产卵场位于马里亚纳群岛的西部海域(大致为15°N, 140°E)[16]，孵化后的柳叶鳗(Leptocephalus)随海流漂流至东亚各国沿海，进而抵达河口水域，变态为玻璃鳗(Glass eel)[17]。而后经过复杂的变态发育，成体在淡水河流中生活。从1980年代开始，日本鳗鲡的种群衰退不断加剧。玻璃鳗(Glass eel)、黄鳗(Yellow eel)以及银鳗(Silver eel)的野生群体资源量均出现了严重衰退[18]。在2016年公布的IUCN濒危物种红色名录中，日本鳗鲡已经被列为“濒危”[19]。虽然目前已经可以人工育成日本鳗鲡的玻璃鳗幼体[20]，但是由于尚未攻克玻璃鳗幼体饵料的难题，因此人工育成的幼苗成活率低，活性弱，尚无法实现产业化生产，日本鳗鲡的养殖仍然主要依靠从自然海域捕获野生幼苗。随着日本鳗鲡野生资源的衰退，鳗鲡市场受到巨大的冲击。目前研究表明，水利设施建设、过度捕捞、环境栖息地的破坏，以及全球气候变化等都可能是造成鳗鲡种群衰退的原因，然而这些因素对鳗鲡资源量影响的强度和机制尚不明确[21]。日本鳗鲡作为随机交配物种[22-25]的典型代表，各地理群体均为同一基因库随机分配，并不存在中性遗传分化。因此在离散位点上检测到的群体间遗传分化的信号是其单一世代内受自然选择作用的结果。这使其成为研究本地适应性问题的理想材料。

在本研究中，我们利用酶切位点相关DNA测序(Restriction site-associated DNA sequencing, RAD-Seq)技术开发日本鳗鲡多态性离散SNP标记，并通过KASPar(KBiosciences Competitive Allele-Specific PCR genotyping system)技术进行SNP位点的多态性验证及在两个群体中开展初步群体遗传学分析。相关SNP标记可以用于未来日本鳗鲡群体遗传结构的解析以及种群适应本地环境机制的研究。

1 材料与方法

1.1 样品采集和和DNA提取

在福建宁德市三沙镇(SS)和辽宁丹东市(DD)向当地从事鳗苗捕捞作业的渔民购买日本鳗鲡的玻璃鳗样品(SS，2014年3月，36尾；DD，2014年5月，30尾)，保存于95%的酒精中。DNA提取采用传统方法。首先用蛋白酶K消化肌肉组织，然后使用苯酚-氯仿法抽提DNA。

1.2 SNP开发

从两个群体各选取12尾玻璃鳗样品进行RAD-seq。文库构建和RAD-seq的具体实施参照了张柏东等[26]发表的论文。每个样本RAD-seq获得的总读长至少达到1×基因组(约1.2G)深度。原始数据按照如下标准进行过滤：1.去除有接头污染的序列；2.去除未识别碱基大于或者等于10%的序列；3.去除phred值小于5的碱基占总序列长度大于或者等于50%的序列；4.去除文库构建过程中PCR扩增产生的重复序列；5.去除序列上没有EcoRI酶切位点(AATTC)的序列。将过滤后的序列根据个体特异的标签序列进行归类。使用BWA 0.5.9将序列比对到日本鳗鲡的基因组草图上[27]，参数使用默认值(Mismatch penalty 4；Gap open penalty 6)[28]。使用SAMTOOLS[29]获取SNP位点的信息。

使用ARLEQUIN 3.5[30]筛选可能受到自然选择作用的位点。该软件可以绘制遗传距离FST值关于杂合度的散点图，将其与模拟的中性分布进行对比来进行离散位点的筛选。中性分布使用双岛屿模型(Two-island model)进行模拟，假设自由交配[31]。

1.3 多态性验证

SNP分型由艾吉析科技(上海)有限公司完成。该公司使用KASPar技术对丹东群体的30尾玻璃鳗样品和三沙群体的31尾玻璃鳗样品进行SNP分型。KASPar技术基于竞争性位点特异性PCR扩增技术，将正向引物的3’末端设计在SNP位点上，同一位点的两个不同等位基因被扩增为携带不同荧光颜色的产物。每个SNP位点所使用的特异性引物信息见表 1。

使用ARLEQUIN 3.5计算观测杂合度(Ho)，期望杂合度(He)。使用GENEPOP 4.5.1[32]检验哈温平衡和连锁不平衡，显著性水平用Bonferroni方法进行校正。使用ARLEQUIN 3.5计算两群体间的遗传距离FST值，显著性水平通过10 000次重抽样模拟(Bootstrapping permutations)进行计算。

2 结果

2.1 RAD-seq过滤结果

在RAD-seq过程中，最终产出42.1Gb原始数据，MS25个体的原始数据量达到6×基因组(约1.2G)深度，其余个体的数据量达到1×基因组(约1.2G)深度，Q20比例基本都在90%以上。经过数据过滤后，共产出40.0Gb高质量数据，平均每个样本1.67Gb数据。数据质量较高，数据量较大，满足后续分析需要。

2.2 SNP开发和验证

基于RAD-seq的数据共获得73 557个SNP位点的信息。ARLEQUIN离散位点检验中，使用R软件处理ARLEQUIN软件输出的文本文件，选取FST值达到99.5%分位数以上的位点作为候选位点，共有250个。在这些位点中，有85个SNP位点可以成功地进行功能注释，涉及到多个方面的生物学功能。由于能成功进行功能注释的位点更有可能在生物的本地适应性中发挥作用，因此使用KASPar技术对其中48个能进行功能注释的位点进行SNP分型。有6个位点未分型成功，2个位点的数据缺失率大于15%，最终有40个位点的数据可用。在本研究中，同一位点2个不同等位基因扩增后的产物分别被HEX和FAM标记，根据2种荧光信号的强度进行SNP分型，荧光信号强度较强，分型结果清晰。

2.3 遗传标记信息

在这40个SNP位点中，2个位点具有单态性，AJ_SNP_06位点在丹东群体和三沙群体中均只检测到胞嘧啶脱氧核糖核苷酸(C)一种等位基因，而AJ_SNP_08只检测到腺嘌呤脱氧核糖核苷酸(A)一种等位基因。其余38个位点具有多态性，在丹东(DD)群体中，这38个位点的观测杂合度(Ho)为0～0.407，期望杂合度(He)为0.033～0.484；在三沙(SS)群体中，这38个位点的观测杂合度(Ho)为0.032～0.500，期望杂合度(He)为0.032～0.494。具体信息见表 2。经过Bonferroni方法校正，在丹东(DD)群体中，AJ_SNP_03位点偏离哈温平衡，所有的位点间均不存在连锁不平衡现象。在三沙(SS)群体中，AJ_SNP_25位点和AJ_SNP_35位点偏离哈温平衡，所有位点间均不存在连锁不平衡现象。综上，这38个SNP标记具有较高的多态性，且未检测到偏离哈温平衡或出现连锁不平衡，因此符合群体遗传学分析的需要，可用于解析日本鳗鲡的群体遗传结构。而且这些标记均可定位到日本鳗鲡基因组草图上，研究者可根据自己的需要选择SNP分型方法。丹东群体和三沙群体间的遗传距离FST值为-0.005(P=0.898)，表明这2个群体间并无显著的遗传分化。

3 讨论

遗传多样性是生物与环境长期相互作用的结果，也是生物适应生存环境的基础。一般认为，遗传多样性水平越高，生物应对环境变化的能力越强，反之亦然[33]。生物遗传多样性的水平通过分子标记来评估，目前微卫星标记和SNP标记是最流行的2种分子标记。与微卫星标记相比，SNP标记在分型的准确性以及数据处理的自动化等方面更具优势[34]。本研究开发的38个多态性SNP位点，在丹东(DD)群体中的观测杂合度(Ho)为0～0.407，期望杂合度(He)为0.033～0.484；在三沙(SS)群体中的观测杂合度(Ho)为0.032～0.500，期望杂合度(He)为0.032～0.494。日本鳗鲡的这些SNP标记的杂合度并未出现显著降低，与其他海洋鱼类的杂合度水平相当。在隆头鱼(Ctenolabrusrupestris)中，SNP标记的期望杂合度(He)为0.063～0.495[35]；在大黄鱼(Larimichthyscrocea)中，SNP标记的期望杂合度(He)为0.066～0.503[36]；在欧洲川鲽(Platichthysflesus)中，SNP标记的期望杂合度(He)为0.176～0.600[37]。结果表明，尽管日本鳗鲡的种群资源严重衰退，然而群体遗传多样性水平并未出现显著降低。但是这并不意味着日本鳗鲡的资源状况良好，而有可能与海洋鱼类的生活史策略(Life history strategy)相关。日本鳗鲡为r-对策者，繁殖大量子代，但子代存活率低。因此即使种群遭遇瓶颈，子代可能依然会保持较高的遗传多样性水平。然而，日本鳗鲡的资源衰退已十分严重，若不加强对日本鳗鲡种群资源的保护，其资源状况可能陷入恶性循环，最终导致日本鳗鲡物种的灭绝。

探究生物不同地理群体的遗传结构，不仅可以加深对生物多样性的认识，也是合理有效地管理生物资源的基础。目前对海洋鱼类群体遗传结构的研究大多采用中性遗传标记，这些研究的结果表明，很多海洋鱼类不同地理群体间并不存在显著的遗传分化。使用离散位点标记进行研究，可以揭示更精细的群体遗传结构[38-39]。本研究开发的SNP标记是根据丹东群体和三沙群体RAD-seq的数据开发的两群体间高FST值标记，因此在解析日本鳗鲡群体遗传结构时会具有更高的分辨率。在本研究中，基于38个SNP标记的丹东群体和三沙群体间的遗传距离FST值为-0.005(P= 0.898)，且不显著，表明这2个群体间并无显著的遗传分化。这可能与日本鳗鲡随机交配[22-25]的生物学特征有关。日本鳗鲡性成熟后会洄游到共同的产卵场进行繁殖，不同地理群体是从同一基因库随机分配而来，因此遗传组成并无显著差异。然而，本研究开发的SNP标记是两群体间高FST值的标记，应当具有更高的解析力。这2个群体间之所以出现小而不显著的FST值，可能是由于开发SNP时采用的样本量较小。开发SNP标记时选取的是小样本，丹东群体和三沙群体各12尾，而后采用大样本进行SNP分型。小样本检测到的离散位点可能包含统计学误差，导致同样的SNP标记在大样本中开展群体遗传学分析时不能检测到显著的遗传分化。此外，由于SNP标记只有2个等位基因，需通过使用大量标记来提高其对群体遗传结构的解析力，因此本研究中分析的SNP位点较少也可能是检测到的遗传分化水平低的原因。

除了用于研究群体遗传结构，离散位点标记对于本地适应性研究也十分重要。日本鳗鲡作为随机交配物种[22-25]的典型代表，各地理群体均为同一基因库随机分配，并不存在中性遗传分化。因此在离散位点上检测到的群体间分化的信号为其单一世代内受自然选择作用的结果。这使其成为研究本地适应性问题的理想材料。目前已有研究使用SNP标记对鳗鲡属鱼类的本地适应性开展研究[40-42]，然而SNP标记资源在日本鳗鲡中尚属空白。本研究首次通过二代测序技术开发了一批离散SNP标记，这些SNP标记是基于丹东群体和三沙群体筛选的高FST值的分子标记，有可能与参与日本鳗鲡对本地环境的适应，因此可以作为候选分子标记应用到日本鳗鲡的本地适应性研究中，对深入理解微时间尺度内生物的适应性进化问题具有重要意义。

4 结语

本研究首次通过二代测序技术开发了一批离散SNP标记。这些标记大多符合哈温平衡，且标记间不存在连锁不平衡现象，符合群体遗传学分析的要求，可用于解析日本鳗鲡的群体遗传结构。另外，这些SNP标记是基于丹东群体和三沙群体筛选的高FST值的分子标记，极有可能与参与日本鳗鲡对本地环境的适应，因此可以作为候选分子标记应用到日本鳗鲡的本地适应性研究中。

致谢：感谢李玉龙在数据处理中给予的帮助。

[1] Williams G C. Adaptation and Natural Selection[M]. New Jersey: Princeton University Press, 1966: 20-55.

[2] Hauser L, Carvalho G R. Paradigm shifts in marine fisheries genetics: ugly hypotheses slain by beautiful facts[J]. Fish Fish, 2008, 9: 333-362.

[3] Hemmer-Hansen J, Nielsen E E, Frydenberg J, et al. Adaptive divergence in a high gene flow environment:Hsc70 variation in the European flounder (PlatichthysflesusL.)[J]. Heredity, 2007, 99: 592-600.

[4] Gaggiotti O E, Bekkevold D, Jorgensen H B H, et al. Disentangling the effects of evolutionary, demographic, and environmental factors influencing genetic structure of natural populations: Atlantic herring as a case study[J]. Evolution, 2009, 63: 2939-2951.

[5] Larmuseau M H D, Vancampenhout K, Raeymaekers J A M, et al. Differential modes of selection on therhodopsingene in coastal Baltic and North Sea populations of the sand goby,Pomatoschistusminutus[J]. Mol Ecol, 2010, 19: 2256-2268.

[6] Andersen O, De Rosa MC, Pirolli D, et al. Polymorphism, selection and tandem duplication of transferrin genes in Atlantic cod (Gadusmorhua) - Conserved synteny between fish monolobal and tetrapod bilobal transferrin loci[J]. Bmc Genet, 2011, 12: 51. doi: 10. 1186/1471-2156-12-51.

[7] Allendorf. Genetics and the conservation of natural populations: allozymes to genomes[J]. Mol Ecol, 2017, 26: 420-430.

[8] Moen T, Hayes B, Nilsen F, et al. Identification and characterisation of novel SNP markers in Atlantic cod: Evidence for directional selection[J]. Bmc Genet, 2008, 9: 18. doi: 10. 1186/1471-2156-9-18.

[9] Nielsen E E, Hemmer-Hansen J, Poulsen N A, et al. Genomic signatures of local directional selection in a high gene flow marine organism; the Atlantic cod (Gadusmorhua)[J]. Bmc Evol Biol, 2009, 9: 276. doi: 10. 1186/1471-2148-9-276.

[10] Lamichhaney S, Barrio A M, Rafati N, et al. Population-scale sequencing reveals genetic differentiation due to local adaptation in Atlantic herring[J]. P Natl Acad Sci USA, 2012, 109: 19345-19350.

[11] Limborg M T, Helyar S J, de Bruyn M, et al. Environmental selection on transcriptome-derived SNPs in a high gene flow marine fish, the Atlantic herring (Clupeaharengus)[J]. Mol Ecol, 2012, 21: 3686-3703.

[12] Hess J E, Campbell N R, Close D A, et al. Population genomics of Pacific lamprey: adaptive variation in a highly dispersive species[J]. Mol Ecol, 2013, 22: 2898-2916.

[13] Flatt T. Genomics of clinal variation inDrosophila: disentangling the interactions of selection and demography[J]. Mol Ecol, 2016, 25: 1023-1026.

[14] 张春光, 唐文乔，刘东，等.中国动物志(硬骨鱼纲/鳗鲡目、背棘鱼目)[M].北京：科学出版社，2010. Zhang C G, Tang W Q, Liu D, et al. Fauna Sinica (Osteichthyes: Anguilliformes, Notacanthiformes)[M]. Beijing: Science Press, 2010.

[15] Tesch F W. The Eel: Biology and Management of Anquillid Eels[M]. London: Chapman & Hall, 1977.

[16] Tsukamoto K. Discovery of the spawning area for Japanese eel[J]. Nature, 1992, 356: 789-791.

[17] Tseng M C, Tzeng W N, Lee S C. Population genetic structure of the Japanese eelAnguillajaponicain the northwest Pacific Ocean: evidence of non-panmictic populations[J]. Mar Ecol Prog Ser, 2006, 308: 221-230.

[18] Tatsukawa. Eel Resources in East Aisa. In: Yamauchi K (ed) Eel Biology[M]. 1st edn. Tokyo: Springer-Verlag, 2003: 293-298.

[19] Jacoby D,Gollock M.The IUCN Red List of Threatened species[EB/OL].(2016-03-01)[2017-04-25].http://www.incnredlist.org/details/166184/0: IUCN,2014.

[20] Nomura K, Ozaki A, Morishima K, et al. A genetic linkage map of the Japanese eel (Anguillajaponica) based on AFLP and microsatellite markers[J]. Aquaculture, 2011, 310: 329-342.

[21] Dekker W. Status of European eel stock and fisheries. In: Yamauchi K (ed) Eel Biology[M]. Tokyo: Springer Verlag, 2003: 237-254.

[22] Sang T K, Chang H Y, Chen C T, et al. Population-structure of the Japanese eel,Anguilla-japonica[J]. Mol Biol Evol, 1994, 11: 250-260.

[23] Ishikawa S, Aoyama J, Tsukamoto K, et al. Population structure of the Japanese eelAnguillajaponicaas examined by mitochondrial DNA sequencing[J]. Fisheries Sci, 2001, 67: 246-253.

[24] Han Y S, Hung C L, Liao Y F, et al. Population genetic structure of the Japanese eelAnguillajaponica: panmixia at spatial and temporal scales[J]. Mar Ecol Prog Ser, 2010, 401: 221-232.

[25] Minegishi Y, Aoyama J, Yoshizawa N, et al. Lack of genetic heterogeneity in the Japanese eel based on spatiotemporal sampling[J]. Coast Mar Sci, 2012, 35: 269-276.

[26] Zhang B D, Xue D X, Wang J, et al. Development and preliminary evaluation of a genomewide single nucleotide polymorphisms resource generated by RAD-seq for the small yellow croaker (Larimichthyspolyactis)[J]. Mol Ecol Resour, 2016, 16: 755-768.

[27] Henkel C V, Dirks R P, de Wijze D L, et al. First draft genome sequence of the Japanese eel,Anguillajaponica[J]. Gene, 2012, 511: 195-201.

[28] Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform[J]. Bioinformatics, 2010, 26: 589-595.

[29] Li H, Handsaker B, Wysoker A, et al. The Sequence Alignment/Map format and SAMtools[J]. Bioinformatics, 2009, 25: 2078-2079.

[30] Excoffier L, Lischer H E L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[J]. Mol Ecol Resour, 2010, 10: 564-567.

[31] Beaumont M A, Nichols R A. Evaluating loci for use in the genetic analysis of population structure[J]. P Roy Soc B-Biol Sci, 1996, 263: 1619-1626.

[32] Rousset F. GENEPOP′007: a complete re-implementation of the GENEPOP software for Windows and Linux[J]. Mol Ecol Resour, 2008, 8: 103-106.

[33] Frankham R, Briscoe D A, Ballou J D. Introduction to Conservation Genetics[M]. Cambridge: Cambridge University Press, 2002.

[34] Vignal A, Milan D, Sancristobal M, et al. A review on SNP and other types of molecular markers and their use in animal genetics[J]. Genet Sel Evol, 2002, 34: 275-305.

[35] Jansson E, Taggart J B, Wehner S, et al. Development of SNP and microsatellite markers for goldsinny wrasse (Ctenolabrusrupestris) from ddRAD sequencing data[J]. Conserv Genet Resour, 2016, 8: 201-206.

[36] Wang P P, Xiao S J, Han Z F, et al. SNP discovery in large yellow croaker (Larimichthyscrocea) using Roche 454 pyrosequencing sequencing platform[J]. Conserv Genet Resour, 2015, 7: 777-779.

[37] Pedron N, Morvezen R, Moan A L, et al. New set of candidate gene SNPs and microsatellites to disentangle selective and neutral processes shaping population responses of European flounder (Platichthysflesus) to anthropogenic stress and contrasted environments[J]. Conserv Genet Resour, 2015, 7: 823-826.

[38] Russello M A, Kirk S L, Frazer K K, et al. Detection of outlier loci and their utility for fisheries management[J]. Evol Appl, 2012, 5: 39-52.

[39] Milano I, Babbucci M, Cariani A, et al. Outlier SNP markers reveal fine-scale genetic structuring across European hake populations (Merlucciusmerluccius)[J]. Mol Ecol, 2014, 23: 118-135.

[40] Gagnaire P A, Normandeau E, Cote C, et al. The genetic consequences of spatially varying selection in the panmictic American eel (Anguillarostrata)[J]. Genetics, 2012, 190: 725-736.

[41] Ulrik M G, Pujolar J M, Ferchaud A L, et al. Do North Atlantic eels show parallel patterns of spatially varying selection?[J]. Bmc Evol Biol, 2014， 14： 138.doi:10.1186/1471-2148-14-138.

[42] Pujolar J M, Jacobsen M W, Als T D, et al. Genome-wide single-generation signatures of local selection in the panmictic European eel[J]. Mol Ecol, 2014, 23: 2514-2528.