水仙(Narcissus tazetta var. chinensis)为石蒜科(Amaryllidaceae)多年生鳞茎植物，是欧洲水仙(N. tazetta)的变种，也是世界著名观赏花卉和中国传统十大名花之一，以其独特的形态和香气为世人所熟知。史籍考证与植物学证据表明该变种系地中海沿岸起源，主要归化分布于亚洲东部沿海地区或近海岛屿^[1]，在我国栽培历史悠久，具有重要的观赏、经济与文化价值。随着对水仙种质资源数量与质量需求的日益提高，前人已在抗性生理^[2-3]、矮化机理^[4]、组培扩繁^[5]、功能基因^[6]等方面取得了一系列的研究。

掌握植物遗传多样性格局，是开展资源保护与利用的前提基础。探究植物群体遗传结构在时间和空间上的遗传变异特征及分布格局，在反映物种进化历史、生境适应等方面具有重要作用^[7-9]。群体研究主要包括传统的形态学和现代分子标记2种手段, 分子标记具有不受环境或植株发育阶段影响、规避测量误差、可重复性高等优势。近10年来，随着分子生物学发展，国内外有关水仙属遗传多样性的研究，已由早期的AFLP、ISSR、RAPD等标记类型^[10-11]逐渐过渡到DNA测序。其中，基于叶绿体DNA (cpDNA)标记和序列分析具有不受选择影响、多态性丰富等优点；同时，大多数被子植物cpDNA为母系遗传，比双亲遗传的nrDNA在亚种群间表现出更大的遗传分化，能产生更清晰的种群历史踪迹^[12]，因此，当前cpDNA标记已成为研究植物群体遗传与谱系地理的主流手段之一，广泛应用于植物种间较低分类阶元和种内谱系地理学研究。

在野外调查中，笔者发现当前我国水仙归化地普遍存在生境片断化严重、种群规模和数量锐减等问题，有关其在区域尺度上的遗传多样性现状尚未报道。鉴于此，本研究采用2个水仙高变区的cpDNA标记，首次对我国水仙集中分布的东部沿海地区的5个主要归化群体的序列变异进行分析，旨在了解以水仙归化群体为代表的海岛植物地理隔绝状态下的群体遗传多样性水平，以及分布格局的可能形成机制，以期为海岛特色种质资源的科学评价、合理保护与经营利用提供分子遗传学依据。

2017年冬季，个体样本自北向南分别采自上海、浙江、福建的5个归化水仙代表群体，覆盖了水仙在中国东部沿海的主要分布区，共计49个样本，记录样品编号、采集地点、经纬度、花型等基本信息(表 1), 其中平潭岛(PTD)重瓣水仙为逸生的栽培品种‘玉玲珑’。群体内个体的采集距离至少相隔10 m以上，选择当年生新发的幼嫩叶片经变色硅胶干燥处理后，置于-20℃冰箱备用。活体材料引种至南京中山植物园，凭证标本存于南京林业大学树木标本馆(NF)。

表 1(Table 1)

表 1 5个水仙归化群体的采集信息 Table 1 Information of collected species for five naturalized Narcissus tazetta var. chinensis populations 群体 Population 采样地和代码 Sample location and code 个体数 Number 经度(E) Longitude 纬度(N) Latitude 花型 Flower type 采集人 Collectors 采集日期 Collection date POP 1 浙江大檑山屿 Daleishan Island, Zhejiang (DLSY) 10 121.09° 27.49° 单瓣 Single 段一凡, 朱弘 Duan Yi-fan, Zhu Hong 2017.01 POP 2 浙江南麂岛(三盘尾) Nanji Island, Zhejiang (NJD) 10 121.10° 27.46° 单瓣 Single 段一凡, 朱弘 Duan Yi-fan, Zhu Hong 2017.01 POP 3 上海崇明岛 Chongming Island, Shanghai (CMD) 10 121.77° 31.59° 单瓣 Single 王贤荣, 段一凡 Wang Xian-rong, Duan Yi-fan 2017.02 POP 4 福建平潭岛 Pingtan Island, Fujian (PTD) 9 119.77° 25.44° 重瓣 Double 王贤荣, 王华辰 Wang Xian-rong, Wang Hua-chen 2017.02 POP 5 福建漳州 Zhangzhou, Fujian (ZZ) 10 117.52° 24.54° 单瓣 Single 王贤荣, 王华辰 Wang Xian-rong, Wang Hua-chen 2017.02

