喜马拉雅—横断山区优越虎耳草谱系地理学研究

贾留坤 焦秀洁 刘海瑞 李 彦 更吉卓玛 王智华 高庆波

(1.中国科学院西北高原生物研究所中国科学院高山植物适应与进化重点实验室，西宁 810001； 2.中国科学院大学，北京 100039； 3.中国科学院西北高原生物研究所信息与学报编辑部，西宁 810001； 4.中国科学院西北高原生物研究所青海省作物分子育种重点实验室，西宁 810001)

自Avise于1987年提出谱系地理学概念(phylogeography)[1]，经过40多年发展，该学科已成为生物地理学的一个重要分支。谱系地理学借助分子生物学技术，广泛应用直接测序来检测核DNA 转录间隔区(ITS)和叶绿体DNA非编码区的序列变异[2]，将微观进化与宏观进化相结合，研究物种现有遗传分布格局的历史成因，从进化历程、时间和空间尺度的动态变化等方面重建生物区系进化历史[3～4]。

青藏高原是世界上平均海拔最高的地区，该区域地理结构复杂、生境多样、生态位分化明显，具有较高的生物多样性，是世界生物多样性热点地区[5]，且拥有大量的特有种，是世界上特有种密度最高的地区之一[6]。青藏高原的快速隆升被认为对高原物种分化产生了重要影响，使得该区域成为研究高山植物适应性进化的理想地区，相关研究也发现在该地区的特色类群如虎耳草属(SaxifragaL.)[7]、马先蒿属(PedicularisL.)[8]、风毛菊属(SaussureaDC.)[9]、红景天属(RhodiolaL.)[10]中均能检测到明显的快速辐射分化。由于复杂的地理拓扑结构和独特的地形，青藏高原在末次盛冰期受到的影响与欧洲和北美大陆不同，未形成统一的大冰盖[11]，使得该地区植物对于第四纪冰期的响应呈现出不同的谱系地理模式。“退缩—回迁”模式认为，高原台面上的居群在冰期时退缩到高原边缘，冰期后或间冰期再由避难所迁移至高原台面，如祁连圆柏(Juniperusprzewalskii)[12]、青海云杉(Piceacrassifolia)[13]、长花马先蒿(Pedicularislongiflora)[14]、窄叶鲜卑花(Sibiraeaangustata)[15]等类群的谱系地理研究结果符合该模式。通常，符合“退缩—回迁”模式的物种，一般在高原边缘的避难所内拥有较多的单倍型和较高的遗传多样性，由于回迁过程中受到奠基者效应影响，高原台面上仅固定少数单一的单倍型，并呈现出较低的遗传多样性。“微型避难所”模式认为，某些耐寒的植物在冰期时并未退缩到高原边缘，而是在高原台面若干个相互隔离的微型避难所内保留下来，在冰期后或间冰期时，从避难所向周边扩散，并由于地理阻隔而发生不同程度的异域分化，该模式的代表物种包括菊叶红景天(Rhodiolachrysanthemifolia)[16]、唐古特虎耳草(Saxifragatangutica)[17]、西藏沙棘(Hippophaetibetana)[18]等。拥有微型避难所的物种一般拥有较多的特有单倍型和多个随机分布的遗传多样性较高的居群。因此，研究青藏高原地区高山植物的进化历史不但可以为该地区植物遗传资源的保护和利用提供理论基础，也可为高原大规模隆升的时间尺度和第四纪冰期是否存在‘大冰盖’等科学问题提供间接证据。