表 1 5个水仙归化群体的采集信息 Table 1 Information of collected species for five naturalized Narcissus tazetta var. chinensis populations

植物基因组DNA采用柱式抽提纯化法试剂盒(Tiangen, 北京)提取, 用1%琼脂糖凝胶电泳和Tanon-2500全自动数码凝胶图像分析仪(Tanon, 上海)检测DNA质量。参考已报道的水仙属候选DNA条码^[13]进行筛选，预试验结果表明，水仙cpDNA编码区的matK (F: 5ʹ-CGATCTATTCATTCAATATTTC-3ʹ, R: 5ʹ-TCTAGCACACGAAAGTCGAAGT-3ʹ)与非编码区的trnH-psbA (F: 5ʹ-GTTATGCATGAACGTAATGCTC-3ʹ，R: 5ʹ-CGCGCATGGTGGATTCACAATCC-3ʹ)正反向引物片段在种内群体水平表现出较高的变异，且测序成功率最高，故选择上述序列进行研究。PCR扩增体系总体积为25 μL:包含2×Taq PCR MasterMix (TsingKe, 南京) 12.5 μL, 上下游引物各1.0 μL (GenScript, 南京)，DNA模板2.0 μL, 最后用ddH₂O 8.5 μL补齐。matK的PCR扩增程序为94℃预变性3 min，然后94℃变性30 s，41℃退火45 s，72℃延伸1 min，循环35次，最后72℃延伸10 min；trnH-psbA的PCR扩增程序为94℃预变性3 min，然后94℃变性30 s，47℃退火45 s，72℃延伸1 min，循环35次，最后72℃延伸10 min。产物经1%琼脂糖凝胶电泳检测合格后送至南京金斯瑞生物科技公司进行纯化与测序，样本基于ABI 3730XL测序平台进行正向测序。