虎耳草属(SaxifragaLinnaeus)约有500种，主要分布于北半球高山地区。该属起源于北美的落基山，向南扩散到南美洲的安第斯山和火地岛，在北部经白令海峡扩散到欧亚大陆[19]，分布范围包括欧洲南部的阿尔卑斯山、高加索山，青藏高原地区，北极圈附近也有分布。我国有虎耳草属植物216种，主要分布在青藏高原—喜马拉雅地区[20]。为了适应高原地区复杂的生态环境，该属植物在生活史、形态、生理特征等方面表现出极高的多样性。已有研究表明，虎耳草属植物的分化与青藏高原的隆升以及该地区复杂的地理结构有关，因此该属植物被广泛应用于谱系地理学和系统发育学研究，来揭示北极和高山地区植物的分化式样[21]。优越虎耳草(SaxifragaegregiaEngl.)隶属于虎耳草科(Saxifragaceae)虎耳草属(SaxifragaLinnaeus)山羊臭组(sect.CiliataeHaw.)，为多年生草本，植株高9～32 cm。基生叶具长柄，叶片多为心形、心状卵形，边缘具卷曲长腺毛，最上部叶片披针形至长圆形。多歧聚伞花序伞房状，花梗被短腺毛。萼片在花期反曲，卵形至阔卵形，背面和边缘具腺毛。花瓣黄色，花丝钻形。子房近上位，卵球形，花期7～9月。主要分布于青藏高原边缘地区密闭的林下、灌丛，海拔2 800～4 500 m。主要范围包括甘肃南部、青海东部和南部、四川西部、云南西北部和西藏东部。模式标本采自青海、甘肃交界的大通河流域[22]。

1 材料与方法

1.1 实验材料

研究材料采集于2003～2012年，共计13个居群、151个个体。采集地点主要集中在青藏高原东南部的横断山地(表1，图1)，包括青海省南部一个居群、西藏东部两个居群和四川西部10个居群。每个居群随机采集3～24个个体，个体之间间隔10米以上，选取长势较好植株的新鲜幼嫩叶片，随即用硅胶干燥，带回实验室后放于-20℃冰箱中保存。凭证标本馆藏于中国科学院西北高原生物研究所青藏高原生物标本馆(HNWP)，并由Richard J. Gornall教授和高庆波副研究员共同鉴定。

表1优越虎耳草13个居群的采集信息

Table1Samplinginformationof13populationsofS.egregia

居群Pop.凭证标本Voucher采样地Location地理坐标Coordinates海拔Altitude(m)采样数No.P1chen03026青海班玛Banma,QinghaiN32°48'34″;E100°56'50″368019P2chen03075四川色达Seda,SichuanN31°52'02″;E100°41'23″336024P3chen03084四川康定Kangding,SichuanN30°15'17″;E101°30'37″355010P4chen06111四川红原Hongyuan,SichuanN31°52'46″;E102°40'22″397015P5chen06242四川乡城Xiangcheng,SichuanN30°41'11″;E97°15'42″338019P6chen06290四川雅江Yajiang,SichuanN30°09'38″;E100°40'45″429022P7chen2007049西藏类乌齐Leiwuqi,XizangN31°14'23″;E96°33'56″38208P8chen2007202西藏昌都Changdu,XizangN31°24'27″;E97°28'13″36709P9chen2007219四川德格Dege,SichuanN31°57'21″;E98°51'48″41108P10chen2012252四川理塘Litang,SichuanN30°15'31″;E100°15'47″38903P11chen2012268四川道孚Daofu,SichuanN30°43'02″;E101°21'43″37706P12chen2012288四川炉霍Luhuo,SichuanN31°37'29″;E100°43'12″34704P13chen2012293四川壤塘Rangtang,SichuanN32°07'51″;E100°00'08″31404

图1 优越虎耳草cpDNA联合片段19个单倍型地理分布图 饼状图表示每个居群中单倍型的频率Fig.1 Geographic distribution pattern of the 19 detected cpDNA haplotypes in S.egregia Pie charts show the proportions of haplotypes within each population.

1.2 DNA提取与PCR扩增

采用改良的 CTAB法[23]，从硅胶干燥叶片中提取总DNA。用NanoDrop 2000(Thermo Scientific)逐一检测每个DNA样品的质量与浓度，并将所有样品稀释到50 ng·μL-1，用于随后的PCR扩增。采用通用引物分别对所有个体的cpDNArpl20-rps12和trnL-trnF两个片段进行PCR扩增[24～25]。PCR扩增反应体系为50 μL：包含5 μL的10 × PCR Buffer(含1.5 mmol·L-1Mg2+)，2 μL的10 mmol·L-1dNTPs，5 pmol·L-1的正反引物各2 μL，0.4 μL的Taq DNA聚合酶(5 U·μL-1)，0.4 μL(约20 ng)的DNA模板，用双蒸水补足到50 μL。PCR扩增反应程序为：94℃预变性5 min；94℃变性1 min，55℃退火50 s，72℃延伸1 min，32循环；72℃延伸10 min。取2 μL扩增产物进行1.5%琼脂糖凝胶电泳，检测扩增产物的产量和质量。剩余的48 μL PCR扩增产物直接送北京三博远志生物技术有限责任公司进行割胶纯化并进行双向测序。

1.3 数据分析

测序所得序列首先用Chromas软件和手工进行核对，再将cpDNA两个片段(rpl20-rps12和trnL-trnF)进行串联合并，用MEGA6.06软件进行对位排列，并加以手工校正[26]。用DnaSP 5.1软件统计单倍型[27]。居群内平均遗传多样性(average gene diversity within populations，Hs)，总的遗传多样性(total gene diversity，Ht)以及居群遗传分化系数Gst和Nst均在软件PERMUT中计算[28]。

用Arlequin 3.5.1.3软件计算居群的遗传多样性(gene diversity，He)及核苷酸多样性(nucleotide diversity，π)[29]，并用AMOVA(analysis of molecular variance)软件检测遗传变异在居群间和居群内的分布，并对Fst进行评价，以进一步揭示居群的分化程度(1 000次重复的显著性检验)[30]。用软件DnaSP 5.1进行歧点分布分析，并在Arlequin 3.5.1.3软件包中计算中性检验Tajima’s D和Fu’sFs值[31～32]。

以最大简约法为原则，通过NETWORK 4.6.1.3构建单倍型之间的系统发育关系图[33]。构建网络图时包含了所有信息位点，插入/缺失和点突变设置相同的权重。在单倍型中央链接网络图中，年轻的单倍型发生在网络图的边缘，古老的单倍型多存在于网络图的内部[34]。

2 结果与分析

2.1 单倍型变异位点和地理分布格局

对优越虎耳草13个居群、151个个体cpDNA片段rpl20-rps12和trnL-trnF的测序结果串联合并，对位排列后序列长度为1 599 bp。通过DnaSP软件分析，共得到19个单倍型(图1)，检测出36个变异位点(表2)。19种单倍型的频率变化范围为0.007～0.318；其中频率最高的单倍型H2，在48个个体中出现；单倍型H6次之，在21个个体中出现；单倍型H11、H14、H15、H19分别只在一个个体中出现；其余单倍型则在2～15个个体中出现。在19个单倍型中，12个为居群特有单倍型。这些特有单倍型并未固定于某个区域，而是在优越虎耳草自然分布范围内随机分布。每个居群的遗传多样性(He)、核苷酸多样性(π)、单倍型组成及频率见表3。

2.2 单倍型间的系统关系

最大简约法(Maximum parsimony,MP)、最大似然法(Maximum likelihood,ML)和贝叶斯法(Bayesian inference,BI)构建的系统发育树均没有很好地揭示出优越虎耳草各单倍型之间的系统发育关系，可能是由于优越虎耳草各单倍型之间的分歧时间较短，没有积累足够的变异所致。NETWORK软件能更好的揭示分化较浅的序列之间的系统发育关系，于是我们利用最大似然法对优越虎耳草19个单倍型构建中央链接网络图(图2)。结果显示，出现频率最高的单倍型H2位于网络图的中心位置，推测H2为古老的单倍型。居群特有单倍型一般位于网络图的边缘，可能是通过近期辐射分化衍生而来的年轻单倍型。

2.3 居群的遗传结构