将获得的单基因序列分别在BioEditor 7.0.9软件^[14]中的ClustalW模块进行联配与手工校正。使用SequenceMatrix 1.7.8软件包^[15]进行双基因序列的拼接。利用DNASP 5.10.1软件^[16]进行联合序列的特征分析，包括统计单倍型(haplotype)数量(h)、平均核苷酸差异数(average number of nucleotide, K)、单倍型多样性(haplotype diversity, H_d)、核苷酸多样性(nucleotide diversity, π)、遗传分化系数(G_st)、基因流(N_m)等遗传多样性指数以及种群历史动态的中性检验(neutrality test)和失配分析(mismatch distribution), 群体间基因流采用N_m=(1/G_st-1)/2^[20]进行估算。使用Arlequin 3.1软件^[17]进行分子方差分析以判别变异的来源。使用PermutCpSSR 2.0软件包(http://www.pierroton.inra.fr/genetics/labo/Software/Permut/)中的U统计法检验(1 000次重复置换)比较遗传分化系数G_st和N_st, 以检验水仙是否存在谱系地理结构。使用GenAlEx 6.503软件包^[18]在Microsoft Excel 2007中进行Mantal检验，分析水仙群体遗传距离与地理距离是否存在地理隔离(isolation by distance, IBD)作用。在Network 5.0.1.0软件(http://www.fluxus-engineering.com/sharenet.htm/)中基于最大简约性原则的中介邻接网络法(median-joining network, MJ)获得单倍型中介网络图，并在Adibe Illustrator CS6软件中进行优化布局与重新绘制。并利用MEGA 5.0软件^[19]构建以邻接法(neighbor-joining, NJ)为算法的单倍型系统发育树, 并采用bootstrap法(1 000次重复置换)以检验聚类树各分支置信度。

两个cpDNA片段拼接后的序列总长度为1 443 bp, G+C含量为33.10%，共检测出8个变异位点(7个简约信息位点和1个自裔位点)；共定义了6个单倍型(haplotypes, Hap1~Hap6)，其中matK对应的变异位点为673、767、801、807和809，trnH-psbA对应的变异位点为906、1 333和1 383。根据溯祖理论, 贡献单倍型(Hap2)位于单倍型中介网络中心, 应为5个地理群体的共有古老单倍型，而私自有单倍型则分散在各个群体，应为近期辐射分化出来的年轻单倍型(图 1, 表 2)。

图 1(Fig. 1)

图 1 中国东部沿海5个水仙归化居群的叶绿体单倍型分布 Fig. 1 Distribution of cpDNA haplotypes of five naturalized Narcissus tazetta var. chinensis populations in eastern coast of China

表 2(Table 2)

表 2 基于叶绿体matK+trnH-psbA序列的水仙6个单倍型中的变异位点 Table 2 Variable sites among the six haplotypes of Narcissus tazetta var. chinensis based on sequences of chloroplast matK+trnH-psbA 单倍型 Haplotype matK trnH-psbA 分布数量 Number 673 767 801 807 809 905 1333 1382 Hap1 A T G T G C T G 2 Hap2 G * * * * * * * 39 Hap3 G * * * * G * * 4 Hap4 G C A A T * * * 1 Hap5 G C A A T * A T 1 Hap6 G * G * * * A * 2 *:该碱基变异位点与Hap1相同。 *: Variation site of the base is the same as Hap1.

表 2 基于叶绿体matK+trnH-psbA序列的水仙6个单倍型中的变异位点 Table 2 Variable sites among the six haplotypes of Narcissus tazetta var. chinensis based on sequences of chloroplast matK+trnH-psbA

从表 3可见，群体的总体单倍型多样性(H_d)为0.363±0.086，核苷酸多样性(π)为0.000 49±0.000 19；按遗传多样性水平，5个水仙归化群体依次为大檑山屿(DLSY) > 平潭岛(PTD) > 南麂岛(NJD)=漳州(ZZ) > 崇明岛(CMD)。

表 3(Table 3)

表 3 5个水仙归化群体的遗传多样性 Table 3 Genetic diversityof five naturalized Narcissus tazetta var. chinensis populations 群体 Population 单倍型数(h) Number of haplotype 多态性位点数 Number of polymorphic sites (S) 平均核苷酸差异数(K) Average number of nucleotide differences 单倍型多样性 Haplotype diversity (Hd) 核苷酸多样性 Nucleotidediversity (π) DLSY 3 2 0.711 0.622±0.138 0.000 49±0.000 14 NJD 2 1 0.356 0.356±0.159 0.000 25±0.000 11 CMD 1 0 0.000 0.000±0.000 0.000 00±0.000 00 PTD 2 6 1.333 0.378±0.181 0.000 92±0.000 69 ZZ 2 1 0.356 0.356±0.159 0.000 25±0.000 11 总和Total 6 8 0.711 0.363±0.086 0.000 49±0.000 19

表 3 5个水仙归化群体的遗传多样性 Table 3 Genetic diversityof five naturalized Narcissus tazetta var. chinensis populations

从表 4可见，两两群体间的基因流较强, 且所有群体间的基因流都 > 4；分子方差分析(AMOVA)结果表明，群体内变异占了总变异的91.98%，为遗传变异的主要来源；固定指数(F_st)为0.080 22，群体间的遗传分化较弱(表 5)；通过1 000次重复计算获得的群体遗传分化系数(N_st=0.020 < G_st=0.031; P < 0.05), 表明水仙在物种水平上不具有明显的谱系结构。

表 4(Table 4)

表 4 5个水仙归化群体间的遗传分化系数(Gst, 对角线下)及基因流(Nm, 对角线上) Table 4 Genetic differentiation coefficient (Gst, under diagonal) and Gene flow (Nm, above diagonal) among five naturalized Narcissus tazetta var. chinensis populations 群体代码 Population code DLSY NJD CMD PTD ZZ DLSY - 12.000 4.222 4.385 4.000 NJD 0.040 - 4.222 4.454 4.000 CMD 0.106 0.106 - 4.808 4.222 PTD 0.102 0.101 0.094 - 5.445 ZZ 0.111 0.111 0.106 0.084 -

表 4 5个水仙归化群体间的遗传分化系数(Gst, 对角线下)及基因流(Nm, 对角线上) Table 4 Genetic differentiation coefficient (Gst, under diagonal) and Gene flow (Nm, above diagonal) among five naturalized Narcissus tazetta var. chinensis populations

水仙归化群体间的遗传距离(G_st)与地理距离(K_m, 自然对数转化)矩阵的Mantel检验结果表明, 两者之间具有显著的线性关系(r=0.929, P=0.02 < 0.05) (图 2)。

图 2(Fig. 2)

图 2 水仙归化群体遗传距离与地理距离的Mantel检测 Fig. 2 Mantel test between genetic distance and geographical distance of naturalized Narcissus tazetta var. chinensis populations

选择中性检验和失配分析推测水仙群体的扩张事件。通常用Tajima的D指数^[21]可以反映较大时间尺度上的种群事件，而Fu的D^*指数^[22]对近期的事件比较敏感。本研究结果表明，中性检验均无显著性意义(D=-0.165 395, P > 0.05; D^*=0.590 83, P > 0.10)。同时，错配分布分析(MDA)中，观测值与期望值比较得到的双峰曲线结果也表明水仙群体背离了快速扩张模型的假设，即群体在进化过程中未经历过快速扩张(图 3)。

图 3(Fig. 3)

图 3 水仙49个个体错配分布分析 Fig. 3 Mismatch distribution analysis for 49 individuals of Narcissus tazetta var. chinensis

基于遗传距离邻接法构建的单倍型系统发育树(图 4)可知，与中介网络图结果一致(图 1)，水仙6个单倍型明显被划分为2个谱系分支，Hap4与Hap5聚为一大支，其余单倍型聚为另一大支，靴带支持率(bootstrap)分别为80%和88%。其中Hap1、Hap3与Hap2又单独聚为一小支, 从而表现出更近的亲缘关系。

图 4(Fig. 4)

图 4 水仙6个单倍型的系统发育树。分支上数值为1 000次重复得到的靴带支持率。 Fig. 4 Phylogenetic tree of six haplotypes of Narcissus tazetta var. chinensis. Data above branches indicate the bootstrap values based on 1 000 replicates.

开展遗传变异分析有助于评价物种的进化潜力与环境适应能力，也是开展保护生物学研究的核心^[23-25]。基于cpDNA序列变异的结果，从整体上看, 我国归化水仙的遗传多样性水平较低，这与前人应用RAPD等其他标记手段的研究结果一致^[10]; Mantel检测的相关系数较高且达到显著水平(r= 0.929, P < 0.05)，符合地理隔离IBD模型，归化群体遗传距离更受地理距离影响显著，表明海岛间的地理隔离已经成为该种基因流扩散的主要屏障。细胞染色体研究表明，中国水仙主要为同源三倍体^[26], 在自然状况下水仙只能采取以母鳞茎为中心的侧球无性繁殖方式，大多呈聚集分布，缺乏长距离的传播机制^[27]，使得居群内交流远大于居群间。本研究中的分子方差分析结果亦显示水仙群体变异主要来源于种群内(91.98%)。此外，水仙栖息地严酷生境尤其是土壤的贫瘠及基盐胁迫无疑降低了水仙种群的生态适应性，而其物候花期寒冷的气温与强烈的海风等因素或可导致长距离的传粉者缺失, 奠基者效应的出现大大降低了群体遗传多样性。

植物在长期的迁移演化过程中为适应不同生境形成了不同的地理生态类群，生态隔离(isolation by ecology, IBE)作用往往可以加速物种的遗传分化^[28]。前人的研究均表明归化水仙在我国存在一定的遗传变异特征^[29]，陈晓慧等^[30]基于居群统计，分析比较了我国沿海4个岛屿8个群体的形态学特征差异, 认为南麂产的水仙具有植株小、花朵大的特点，是一种珍稀遗传资源。而结合本文的cpDNA单倍型分布来看(图 1)，各群体均含有1个以上的私有单倍型，暗示我国沿海岛屿的环境异质性较高，水仙居群可能存在局域适应(local adaptation)现象；水仙属植物通常性喜温暖、湿润及沙质土壤，在我国东部归化群体中，大檑山屿(DLSY)含有2个私有单倍型导致其遗传多样性综合排序第一，这可能是由于其地处南麂列岛国家级自然保护区核心区的地理优势，在区系划分上属于华东植物区系南北过渡特征，亚热带海洋性季风区及台湾暖流带来的水热气候条件相对充沛^[31]；基于群落生态学研究亦表明, 大檑山屿具有适宜水仙生长的排水较好、土壤盐分适中与温度适宜的微生境^[27]，可能使得大檑山屿群体产生了适应性分化。有意思的是，该群体与紧邻的南麂岛(NJD)共享一个单倍型(Hap3)，且基因流(N_m=12.000)最高、遗传分化(G_st=0.040)最小，相较其它群体而言显示出最近的亲缘关系；而虽然在物种水平上遗传多样性较低，但水仙各群体间的平均基因流均 > 4，说明群体间的基因流持续发生(表 5)，结合单倍型系统发育分析，聚类结果也与地理距离不一致(图 4)。根据实际调查走访得知，造成上述现象很可能与历史上海岛渔民人为多次移栽事件密切相关，即群体间可能经常性交换迁徙个体; cpDNA单倍型中介网络图还显示，平潭岛(PTD)与其余群体的亲缘关系相对最远，该群体是归化水仙中为数不多的重瓣类型‘玉玲珑’，其副花冠和雄蕊花瓣状，黄白相间，与大部分单瓣、纯白花色类型截然不同^[32], 有关其花色及花形变异的生理与分子机制还有待进一步研究，但无疑反映了人类经济活动, 对水仙遗传分布格局带来的持续影响。

表 5(Table 5)

表 5 5个水仙归化群体的分子方差(AMOVA)分析 Table 5 Analysis of molecular variance (AMOVA) for five naturalized Narcissus tazetta var. chinensis populations 变异来源 Source of variation 自由度 Degree of freedom 总方差 Sum of squares 变异成分 Variance component % 固定系数 Fixation index (Fst) 群体间Among populations 4 2.461 0.028 94 8.02 0.080 22 群体内Within populations 44 14.600 0.331 82 91.98 - 总和Total 48 17.061 - - -

表 5 5个水仙归化群体的分子方差(AMOVA)分析 Table 5 Analysis of molecular variance (AMOVA) for five naturalized Narcissus tazetta var. chinensis populations

目前在我国水仙市场上，栽培品种中除了‘玉玲珑’，还有‘金三角’、‘金边’和‘云香’等，但都分布不均，栽培面积极少且产地不明，品种类群十分单一；而且由于水仙很难通过有性繁殖同其他物种进行杂交获得新的基因，长期的无性繁殖，势必导致种群衰退、品质降低等问题，这将严重制约水仙花卉产业的可持续发展^[33]。因此，开展我国归化水仙资源种群动态和遗传变异的研究进而获得基础的遗传学背景信息是非常必要的。本研究综合利用了植物分子谱系地理学的方法，首次对我国东部沿海水仙5个代表归化群体49个样本进行遗传多样性检测及分析，初步阐明了该物种在区域尺度上的遗传多样性水平与空间分布，这无疑有助于下一步水仙衰退种群的恢复保育工作的开展。同样以南麂岛产水仙尤其是大檑山屿为代表，基于群体遗传cpDNA标记的分子证据，该归化群体应为中国东部沿海一个特殊且重要的基因型或生态型分化类群，目前虽然已经列入核心保护区，但由于靠近旅游区仍然面临人类经济活动的日益影响，应结合沿海岛屿特殊的生态系统进行综合管理，做好持续的宣传、保护与管理工作。短期而言，就地保护可能是最为有效的措施，可以在保护区内增加护栏、控制伴生入侵杂草及定期巡护工作；长期而言，一方面通过原地建立固定样地，开展对种群及物候的动态监测，另一方面通过离体培养、人工种苗扩繁、自然回归试验等措施，从而促进其种群的正向更新。同时也应发挥其重要的科学价值，将该区域优先选为核心种质苗圃或遗传基因库建设的遗传材料来源地，或可从中选择适宜的亲本以提高杂交育种效率。总之, 本研究可为海岛特色种质资源的综合评价、科学保护与开发利用提供理论依据。