利用PERMUT软件对优越虎耳草19个单倍型进行分析得到居群内平均遗传多样性(Hs)为0.466，总的遗传多样性(Ht)为0.868，居群间遗传分化Gst值和Nst值分别为0.438和0.463。使用U-统计方法对优越虎耳草整个分布区单倍型变异的地理分布模式进行检验后发现Nst大于Gst，但不显著(P>0.05)，表明优越虎耳草在其现有分布区域内没有明显的谱系地理学结构。分子变异分析(AMOVA)显示，居群内遗传变异为57.37%，居群间遗传变异为42.63%(表4)，表明优越虎耳草cpDNA联合片段的遗传变异主要存在于居群内(Fst=0.426；P<0.001，1,000次重复的显著性检验)。

表2 优越虎耳草2个cpDNA片段19种单倍型进行序列比对后的变异位点

注：序列顺序为5′-3′ 表中的数字表示：1. TTATTATTATGTATACGTA; 2. AAATACTATATAAAAAT; 3. AATTATTA; 4. ACTAAAAAAAAA; 5. AAAAAAG; 6. AA; 7. AAAA; 8. TTTTT; 9. TTTT; 10. ATTAATC; 11. TTTGACTCCTTTTGACTCCT; 12. AACTAT; 13. CCTATACTTCAATT; 14. CCTATACTTCAAATT; 15. TATATA; 16. TTCCTGCAT; 17. CCCTA; 18. TCCTA; 19. ATTCATATTCTATTTCTATCCT; 20. ATTCATATTCTATTTCTATCCTA; 21. ATTCATATTCTATTTCTAGCCT; 22. ATTCATATTCGATTTCTATCCT

Note:Sequences are numbered from 5′-3′ end Representations of No. in the table are as follows: 1. TTATTATTATGTATACGTA; 2. AAATACTATATAAAAAT; 3. AATTATTA; 4. ACTAAAAAAAAA; 5. AAAAAAG; 6. AA; 7. AAAA; 8. TTTTT; 9. TTTT; 10. ATTAATC; 11. TTTGACTCCTTTTGACTCCT; 12. AACTAT; 13. CCTATACTTCAATT; 14. CCTATACTTCAAATT; 15. TATATA; 16. TTCCTGCAT; 17. CCCTA; 18. TCCTA; 19. ATTCATATTCTATTTCTATCCT; 20. ATTCATATTCTATTTCTATCCTA; 21. ATTCATATTCTATTTCTAGCCT; 22. ATTCATATTCGATTTCTATCCT

表3 优越虎耳草13个居群的遗传多样性、核苷酸多样性和单倍型组成

图2 优越虎耳草19个单倍型的中央链接网络图Fig.2 Median-joining network for the 19 detected haplotypes in S.egregia

表4优越虎耳草13个居群的分子遗传变异(AMOVA)分析结果

Table4Analysisofmolecularvariation(AMOVA)forthe13populationsofS.egregia

变异来源Source of variation自由度Degrees总方差Sum of squares变异成分Variation component变异比例Variation(%)居群间Among populations1275.2100.4970942.63居群内Within populations13892.3020.6688557.37总计Total150167.5121.16595固定系数Fixation index(Fst)Fst=0.426(P< 0.001)

2.4 居群的历史动态

在优越虎耳草13个居群中，有11个存在居群内遗传分歧(intra-population divergence)，对所有测序成功的个体和具有居群内遗传分歧的居群进行中性检验，在整个分布范围内，Tajima’s D和Fu’s Fs值均不显著，表明优越虎耳草在整体水平上没有经历过近期扩张。在居群水平上，只有居群P3的Tajima’s D为显著负值，其余居群的两个中性检验值均不显著，表明优越虎耳草在居群水平上也没有经历过近期扩张(表5)。利用Arlequin软件进行歧点分布分析，结果表明优越虎耳草在整个分布范围内SSD值(0.041 8)和Raggedness index值(0.039 0)均显著，表明该物种在整个分布范围内未经历过近期扩张，这与中性检验结果一致。在居群水平上，居群P1、P2、P3、P5、P7、P9和P11实际观测到的歧点分布与期望的居群扩张模型可以显著分开，表明这7个居群没有经历过大规模的居群近期扩张；居群P4、P6、P8、P10则不能排除经历过近期居群扩张(表5，图3)。

表5优越虎耳草151个个体及具有居群内遗传变异的11个居群cpDNA联合片段中性检验、歧点分布分析结果

Table5Resultsoftheneutralitytestsandmismatchdistributionanalysisforthe151individualsandthe11populationswhichhaveintra-populationdifferentiation,usingthecpDNAcombinedfragmentsofS.egregia

Population codeTajima’s DFu’s Fs*SSDRaggedness indexTotal individuals-1.3483218.915720.0418*0.0390*P11.8196829.364680.24610.4270*P22.0853721.079840.5444*0.7732P3-1.66706*8.858420.0577*0.7200P4-0.02381-0.064330.03330.2259P5-1.0109011.459430.5702*0.1128P6-0.641126.548050.01030.5634P70.4142110.019810.2483*0.6939P80.156471.724670.08250.1975P9-1.310095.547260.5332*0.4834P100.000002.210060.30170.6667P110.8505710.327100.5689*0.7867

注：*表示P>0.05

Note:*signP>0.05

图3 优越虎耳草151个个体cpDNA联合片段歧点分布分析结果Fig.3 Mismatch distribution analysis for the 151 individuals based on the cpDNA combined fragments of S.egregia

3 讨论

3.1 优越虎耳草的居群遗传结构

通过遗传变异分析发现，优越虎耳草总的遗传多样性较高(Ht=0.868)，而居群内平均遗传多样性较低(Hs=0.466)。居群遗传结构符合青藏高原地区植物普遍的遗传变异类型，如长花马先蒿(Pedicularislongiflora)[14]、露蕊乌头(Aconitumgymnandrum)[35]、祁连圆柏(Juniperusprzewalskii)[12]、西藏沙棘(Hippophaetibetana)[18]、西川红景天(Rhodiolaalsia)[36]、银露梅(Potentillaglabra)[37]、Juniperustibeticacomplex[38]和沙棘(Hippophaerhamnoides)[39]等，均具有较高的总遗传多样性和较低的居群内平均遗传多样性。cpDNA属于母系遗传，基因流主要通过种子传播。优越虎耳草为低矮草本，多生于郁闭林下和灌丛，且种子没有有利于传播的附属结构，传播距离较近，使得基因流主要发生在居群内，居群间基因流较低，最终使种内遗传分化增大，总的遗传多样性也较高，而居群内遗传多样性较低。在优越虎耳草中只有单倍型H2和H4在居群中广泛分布，并有12个为居群特有单倍型，占比达63.2%。通常，特有单倍型比例较高的物种往往不能揭示出明显的谱系地理格局[36～37,39]。本研究对优越虎耳草的两个遗传分化系数Gst和Nst进行1 000次的置换检验，发现Nst>Gst(P>0.05)，但不显著，表明优越虎耳草在其分布范围内不存在明显的谱系地理格局。此外，中央链接网络图显示，居群特有单倍型一般位于网络图的边缘，可能是通过近期辐射分化衍生而来的年轻单倍型，较高比例的年轻特有单倍型可能掩盖了优越虎耳草在整个分布范围内的谱系地理结构。青藏高原地区特有单倍型比例较高的其他物种，如西川红景天(Rhodiolaalsia)[36]，唐古特虎耳草(Saxifragatangutica)[17]等也未能揭示明显的谱系地理结构。

3.2 冰期避难所及居群进化历史的推测

通过对已有研究进行总结，发现与喜马拉雅—横断山区复杂的地质历史以及气候变化相对应，该地区植物具有复杂多样的进化历史[40]，其中可能存在的主要谱系地理模式有：冰期高原台面的种群退缩到高原边缘的避难所模式，其遗传结构表现为遗传多样性由高原边缘向高原台面递减，符合该进化模式的物种其居群历史一般是在冰期气候寒冷的时候由高原台面向高原东南部、东部的低海拔地区迁移，并在此得以保存，在间冰期或冰期后气候回暖时再向高原台面扩散，回迁过程中受到奠基者效应的影响，最终形成今天的谱系地理结构；冰期高原台面的种群退缩到本地避难所(含微型避难所)的模式，其遗传结构往往表现为特有单倍型和具有较高遗传变异的居群在物种分布范围内均匀分布，符合该进化模式的物种多为某些耐寒植物，第四纪冰期时并未迁移到高原边缘，而是在高原台面若干个彼此隔离的避难所保留下来，在间冰期或冰期后气候回暖时从避难所向周边扩散，并且发生不同程度的异域分化。

本研究发现优越虎耳草符合“微型避难所”假说。大量的特有单倍型在整个分布范围内随机分布，遗传多样性和核苷酸多样性较高的居群也呈均匀分布，没有出现遗传多样性较高的某一区域，也没有呈现遗传多样性由高原边缘向高原台面递减的趋势。根据歧点分布分析，在整体水平上SSD值和Raggedness index值均显著，在居群水平上，P1、P2、P3、P5、P7、P9和P11也表现出显著的SSD值和(或)Raggedness index值，可以将实际观测到的歧点分布与期望的居群扩张模型显著分开，表明居群没有经历过大规模的近期扩张现象。在优越虎耳草整体水平和居群水平进行的中性检验结果显示，Tajima’ D和Fu’s Fs值均不显著，与歧点分布分析结果相似。综合以上证据，我们推测优越虎耳草在第四纪冰期时存在多个微型避难所。值得注意的是P2、P12、P13的单倍型多样性明显低于其他居群，只有P2具有显著的SSD值，没有经历近期的扩张，即存在剩余两个居群是由附近居群扩散而来的可能性，在扩散的过程中由于奠基者效应，使得只有某一单倍型被保留了下来。或是这3个居群冰期时在该地保留下来，冰期结束后的扩散过程中经历了严重的瓶颈效应，在此过程中只有某一单倍型被保留，形成了今天的单倍型组成。

横断山区高山和大河交错分布，地理拓扑结构复杂，被认为是第三纪重要的物种分化中心，也被认为是很多植物在冰期的重要避难所[13,41]。而优越虎耳草主要分布在喜马拉雅—横断山区，生于温暖潮湿的林下。第四纪冰期时横断山区虽然没有被大规模的冰川覆盖，但寒冷干旱的冰期与温暖湿润的间冰期交替发生，气候的剧烈波动导致横断山区木本植物分布范围的多次缩小和扩张。我们推测生于林下的优越虎耳草的分布范围也必然随之发生变化，在此过程中居群不断片段化，最终在相互隔离的“高山岛屿”中发生异域分化，使该物种产生了大量特有单倍型。

优越虎耳草的主要生境为郁闭的林下及灌丛，然而到目前为止，青藏高原地区草本植物的谱系地理学研究中，对生境为灌丛或森林的涉及很少。而据我们的推测，优越虎耳草应该与同域分布的乔木或灌木具有相同或相近的谱系地理历史。Opgenoort利用cpDNA作为分子标记对与优越虎耳草分布区重叠的Juniperustibeticacomplex(刺柏属)进行谱系地理学研究，发现其特有单倍型比例(64.5%)、特有单倍型和具有较高遗传多样性的居群的分布与我们的研究结果相似[38]。Opgenoorth据此推测末次盛冰期时J.tibeticacomplex在青藏高原存在多个微型避难所，有些微型避难所的海拔超过了3 500 m，这也是已知最高的末次盛冰期树线[38]。此外，西藏沙棘的谱系地理学研究揭示了相似的结果，该物种在第四纪冰期时同样存在多个微型避难所，某些避难所海拔甚至超过了4 000 m[18]。由此我们可以推测，在第四纪冰期时如果某些乔木和灌木能在高原台面的多个微型避难所生存下来，那么生于林下的优越虎耳草也可能因为这些乔木或灌木提供的微型避难所而在冰期时保留下来。此外，另一个生长于林下和灌丛的物种—菊叶红景天(Rhodiolachrysanthemifolia)具有几乎完全相同的谱系地理结构。然而生于森林、灌丛下的草本植物与乔木和灌木之间是如何共同演化的，尚需进一步研究。