國立臺灣大學生物資源暨農學院森林環境暨資源學系 碩士論文
School of Forestry and Resource Conservation College of Bioresources and Agriculture
National Taiwan University Master Thesis
氣候變遷對台灣高山植物分布及遺傳多樣性影響之評估
──以山薰香屬為例
Climate Change Impact Assessment on Species Distribution and Genetic Diversity of Alpine Plants of Taiwan
──Insight from Oreomyrrhis Clade of Chaerophyllum
許正德 Cheng-Te Hsu
指導教授﹕鍾國芳 博士 Advisor: Kuo-Fang Chung, Ph.D.
中華民國 102 年 7 月
July 2013
謝誌
乍入科學殿堂，諸多碰撞試誤，完全不懂的概念和方法便是藉由老師的教 導、同學的砥礪、以及自我的醒悟，大抵這三年就是這麼交代的，眾多感謝和感 動以下一一敘說。
首先感謝恩師鍾國芳老師在各方面的幫助和指導，包括提供豐富的材料及經 費，讓我的研究得以進行，並且亦師亦友的協助我做學問的道理，畢生難忘。接 著感謝口試委員江友中老師及丁宗蘇老師在論文修改以及分析方法的指正。再 者，感謝一同努力的研究室同窗們，你們的支持與包容是我前進的動力！感謝婷 婷在研究之初的提點。謝謝政道在撰寫程式上的指點及協助，開啟我對於寫程式 的興趣。旨价包容我的大大小小的懶惰，並且鼓勵我邁出每次的猶豫，以及每次 難能可貴的激辯。謝謝章璿、筱蕾在研究上的討論以及不斷供應糧食。感謝瀚嶢 共同承擔六個學期樹木學加上植物分類學實習的助教工作，並且互相學習與進 步。謝謝建華多次的擔任採集行程的安全駕駛。以及感謝慧舟（小舟）、永信、
姿麟、雅青、佩儀、馨慧和柏璋在研究上和精神上的鼓勵。另外要致謝於大學的 同窗們，感謝彥睿（小班）在研究內容上的眾多討論，以及數次的晚餐、宵夜對 談。與俊佑、劉鎮共同奮戰，讓我在夜晚的系館不孤獨。要感謝的人太多，若有 疏漏還請見諒。
最後，要感謝我的家人背後的支持，即便諸多不諒解，但多半屬擔心，謝謝 人生之路上的養育與栽培。
然而這不是個結束，研究之路迢迢，且戰且走，自然社會中還有許多的奧秘 待我們去發現，希望自己仍保有簡單的好奇心，步步省思大自然傳達的訊息。因 此，感謝山神！
I
摘要
全球氣候的持續暖化將迫使高山植物往更高海拔遷徙，導致其分布面積大幅 縮減，因此評估指出高山植群帶是全球氣候變遷趨勢下最脆弱的生態系之一，然 而物種之存續尚繫於種內之遺傳多樣性及演化潛力以因應環境的改變，但鮮少研 究探討高山植物分布面積縮小對物種遺傳多樣性的影響。本研究以繖形科山薰香 屬植物為例，根據物種採樣點座標及 WORLDCLIM 氣候資料，以 MAXENT 程 式建立物種分布模型，並將之投影至上次最大冰期（Last Glacial Maximum, LGM）
的氣候模型及兩組預測的未來氣候情境，計算預測山薰香及台灣山薰香的分布，
評估其脆弱度，並整合物種分布預測與親緣地理資料，評估在不同暖化情境下遺 傳多樣性的喪失，其中南湖山薰香由於稀有而資料甚少，並未納入分析。預測結 果顯示，山薰香及台灣山薰香在 LGM 時族群較現今大幅擴張，而在未來不同的 暖化情境下，兩物種分別向上遷徙了 18–497 公尺及 17–476 公尺，分布面積則分 別縮小了 44.7–85.8％及 49.2–99.9％，且以東北部山區的脆弱度最高而漸次向南 遞減。根據分布面積喪失的預測，山薰香及台灣山薰香在葉綠體 atpB-rbcL 基因 間片段之基因單型的多型性在未來分別將有 5.3–21.1％及 0–80％的喪失。由於東 北部同時為遺傳多樣性的熱點，未來保育應加強東北部族群的監測，並針對即將 喪失基因單型的族群進行採樣與保種，以維持物種之遺傳多樣性。
關鍵字
物種分布模型、MAXENT、全球暖化、山薰香屬、繖形科、脆弱度、保育遺 傳學
II
Abstract
Studies have predicted that alpine plants will be forced to migrate into higher ele- vations as current trend of global warming continues, resulting in drastic deduction in their distribution range. Consequently alpine vegetation is regarded as one of the most vulnerable ecosystems threaten by global climate changes. However, species’ survivor- ship and evolvability to cope with environmental changes also depend on genetic diver- sity. Therefore it is also important to project how shrinking in distribution range would impact extent of genetic diversity. Based on coordinates of collecting data and climate data of WORLDCLIM, we employed the program MAXENT to construct the distribu- tion models of two endemic species of Apiaceae in high elevation of Taiwan,
Chaerophyllum involucratum and C. taiwanianum, and used the models to project their
past and future distributions during the last glacial maximum (LGM) and under the fu- ture climatic scenarios (i.e. B2B and A2A) predicted by Intergovernmental Panel on Climate Change (IPCC). Changes of projected distribution were used to evaluate the vulnerability of the species. The distributions were also integrated with phylogeographic data to predicted loss of genetic diversity under B2B and A2A. Based on the projections, distributions of both species were much wider during LGM compare to their present distributions. Under B2B and A2A, 18–497 m upward migration and 44.7–85.8%shrinking in distribution are projected in C. involucratum, while 17–476 m and 42.2–
99.9% in C. taiwanianum are estimated. Populations in the northeastern part of the is- land are most vulnerable, with vulnerability decreasing southward. Based on deduction of distribution range, C. involucratum is projected to lose 5.3–21.1% of haplotypes in the chloroplast atpB-rbcL spacer, while 0–80% loss was predicted in C. taiwanianum.
Because the northeastern populations also present hotspots of genetic diversity, their
III
high vulnerability urges attention to monitor whether populations there would have shrunk as predicted. Furthermore, ex situ conservation such as germplasm preservation in seedbank should focus on populations harboring haplotypes that are projected to go extinct in the near future.
Keywords
species distribution model, global warming, conservation genetics, Oreomyrrhis, vulnerability, Apiaceae, MAXENT
i
目錄
摘要 ... I Abstract ... II
目錄 ... i
圖目錄 ... iii
表目錄 ... iv
附錄目錄 ... v
前言 ... 1
I. 氣候變遷對高山生態系的衝擊 ... 1
II. 脆弱度評估 ... 2
III. 物種分布預測模型 ... 3
IV. 台灣高山植物之親緣地理 ... 4
V. 台灣山薰香屬簡介 ... 5
研究目的 ... 7
材料與方法 ... 8
I. 研究材料 ... 8
A. 生物資料取得 ... 8
B. 環境資料取得 ... 10
II. 資料分析 ... 11
A. 物種及環境資料前處理 ... 11
B. 族群遺傳變異分析 ... 12
C. 巢狀支序親緣地理分析 ... 12
D. 物種分布預測模型建立 ... 13
E. 多樣性評估分析 ... 14
ii
F. 脆弱度評估 ... 15
結果 ... 17
I. 族群遺傳變異分析、親緣網狀階層分析 ... 17
II. 物種分布預測模型時間序列 ... 18
A. 環境因子選取 ... 18
B. 預測分布模型 ... 19
III. 多樣性評估分析 ... 19
A. 分布面積動態 ... 19
B. 遺傳多樣性喪失評估 ... 20
IV. 脆弱度評估 ... 21
A. 脆弱度指數 ... 21
B. 海拔梯度機率頻度分布 ... 21
C. 棲地破碎程度分析 ... 22
討論 ... 23
I. 台灣山薰香屬遺傳多樣性成因 ... 23
II. 物種分布模型 ... 24
A. 資料誤差 ... 24
B. 模型結構 ... 24
III. 多樣性評估分析 ... 25
結論 ... 27
參考文獻 ... 28
iii
圖目錄
圖一、氣候變遷衝擊下的物種脆弱度評估架構（Williams et al., 2008） ... 35
圖二、本研究方法的流程圖 ... 36
圖三、本研究使用的生物氣候資料集。 ... 37
圖四、山薰香的 atpB-rbcL 基因單型地理分布 ... 38
圖五、山薰香的親緣網狀階層分析 ... 39
圖六、台灣山薰香的基因單型地理分布 ... 40
圖七、台灣山薰香的親緣網狀階層分析 ... 41
圖八、山薰香（左）及台灣山薰香（右）的現今平均機率預測分布。 ... 42
圖九、山薰香遺傳多樣性喪失比例，以採集座標計量 ... 43
圖十、山薰香遺傳多樣性喪失比例，以基因單型計量 ... 44
圖十一、台灣山薰香遺傳多樣性喪失比例，以採集座標計量 ... 45
圖十二、台灣山薰香多樣性喪失比例，以基因單型計量 ... 46
圖十三、山薰香脆弱度指數空間分布 ... 47
圖十四、台灣山薰香脆弱度指數 ... 48
圖十五、山薰香海拔梯度機率頻度分布 ... 49
圖十六、台灣山薰香海拔梯度機率頻度分布 ... 50
圖十七、山薰香各時期二元預測分布的內緣比（Interior-edge ratio） ... 51
圖十八、台灣山薰香各時期二元預測分布的內緣比（Interior-edge ratio） ... 52
圖十九、間冰期避難所模型 ... 53
圖二十、山薰香在 SRES 未來氣候 A2A 假設情境下，CGCM2 前、中及後期的親 緣網狀關係圖破裂情形 ... 54
iv
表目錄
表一、台灣山薰香屬葉綠體 DNA atpB-rbcL 內轉錄區間進行聚合酶連鎖反應使用
的引子 ... 55
表二、19 項生物氣候環境因子代號表及描述 ... 56
表三、山薰香及台灣山薰香基因單型之地理分布 ... 57
表四、山薰香各族群遺傳變異總表 ... 58
表五、山薰香族群間遺傳分化指數（FST） ... 59
表六、山薰香網狀階層分析 ... 60
表七、台灣山薰香屬植物之 19 項生物氣候環境因子 ... 62
表八、山薰香的 19 個生物氣候環境因子皮爾遜積差相關係數矩陣 ... 63
表九、山薰香的 19 個生物氣候環境因子皮爾遜積差相關係數矩陣 ... 64
表十、預測模型評估 ... 65
表十一、山薰香多樣性喪失分析 ... 66
表十二、台灣山薰香多樣性喪失分析 ... 67
表十三、未來氣候假設情境 CGCM2 氣候模型下，山薰香脆弱度指數程度之覆蓋 面積表 ... 68
表十四、未來氣候假設情境 HADCM3 氣候模型下，山薰香脆弱度指數程度之覆 蓋面積表 ... 69
表十五、未來氣候假設情境 CGCM2 氣候模型下，台灣山薰香脆弱度指數程度之 覆蓋面積表 ... 70
表十六、未來氣候假設情境 HADCM3 氣候模型下，台灣山薰香脆弱度指數程度 之覆蓋面積表 ... 71
v
附錄目錄
附錄一、標本資訊 ... 72
附錄二、本研究 R 環境下的示範操作--PD_DL.demo.R ... 86
附錄三、R code: PNCA() ... 92
附錄四、R code: PNCA_cutloop() ... 97
附錄五、R code: PNCA_plot() ... 98
附錄六、R code: geoDis_format() ... 103
附錄七、R code: geoDis_output() ... 111
1
前言
I.
氣候變遷對高山生態系的衝擊全球氣候變遷對於物種的影響如同溫水煮蛙，過程漸變而難以回復。近百年 的氣象資料顯示，全球氣溫增加了 0.6℃，而國內的氣象資料更指出臺灣百年來 增溫了 1.4℃（許晃雄等，2012），此暖化的趨勢同時伴隨極端氣候如旱澇、熱浪 或暴雨機率的增加，緩慢而加速地衝擊脆弱的生態系統。高山環境在地質歷史中 時常扮演避難所的角色，同時也是維護生物多樣性重要的棲地，研究顯示，全球 氣候變遷對於高山生態系的衝擊甚劇（Sala et al., 2000），例如阿爾卑斯山區的物
種的觀察結果中，在過去百年內物種最適生存海拔提高約 29 公尺（Lenoir et al.,
2008），依據海拔播遷理論（Ohsawa and Ide, 2008），當受到氣候變遷衝擊時，高 山物種可透過緯度或者海拔遷移，但緯度上的遷移可能受限地理屏障難以達成，
因此迫使族群傾向海拔上的遷移，然而高海拔頂峰區域的面積相對狹小，族群擴 張的範圍受限，區域的族群可能縮小，甚至低於最小活力族群而導致族群滅絕，
因此受到氣候變遷衝擊之下，高海拔嗜冷物種的處境將更加危急，必須急切了解。
檢視台灣的情況，依據台灣植物誌第二版紀錄，台灣維管束植物的物種多樣 性高達 4000 餘種，物種數目隨海拔高度遞減，但特有率卻隨之遞增，海拔 2500 公尺以上特有率高達 60％（Hsieh, 2002），顯見台灣高山物種在全球尺度的特殊 性，然而對於氣候變遷可能帶來的物種滅絕危機，台灣的高山物種將何去何從將 是亟欲關切的課題（李彥希、顏宏旭，2011），於台灣中部高山植物群落多樣性 研究中（Chou et al., 2011），觀察過去近三十年的植物多樣性變化，發現部分物種 有往高海拔遷移的趨勢，在全球增溫的條件下，甚至有滅絕的危機。
2
II.
脆弱度評估根據跨政府氣候變遷委員會（Intergovernmental Panel on Climate Change, IPCC）對脆弱度的定義（Parry et al., 2007; Williams et al., 2008），脆弱度為針對 一個系統的負向效應的感受性（susceptibility），其內含歸類為三要素：（一）暴露 度（exposure）：衝擊強度；（二）敏感度（sensitivity）：目標主體受到的衝擊程度；
（三）調適能力（adaptive capacity）：目標主體的適應能力，簡單以彈力球受到 碰撞回彈的例子比喻，彈力球代表欲評估脆弱程度的主體，暴露度即彈力球受到 的碰撞強度，敏感度及彈力球受到碰撞後彈向地面的速度，而適應能力即彈力球 本身的回彈高度。
在評估氣候變遷造成人為活動的系統中，如水資源、經濟活動等方面應用甚 多，各要素的界定及計量方法較為明確且成熟，然而應用在評估生物多樣性的自 然系統中，諸多要素的不確定性以及相互關係不明確，且於不同分類群間差別迥 異，使得不容易系統性評估自然系統。鑑於分子標記的引入，對於物種內遺傳訊 息的掌握提升，加上地理資訊系統的應用，增進物種空間特性以及棲地環境的理 解，並且整合生物學上的知識，Williams et al.（2008）提出整合性的架構，評估 氣候變遷衝擊下的物種脆弱度（圖一），基於最小化全球生物多樣性的喪失，進 行環境經營管理規劃之前，首先必須實行評估物種脆弱度以設定管理的優先權，
將脆弱度分為外在和內在因素兩部分：外在因素指的是暴露度，影響因子包含氣 候變遷帶來的區域性衝擊，以及微環境或地形的緩衝能力；內在因素指的是物種 的敏感度，與物種特性有關，包括生態、生理及遺傳多樣性，並將物種特性整合 為兩個評估指標：恢復力（resilience）和調適能力，恢復力包含生活史特性、族 群動態、傳播能力、空間分布情形；調適能力強調物種潛在的適應性，包含基因 多樣性、親緣地理多樣性（phylogeographic diversity）、以及可塑性。藉助脆弱度 評估推論潛在的衝擊，以便規劃、實行保育措施，達到降低氣候變遷對物種和其
3
棲地衝擊的目的，並監測實際的衝擊，偵測物種的後續狀態。
本研究參考 Bálint et al.（2011）、Pauls et al.（2013）及 Pfenninger et al.
（2012），整合族群遺傳分析方法及物種分布預測模型，實踐物種脆弱度的評估，
特別是在於這類研究多採用 AFLP 或微衛星資料進行評估，然而此方法為採用序 列資料的先例。以往族群遺傳研究長於遺傳多樣性的計量及推論演化歷史，進而 提出族群分化成因以及經營管理建議，但多拙於量化空間資訊，在此加入物種分 布預測模型，有效量化物種與環境因子的相關性，並且提供一快速的方法模擬物 種潛在分布，在氣候變遷衝擊的前提下，結合兩者獨立研究方法，強化了評估物 種脆弱度精準度。以 Bálint et al.（2011）的案例而言，評估歐洲九種水生昆蟲，
藉由判定不同種內變異，以代表種內多樣性的單位，結合物種分布預測模型進行 遺傳多樣性喪失評估，相較以形態差異為主的物種為單位的評估，解析度有明顯 的提升，且結果顯示面積的喪失與遺傳多樣性的喪失之間具有明顯的正相關性。
III. 物種分布預測模型
近年來地理資訊的資料公開與完備，使得以物種分布模型（Species Distribu- tion Model, SDM）預測物種分布可行，隨之為物種空間分布的假說驗證提供極佳 的研究平台。物種分布模型基於物種生境概念（niche concept, Hutchinson, 1957;
Soberon and Nakamura, 2009），透過數學模型結構建構物種的潛在分布。一般進 行建模的流程，首先將原始資料分為訓練子集及驗證子集，以訓練子集建立模型，
驗證子集用來評估模型優劣。然而模型結構所適用的資料形態及預測效力不一，
Elith et al.（2006）比較 16 個模型結構的預測效力，模型所適用的資料形態分為 出現/未出現（Present/Absent）或僅出現（Present-only）資料，預測效力結果前三 名為 MARS、BRT 及 MAXENT，其中 MAXENT 對於狹域分布物種預測結果佳，
4
且適用僅出現資料，加上具有方便的使用者介面便於參數調整，故諸多研究中採 用之。
MAXENT（Phillips et al., 2006）為僅出現模型，應用最大熵模型（maximum entropy model），取出物種分布與背景隨機點之環境變量，比較兩者環境變量空間，
在滿足一定條件下，運用機器學習（machine learning）方法估計最大熵的機率分 布，藉此建立物種預測分布（Elith et al., 2011）。
模型建立之後，以評估模型判定模型預測效力，首先建立描述實際分布與預 測分布之間正確與否的誤差矩陣（confusion matrix），分為四項：真肯定（true positive）、假肯定（false positive）、真否定（true negative）、及假否定（false neg- ative），接著建立操作反應特徵（receiver operating characteristics, ROC）曲線下面 積（area under curve, AUC）、max Kappa 或 TSS 等評估指標判斷模型預測準確程 度，其中 AUC 值對於資料數量的獨立特性，最常用為判斷基礎，AUC 大於 0.5 表示具預測能力，AUC 值介於 0.5 與 0.7 之間表示預測能力差，介於 0.7 與 0.9 之 間表示預測能力佳，大於 0.9 表示預測能力極佳（McClish, 1989）。
接著依各研究目標所需，需選定閾值將機率預測分布轉為二元預測分布，Liu
et al.（2005）比較 12 種閾值選取方法，以 sensitivity、specificity、overall predic-
tion success 及 Kappa 指標評定優劣，表現較佳的為 prevalence、average predicted probability/suitability 、 sensitivity-specificity sum maximization 、 sensitivity- specificity equality 和 the shortest distance to the top-left corner in ROC plot。選取適 當的閾值有助於以進行面積計算及後續評估分析。IV. 台灣高山植物之親緣地理
在國內探討高山物種的遺傳多樣性及其演化歷史的相關研究，僅有兩個木本
5
植物類群：玉山杜鵑複合種群、台灣冷杉，以及一個草本植物：碎雪草屬，此少 數的親緣地理的研究推論並不一致。
玉山杜鵑（Rhododendron pseudochrysanthum Hayata）複合種群親緣地理的研 究中（黃啟俊, 2006; Chung et al., 2007），選用葉綠體 DNA atpB-rbcL 基因間片段 及核基因 ITS 片段，推測起源於紅星杜鵑，推斷族群在冰河期由北向南遷移而廣 泛分布於中、低海拔區域， 隨後間冰期增溫棲地逐漸破碎種化。
台灣冷杉（Abies kawakamii (Hayata) Ito）的族群遺傳研究（Shih et al., 2007）
檢視三個葉綠體 DNA 片段：trnR-trnG、rbcL、及 trnL-trnF 以及核基因 GapC 內 插子片段，葉綠體的片段中偵測到族群幾乎無變異存在，於核基因 gapC 偵測到 廣布的基因單型，其中雪山族群與其他族群差異較大，但無明顯的族群遺傳結構 存 在 ， 由 孢 粉 證 據 顯 示 台 灣 冷 杉 過 去 歷 史 分 布 可 能 擴 及 現 今 中 海 拔 的 範 圍
（Tsukada, 1966），現今分布相較退縮至高海拔區域，設想與冰河歷史可能有關連，
然而從遺傳資訊推論並無顯著的選汰壓力。
玉山碎雪草（Euphrasia transmorrisonensis Hayata）複合種群的研究中（吳明
洲, 2004; Wu et al., 2005），選用葉綠體 DNA trnL 內插子及 trnL-trnF 片段以及核 基因 ITS 片段，結果呈現明顯的族群遺傳結構，發現了三個支序：北部、北部和 東部、以及中部和南部區域，且於核基因 ITS 中發現在北部及東部支序與其餘兩 者有部分族群雜交的現象，並且推論是由於冰河期的壓力，使得族群由較低海拔 的北部和東部支序向高海拔擴展，其中均衡選汰可能扮演重要的演化動力
V.
台灣山薰香屬簡介山薰香屬（Oreomyrrhis Endl.）為分布在環南太平洋的植物，全世界約 25 種，
其分布在東太平洋區域北起中美洲的山區，往南沿南美洲安地斯山脈分布至阿根
6
廷北部以及智利南端及福克蘭群島，大洋洲區域的類群分布於紐西蘭、澳洲東南 沿海區域、塔斯瑪尼亞島及新幾內亞山區，在西太平洋區域分布於婆羅洲北端的 神山及台灣的高山區域。其分布的氣候區屬於亞高山、高山或亞極地，且常見於 森 林 底 層 開 闊 潮 濕 環 境 或 山 頂 岩 屑 地 的 棲 地 類 型 ， 為 典 型 的 高 山 植 物 類 群
（Mathias and Constance 1955; Chung et al., 2005）。Chung（2007）以 DNA 序列指 出本屬為一單系群（monophyly），由細葉芹屬（Chaerophyllum L.）演化出，因此 維持山薰香屬的地位將使細葉芹屬成為並系群（paraphyly），故 Chung（2007）將 山薰香屬（Oreomyrrhis）併入細葉芹屬。本論文因此仍以山薰香屬稱之。
Chung et al.（2005）及 Chung（2007）研究指出山薰香屬現今環南太平洋的 不連續分布為近代長距離傳播事件所生成的地理分布，根據分子定年，山薰香屬 起源於上新世至更新世間的 Sahul 大陸（約為今日新幾內亞和澳洲大陸），經歷一 次獨立事件拓殖到紐西蘭，以及至少兩次長距離傳播事件到達美洲、新幾內亞及 台灣（Chung, 2006）。
台灣的山薰香屬植物包含三個物種：山薰香（Chaerophyllum involucratum (Hayata) K. F. Chung）、台灣山薰香（Chaerophyllum taiwanianum (Masam.) K. F.
Chung）、以及南湖山薰香（Chaerophyllum nanhuense (Chih H. Chen & J. C. Wang) K. F. Chung），三種皆分布於台灣島內海拔 2500 公尺以上的山地區域，涵蓋雪山 山脈（雪山與大霸尖山）、中央山脈（南湖大山、合歡山、奇萊山、能高山、秀 姑巒山、向陽山與北大武山）及玉山山脈（玉山、鹿林山、塔關山與關山嶺山）
等主要山系，山薰香為廣布種，台灣山薰香侷限分布於北部及中部山系，南湖山 薰香則僅出現於南湖大山。然而台灣山脈多南北走向，河川東西向切割山脈，地 形突起度高，因此各族群間呈間斷分布。
7
研究目的
(1) 探討台灣山薰香屬植物之族群遺傳多樣性的現狀及形成原因。
(2) 建立台灣山薰香屬之物種分布，並預測氣候變遷衝擊之物種分布。
(3) 評估台灣山薰香屬之族群遺傳變異在全球氣候變遷衝擊的影響程度。
8
材料與方法
本研究方法流程參考圖二。
I. 研究材料
研究材料蒐集分為兩部分，一為植物體採集及遺傳資料蒐集，二為環境資料 取得。
A. 生物資料取得
山薰香屬植物的資料大多為鍾國芳老師於 2002 年至 2009 年採集之樣本，部 分樣本取自標本館標本及個人贈與，其餘為作者 2011 年至 2012 年採集，相關標 本資訊參考附錄一。
1. 野外取樣
自生育地採集植物組織置入乾燥劑中乾燥保存，以備進行遺傳分析，另外記 錄採集資訊，包含物種中文名稱、學名、採集者、採集者標本編號、採集地點資 訊、經緯度座標。物種點位資料除野外取樣之外，部分來自於全球生物多樣性資 訊機構（global biodiversity information facility, GBIF），並檢核摒除位置明顯錯誤 的點位。
2. DNA 萃取、擴增及定序
a. DNA 萃取
採用修改的 CTAB 法（Murray and Thompson, 1980）進行 DNA 萃取，操作 步驟如下：
(1) 將 cetrimonium bromide（CTAB）置於 65˚C 水浴槽預熱。
(2) 秤取 0.1g 的乾燥植物組織，加入少許滅菌過的海砂，於研缽中磨碎成
9
粉末狀並置於 2 mL 離心管中。
(3) 加入 4μL β-mercaptoethanal、10μL RNase A 及 800μL 65˚C 的 CTAB，
以試管震盪器震盪離心管，使組織粉末與溶劑均勻混合。
(4) 將離心管置於 65˚C 水浴槽 1 小時，每 15 分鐘將離心管取出搖勻再置 回。
(5) 將離心管移出水浴槽，加入 800μL 的 Chloroform : Isoamyl alcohol 24:1
（CIAA）。
(6) 劇烈搖晃離心管 3 分鐘，使管內溶液呈乳糜狀。
(7) 將離心管置於離心機以 13,500 g 轉速離心 10 分鐘。
(8) 小心取出離心管，以微量滴管移出上清液至新的 1.5mL 離心管，加入 700mL CIAA。
(9) 劇烈搖晃離心管 1 分鐘，再置於離心機以 17,000 g 轉速離心 10 分鐘。
(10) 移出上清液至新的 1.5 mL 離心管，加入 1/10 體積的 3 M Sodium Ace- tate 及 0.54 倍總體積的-20˚C 的異丙醇（isopropenol）。
(11) 輕輕搖勻之後冷藏於-20˚C 2-10 小時。
(12) 以 17,000g 轉速離心 15 分鐘，移去上清液，以 75％酒精清洗沈澱物。
(13) 以 17,000g 轉速離心 1 分鐘，移去上清液，加入 800mL -20˚C 的 95％
酒精清洗沈澱物。
(14) 以 17,000g 轉速離心 1 分鐘，移去上清液，以 800mL -20℃的 99.99％
酒精清洗沈澱物。
(15) 加入 30μL 的 TE buffer 回溶。
(16) 取其中 2μL 進行 1％瓊脂膠（agarose）電泳分析，以 DNA 標準品作 為估計樣品 DNA 的濃度。
b. 聚合酶連鎖反應（polymerase chain reaction, PCR）
本研究取樣葉綠體 atpB 與 rbcL 兩基因片段之間的片段（~740 bp），使用引
10
子如表一。聚合酶連鎖反應於 200μL 離心管中進行，總體積為 25μL，組成包 括：2.5μL 稀釋 10 倍之樣品模版 DNA、2.5μL atpB 端引子、2.5μL rbcL 端引 子、12.5μL 2X Taq DNA Polymerase Master Mix RED（Ampliqon）及 5μL 二次 蒸餾水，溫度循環流程設定參考 Chung（2006）：94˚C（5 分鐘）；30 次循環：
94˚C（30 秒）、54˚C（1 分鐘）、68˚C（30 秒）；68˚C（7 分鐘）；25˚C（1 分鐘）；
4˚C（暫停）。
為確認 PCR 之產物生成，取 2μL 產物與微量染液混合後進行 1％瓊脂膠電 泳分析。
c. 純化及定序
取 10μL PCR 產物送至基龍米克斯生物科技有限公司（系統：ABI 3730XL），
填妥表單完成純化及定序，樣本的定序結果將以*.ab1 格式回傳。
B. 環境資料取得
本研究的環境資料來自於線上的資料庫，每個資料集（表二）包含 19 個生 物氣候圖層（bioclim layers）。資料集的時間尺度包括最近一次冰河極盛期（Last glacial maximum, LGM）、現今、以及溫室氣體排放情境特別報告（Special Report on Emissions Scenarios, SRES）預測的未來氣候假設情境，前兩個時間區段資料集 來源為 WorldClim（http://www.worldclim.org/, Hijmans et al., 2005），其中 LGM 包 含兩個古氣候模型：Community Climate System Model（CCSM）與 Model for In- terdisciplinary Research on Climate（MIROC3.2）；SRES 未來氣候假設情境資料集 來 自 Centro Internacional de Agricultura Tropical （ CIAT ）（ http://www.ccafs-
climate.org/, Ramirez and Jarvis, 2008），選用兩氣候模型 The Second Generation Coupled Global Climate Model（CGCM2）及 Hadley Centre Coupled Model, version 3（HADCM3）在前期（2020 年）、中期（2050 年）與後期（2080 年）三個時間
11
區段的預測。，各時間區段根據溫室氣體排放程度又分為 A2A 和 B2A 兩個假設 情境。總結本研究所採用的資料集，包含 LGM（2 個）、現今（1 個）、SRES 未 來假設情境（2×3×2＝12 個）（圖三），共計使用 15 個資料集，資料解析度為 30 秒，其中 LGM 資料來源提供的解析度為 2.5 分，為便於投影，將之重新取樣為 30 秒（Waltari et al., 2007），而針對全球尺度的預測，考量運算效率採用 2.5 分的 資料解析度。
為方便裁切出研究範圍的環境因子資料層，同時在 WorldClim 下載國家面向 量（polygon）圖層，使用 ArcGis 10（ESRI）的 extract by attribute 工具，將台灣 範圍的面向量抽出，另存成 shapefile 資料層。
II.
資料分析資料分析階段分為四個部分：資料前處理、遺傳分析、物種分布預測模型建 立及多樣性喪失分析，本研究提供 R 環境下的示範（PD_DL_demo.R），參考附錄 二。
A. 物種及環境資料前處理
將 DNA 定序的資料（*.ab1 格式）以 SeqMan ver. 7.1 進行序列判讀並匯出
*.fasta 格式，匯入 MEGA 5（Tamura et al., 2011）進行序列排序（sequence align- ment），並以手動調整序列資料輸出*.nexus 格式。另外，將樣本的地理分群參考 採集資訊輔以經緯度座標判定，製成樣本代號與地理分群的成對表單。
環境資料集以 R 軟體（R Core Team, 2013）中的 dismo、sp、raster 套件處理，
由於儲存空間考量，線上資料庫的資料層採用整數數值，因此參考 WorldClim 網 站的資料格式解釋處理，將部分資料層回復到正確的數值，例如：年均溫（AMT）
原始資料層值域為-497 至 291，需除以 10 轉換為-49.7 至 29.1˚C。接著將資料層
12
裁切成研究範圍的區域，抽取出台灣的面向量圖層，以 crop、mask 兩個函數裁切，
最後，將一個資料集的資料層儲存於同一資料夾，以備後續建模使用。
B. 族群遺傳變異分析
以 Arlequin ver. 3.5.1.3（Excoffier et al., 2005）軟體建立族群之基礎遺傳資訊，
包含基因單型歧異度（haplotype diversity, Hd）及核苷酸歧異度（nucleotide diver- sity, π），並分別估計其標準差，藉此比較各族群的遺傳歧異程度。再者，以 Tajima’s D 及 Fu’s Fs 檢測族群之中性變異假說，並進行 1000 次溯祖模擬計算顯 著性，檢測族群是否偏離隨機基因漂變，藉此推論族群消長歷史，負值代表可能 曾經歷瓶頸效應（bottleneck effect）之後族群擴張或定向選汰（directional selec- tion），正值代表族群可能曾經經歷族群退縮或定向選汰（balancing selection），而 Fu’s Fs 相較 Tajima’s D 對於偵測族群擴張事件較為敏感。最後，建立族群間族群 分化指標（FST）矩陣，並將數值區分為三個等級：低度分化（FST＜0.15），中度 分化（0.15＜FST≦0.25）及高度分化（FST＞0.25），藉此推斷族群間的遺傳相異程 度。
C. 巢狀支序親緣地理分析（Nested Clade Phylogeographic Analysis, NCPA）
以 TCS 1.21（Clement et al., 2000）軟體建構統計最簡約網狀圖（statistical
parsimony networks ）之基因單型親緣網狀圖，依據 Pfenninger and Posada（2002）
的方法移除迴圈，再以 GeoDis 軟體配合推論檢索表（inference key）（Posada et
al., 2000）進一步進行 clade distance（Dc）與 nested clade distance（Dn）關連性
分析與推論族群演化歷史。上述操作中由於軟體所需的資料格式特殊，對於資料 彙整、資料格式建立及製圖，作者在 R 的環境下建立了數個函數以自動化產生：(1) PNCA（附錄三）：將樣本採集資訊、樣本地理分群和 TCS 1.21 產生的
*.graph 及*.graph.log 檔中的資訊整合。
13
(2) PNCA_plot（附錄四）：繪製 TCS 1.21 建構的基因單型親緣網狀圖，端 點（即各基因單型）以所含樣本所在地理分區的圓餅圖表示。
(3) PNCA_cutloop（附錄五）：處理基因單型親緣網狀圖中的迴圈，以便 GeoDis 分析。
(4) geoDis_format（附錄六）：以問答的方式建立 GeoDis 軟體所需的輸入 資料格式。
(5) geoDis_output（附錄七）：將 GeoDis 軟體的分析結果製成利於進行推 論檢索表的輸出格式。
親緣網狀階層分析結果以基因單型親緣網狀圖、基因單型地理空間分布及階 層推論分析表呈現。
D. 物種分布預測模型建立
1. 環境因子選取
整合相關性檢測與 Maxent 3.3.3e（Phillips et al., 2006）內建的環境因子檢測，
以及觀察空間預測結果，選取環境因子組合以建立模型。首先，抽取樣本座標的 現今環境資料層，建立皮爾遜積差相關係數（Pearson's product moment correlation coefficient）矩陣，選取相關性小於 0.8 的環境因子組合，降低環境因子間共線性 的效應（Dormann et al., 2013），避免模型建立時產生過度預測（overestimate）的 情形。採用 R 環境中 cor 函數進行分析。另一部份，藉助 Maxent 3.3.3e 提供的因 子檢測，建立 Jacknife 因子重要性檢測，以及因子反應曲線，前者藉由每次刪去 或保留一個因子觀察模型預測結果，以評估該因子的重要性；後者計算因子邊際 機率分布，便於解釋因子間的相關性。
2. MAXENT 模型建立、評估模型與預測分布
將物種的點位及現今環境資料集匯入 Maxent 3.3.3e，勾選選出之環境因子組
14
合，勾選 create response curves、Do jackknift to measure variable importance，並輸 入檔案輸出路徑，並指定欲投影之時間尺度的資料集資料夾路徑，參數設定
（Settings）：
(1) 勾選 Random seed
(2) Random test percentage 輸入 20，代表隨機抽取百分之二十的物種點位 作為驗證子集，以進行評估模型，其餘為訓練子集，作為建立模型用。
(3) Replicates 輸入 100，代表重複建立模型 100 次，以利後續進行預測模 型挑選。
(4) Replicated run type 選擇 Bootstrap。
(5) 勾選 Write plot data 及 Write background predictions，方便匯入不同軟體 介面使用。
設定完成，模型建立、評估模型與預測分布將自動進行，其摘要表單儲存在 所設定的輸出路徑的 maxentResults.csv 檔案中，包含模型資訊、評估模型 ROC 曲線之 AUC 值及閾值。
在 R 環境下，挑選 100 個預測模型中 Test AUC 值前十名的模型取平均作為 平均機率預測分布，以降低模型的不確定性（Torres et al., 2013），再選擇閾值以
Equal training sensitivity and specificity logistic threshold（Liu et al., 2005）將各時 間尺度之平均機率預測分布重新分類為有（1）、無（0）的二元預測分布。
E. 多樣性評估分析 1. 分布面積動態
將三個時間尺度：LGM、現今和 SRES 未來氣候假設情境 A2A 及 B2A 之二 元預測分布繪圖，並計算相對於現今的預測分布網格面積差異與差異之比例，其 差 異 包 含 穩 定 （ stable ） 區 域 ： 現 今 及 投 影 預 測 分 布 皆 出 現 的 區 域 、 擴 張
15
（expanded）區域：僅出現於投影預測分布的區域，於現今預測分布未出現的區 域、以及縮減（restricted）區域：僅出現於現今預測分布的區域，於投影預測分 布中喪失的區域。
差異比例 投影之二元預測分布面積 現今二元預測分布面積
現今二元預測分布面積
2. 基因單型喪失評估
在各時間尺度的預測分布中，計算於二元預測分布之外的樣本採樣點，統計 喪失的點位以及基因單型，並且計算樣本採樣點以及基因單型喪失的比例。
另一部份，結合親緣網狀階層分析，將喪失的樣本點位納入考慮，重新計算 基因單型親緣網狀圖中各端點圓餅圖的比例，若喪失某一基因單型，則與之相關 的連線便中斷。
F. 脆弱度評估
1. 脆弱度指數（Vulnerability Index, VI）
本研究將藉由脆弱度指數判斷物種潛在易受害的棲地，藉由各時期的機率分 布轉化為脆弱度指數，觀察其出現機率的變化特性，如變化幅度以及空間分布，
用以預測各族群的脆弱程度。實際的施行方法為：將脆弱度指數定義為現今平均 機率分布與各時間尺度之平均機率分布之比值（Matsui et al., 2004; Casalegno et
al., 2010; 孫世鐸, 2011）
，並且重新分類為五個等級加以顏色輔助辨識：安全（0＜VI≦1，綠色）、低度脆弱（1＜VI≦5，藍色）、中度脆弱（5＜VI≦10，黃色）、
高度脆弱（10＜VI≦100，紅色）、以及極端脆弱（100＜VI，紫色），最後在地理 空間上視覺化呈現脆弱度指數的空間分布圖，並且計算面積及出現點在各脆弱度 等級的分布情形。
16
脆弱度指數 現今平均機率分布
各時間尺度之平均機率分布
2. 海拔梯度機率頻度分布
抽取各時間尺度之平均機率預測分布之海拔資料，以每 200 公尺為一級距，
計算各海拔級距上的預測分布機率累積製成直方圖。
3. 棲地破碎程度分析
根據棲地破碎化指標（Bogaert et al., 1999），計算各時間尺度下的內緣比
（interior-edge ratio），在此定義為綴塊（patch）的面積與邊長之比值，藉此評估 棲地破碎化的程度。實際作法則將二元預測分布網格資料轉為面向量圖層，統計 其面積及邊長，並將面積除以邊長累積製成直方圖。
內緣比 綴塊面積
綴塊邊長
17
結果
I. 族群遺傳變異分析、親緣網狀階層分析
此部分分析描述，僅山薰香（Chaerophyllum involucratum）進行 clade dis- tance 與 nested clade distance 關連性分析。
一共取得 269 條長度為 753 個鹼基對之葉綠體 atpB-rbcL 基因間片段序列，
包含山薰香 216 條及台灣山薰香 53 條。
山薰香的資料矩陣中偵測到 9 個多型性位點（polymorphic sites），並區分出 19 個基因單型。樣本在空間上分布於 8 個地理分區（圖四、五）：雪霸（Sp）、南 湖（Nh）、合歡（Hh）、能高（Nk）、秀姑巒（Sgl）、玉山（Ys）、塔關（Tk）及 北大武（Ptw），其中除南湖沒有獨有的基因單型外，其他 7 個地理分區皆具有獨 有的基因單型，其獨有基因單型數目由多至少依序為雪霸（3），秀姑巒、玉山和 塔關次之（2），合歡、能高及北大武各 1 個。此外，7 個共享基因單型各自在地 理空間中不具明顯的空間群聚性。台灣山薰香的資料矩陣中，偵測到 7 個多態位 點，並區分出 5 個基因單型。樣本在空間上分布於 4 個地理分區（圖六）：雪霸、
南湖、玉山及秀姑巒，獨有基因單型在南湖有 2 個及玉山 1 個，2 個共享基因單 型在空間中各分別聚集於北部和中部兩部分（圖七、表三）。
山薰香整體的基因單型歧異度（Hd）為 0.916（s.d.=0.007），各地理分區則介 於 0.399 至 0.788 之間，雪霸最高，南湖次之，最低者為塔關；整體的核苷酸歧 異度（π）為 6.223×10-3（s.d.=3.386×10-3），各地理分區則介於 1.366×10-3 至 9.909×10-3之間，南湖最高，合歡次之，最低者為能高（表四）。
中性假說檢測部分（表四），Tajima’s D 及 Fu’s Fs 不論是族群整體或各地理 分區幾乎都為正值，推論族群可能族群退縮或均衡選汰，但仍有少數地理分區的 族群 Tajima’s D 為負值，表示部分族群可能經歷瓶頸效應之後族群擴張或定向選
18
汰，然而統計上 Tajima’s D 和 Fu’s Fs 皆不顯著，無法推翻中性變異假說，顯示可 能有其他事件發生，如族群二次接觸或族群遷移。
族群間的遺傳分化指數（FST）部分（表五），除了合歡-南湖及北大武-秀姑 巒兩組呈現低度分化，其餘族群間遺傳分化皆為中度或高度分化。
山薰香親緣網狀關係（圖五）部分，整體並無明顯輻射狀結構，TCS 判定的 祖先型為 I，但並非廣布的基因單型。親緣網狀階層分析（表六）總共包含 10 個 1-step clade、4 個 2-step clade 及 2 個 3-step clade，除 clade 1-6 與地理位置相關性 不 顯 著 外 ， 其 餘 支 系 皆 與 地 理 位 置 相 關 。 整 體 推 論 訊 息 呈 現 族 群 連 續 擴 張
（contiguous range expansion），1-step clade 和 2-step clade 出現較多異域隔離
（allopatric fragmentation）及歷史隔離或長距離拓殖（past fragmentation and/or long distance colonization），少數推論訊息為族群連續擴張、基因交流受限或族群 擴散但有部分長距離傳播事件（restricted gene flow/dispersal but with some long distance dispersal）以及距離隔離事件導致基因交流受限（restricted gene flow with isolation by distance）； 3-step clade 的推論訊息皆為族群連續擴張。
II.
物種分布預測模型時間序列A. 環境因子選取
三種台灣山薰香屬的 19 個氣候環境因子如表七，符合典型高山物種的特性。
綜合皮爾遜積差相關係數矩陣（表八、九）、Jacknife 因子重要性檢測和因子反應 曲線，山薰香及台灣山薰香兩組分類群各選出一組環境因子組合，山薰香：年均 溫（AMT）、年降水量（AP）、最乾月降水量（PDM）和最暖四分之一年降水量
（PDrQ）；台灣山薰香：年均溫、月平均氣溫日較差（MDR）和最暖四分之一年 降水量。
19
B. 預測分布模型
藉由選取出的環境因子組合建立模型，以 AUC 值及空間中預測情況判斷模 型優劣，結果如表十，每個模型皆呈現出良好的預測結果（山薰香：0.9877 和台 灣山薰香：0.9995），其機率預測分布反應高山地區機率值越高的趨勢（圖八）。
選擇 Equal training sensitivity and specificity logistic threshold 將 Test AUC 值前 十名機率預測分布處理為二元預測分布，其結果反應真實分布區域，譬如陽明山 區並無分布而北大武山實際有分布。然而台灣尺度中，唯中央山脈中部丹大山區 沒有樣點分布，但卻被預測為有分布區域。
在 LGM 的二元預測分布中，CCSM 及 MIROC3.2 的結果皆呈現分布擴大的 趨勢，但前者較後者預測範圍大了許多，在山薰香的預測中尤為明顯，甚至預測 部分中北部平原地區也有分布。
SRES 未 來 的 二 元 預 測 分 布 中 ， 隨 著 時 序 由 前 期 至 後 期 ， CGCM2 和 HADCM3 皆呈現面積朝向高山集中的現象，整體面積明顯減少，而 CGCM2 較 HADCM3 更甚之，在 A2A 及 B2A 比較上，前者退縮情形較嚴重，減少的面積越 往後期越顯著。
III. 多樣性評估分析
A. 分布面積動態
於山薰香的分布動態結果中（表十一），LGM 的兩組氣候模型預測結果趨勢 類似，100％為穩定區域，126.7-500.9％為擴張區域，並無退縮區域，整體顯示 LGM 時 族 群 擴 張 情 形 。 預 測 未 來 分 布 部 分 ， 整 體 上 CGCM2 改 變 幅 度 較 HADCM3 明顯，而氣候假設情境間的比較，A2A 較 B2A 反應強烈許多。穩定區 域由前期 56.8-88.6％，及至中期降至 29.5-75％，於後期更降到 14.2-55.3％；擴 張區域部分，大部分為 0％，僅少數時間區段預測出擴大的現象，但幅度都不大
20
（0.3-5.3％）；縮減區域在時序上比例逐漸增加，前期 11.4-43.2％，至中期升至 25-70.5％，於後期達到 44.7-85.8％。
台灣山薰香的分布動態部分（表十二），於 LGM 預測結果與山薰香的情況相 同，同樣顯示 LGM 時族群擴張情況。預測未來預測分布部分，綜觀下，資料集 的比較，CGCM2 改變幅度較 HADCM3 明顯，而氣候假設情境間的比較，A2A 較 B2A 反應強烈許多，整體的趨勢與山薰香相似，但分布面積更加狹小，變動程 度更顯著。穩定區域由前期 44.5-100％，到中期降至 6.2-69.6％，於後期甚至降到 0.1-50.8％；擴張區域部分，皆為 0％，除前期 A2A CGCM2 預測到分布面積增加 的情形（141.2％），及前期 A2A 和 B2A HADCM3 預測到輕微的面積增加現象
（3％及 2.5％）；縮減區域趨勢明顯增加，由前期 0-55.5％，至中期升至 30.4-93.8
％，及至後期達到 49.2-99.9％。
B. 遺傳多樣性喪失評估
根據前述未來氣候假設情境時棲地面積縮減情況，判定此分類群為潛在退縮 的類群，因此進一步評估基因單型喪失的程度，以下依計數方法論述，包含點位 座標及基因單型兩部分。
山薰香的結果中（表十一），點位減少比例介於 0-37.5％，並隨著時間增高，
兩組氣候模型差值約 3.2-30.6％，CGCM2 減少比例高於 HADCM3，中後期尤其 明顯。基因單型喪失較點位減少的趨勢緩和，喪失比例介於 0-21.1％，在時間梯 度上越往後期喪失比例越高，並且觀察喪失的基因單型於親緣網狀關係圖的位置 位於末端。
台灣山薰香部分（表十二），點位減少比例介於 0-83％，隨著時間增高，兩 氣候模型差值約-17-54.7％，在 A2A CGCM2 喪失反應強烈。基因單型喪失部分 僅 A2A CGCM2 偵測到基因單型減損，於中期減損 40％，至後期減損到達 80％。
第二部分各地理分區的減少比例中，採集座標與基因單型計量的結果一致，
21
山薰香的喪失比例大致有三個高峰，北部雪霸族群、中部合歡及能高族群以及南 部北大武族群（圖九、十）；台灣山薰香的喪失比例以中部秀姑巒及玉山族群可 能受到較大威脅（圖十一、十二）。
IV. 脆弱度評估
A. 脆弱度指數
觀察山薰香的脆弱度指數空間分布圖（圖十三），CGCM2 較 HADCM3 反應 明顯，A2A 又較 B2A 反應明顯，整體數值分布模式皆由西部向東北和東部遞增，
且分布邊緣脆弱度較高。計量各等級脆弱度指數如表十三、十四，A2A CGCM2 脆弱度分布於暖化前期 81.78％集中於低度脆弱，少部分落在極端脆弱的區域，
及至中期，極端脆弱的區域增至 44.06％，於後期逾半數（63.53％）落於極端脆 弱的區域，統計取樣座標處於脆弱度各等級數量，同樣呈現此趨勢，於後期取樣 座標開始出現於極端脆弱的區域。其餘情境結果趨勢雷同但較為和緩。
台灣山薰香部分，其脆弱度指數空間分布與山薰香呈現趨勢相似但反應更劇 烈（圖十四），CGCM2 較 HADCM3 反應明顯，A2A 也較 B2A 反應明顯，整體 數值分布同樣為西部向東北和東部遞增，分布邊緣的脆弱度高。由表十五、十六 統計各等級脆弱度，觀察 A2A CGCM2 的結果，暖化前期皆處於安全的區域，及 至中期脆弱度分布快速升高，高達 82.96％落在極端脆弱的區域，於後期落在極 端脆弱的區域更升高至 97.37％，統計取樣座標處於脆弱度各等級數量，同樣呈 現此趨勢，於中及後期開始出現於極端脆弱的區域。其餘情境結果趨勢相似但較 和緩，整體脆弱度數值分布較山薰香的來得高。
B. 海拔梯度機率頻度分布
山薰香的現今最適分布海拔約為 2657±295.64 公尺（圖十五），CGCM2 的結 果顯示，最適海拔持續上升，A2A 情境：於前期為 2735±269.56 公尺，中期至
22
3022±192.99 公尺，後期到達 3154±182.48 公尺；B2A 情境：於前期為 2848±
242.09 公尺， 中期降至 2770 ±274.5 公尺，後期升 高到 2996 ±196.86 公尺。
HADCM3 部分，最適海拔同 CGCM2 持續上升的趨勢但較為緩和，A2A 情境：
於前期為 2692±286.66 公尺，中期至 2828±252.99 公尺，後期升至 2944±217.25 公 尺；B2A 情境：前期為 2675±296.81 公尺，中期升至 2764±259.74 公尺，後期到 達 2855±238.64 公尺。
台灣山薰香的現今最適分布海拔約為 2936±252.77 公尺（圖十六），CGCM2 的結果呈現最適海拔上生的趨勢，A2A 情境：雖然前期預測結果下降至 2572±
462.02 公尺，中期升高到 3412±134.1 公尺，後期更到達 3719 公尺；B2A 情境：
前期為 3111±204.48 公尺，中期稍微升高至 3117±204.73 公尺，後期到達 3195±
189.8 公尺。HADCM3 部分，最適海拔同 CGCM2 持續上生的趨勢但較為緩和，
A2A 情境：前期為 2953±241.75 公尺，中期升至 3063±211.21 公尺，後期到達 3114±207.27 公尺；B2A 情境：前期為 2955±241.09 公尺，中期升至 3028±220.96 公尺，後期到達 3079±213.68 公尺。
C. 棲地破碎程度分析
藉助計算內緣比以描述棲地破碎化情形，觀察山薰香及台灣山薰香不同時間 尺度下的直方圖結果（圖十七、十八），由暖化前期至後期，高比值有增多的趨 勢，且低比值的頻率減小，顯示棲地破碎程度加劇並且有部分小棲地喪失。
23
討論
I. 台灣山薰香屬遺傳多樣性成因
由族群遺傳變異分析結果顯示，各地理分區族群的遺傳多樣性程度高，而在 空間中並未發現廣布於全分布區且高頻率的基因單型，於親緣網狀關係圖中雖大 致呈現放射狀的親緣關係，但位於中心可能暗示為祖先型的基因單型則由散於數 個山頭的基因單型組成，結合地理空間分布及親緣關係圖，山薰香具有多起源的 演化歷史，進一步檢視親緣網狀階層分析，越高階層反應的歷史事件皆為族群連 續擴張，顯示山薰香族群於冰河期擴張的可能性，在物種預測分布中同樣支持此 論述，然而初階層反應較近期的歷史事件則有不同的推論，但多呈現族群交流受 限制而有隔離的情況發生，但未排除長距離傳播的可能性，另外由族群間分化指 數結果同樣顯示各族群間的分化程度相當高，反應出族群間遺傳多樣性特殊的性 質。
依據間冰期避難所模型（displacement refugia model, Kropf et al., 2002;
DeChaine and Martin, 2004，圖十九）敘述，相較於分布在低海拔與大陸平地的物 種，生活在高海拔環境的嗜冷物種對冰河週期的溫度變化可藉由在海拔梯度的上 下遷移調適，為維持原先的棲地溫度環境，當冰河期來時，山區的溫度帶向下移 動，嗜冷物種將自原先的山頂棲地向下遷移至山麓，族群可能藉由山麓連接而產 生基因交流，而當冰河退卻進入至間冰期，族群又退卻至山頂，族群間的基因交 流因此中斷，直到冰河期再次來臨，族群又再度相遇。受到如此反覆冰河週期的 影響，進而形成現今高山物種的族群遺傳多樣性及結構。此假說可解釋台灣的山 薰香屬植物親緣關係結構與地理分布間無明顯關連的現象，以及族群間遺傳高程 度分化的情形。
24
II. 物種分布模型
構築於 SDM 之上物種評估，因此預測模型將影響到評估的結果，其影響參 數為（一）資料誤差，包含樣本點及環境資料；（二）模型結構。
A. 資料誤差
台灣山地區域的氣候測站相對較少，環境資料來源為空間內差的結果，使得 資料本身具有空間自相關，山頭間的氣候梯度在 30 秒空間解析度中較難已呈現，
直接影響預測分布細部資訊的判斷，而因子間共線性同樣可能造成模型過度預測 的結果，然而同時不能忽略共線性高的因子可能也是重要的影響因子，因此考慮 因子相關性時盡量降低共線性，但也需考量各因子的反應曲線，以避免刪去重要 的環境因子。在樣本座標取樣上，觀察到部分預測出現區域無採集記錄，如中央 山脈中段丹大山區，原因可能為預測誤差，同時可注意到預測模型具有提出潛在 棲地的能力（Guisan et al., 2006），回饋以取樣上的誤差，並且若能記錄未出現的 資料將對於需要出現/未出現資料的模型結構提供更為精確的資訊。
B. 模型結構
本研究採用 MAXENT 建立預測分布模型，在樣本數小或狹域種時具有不錯 的預測結果（Franklin and Miller, 2009），然而在不同的模型結構，譬如 GLM、
GAM、RF、MARS 等等其他模型結構可能各自反應物種分布與棲地環境之間的 對應關係，整合多種模型結構有助於提高預測分布的準確率，但同時需要解決模 型間不確定性（Marmion et al., 2009）。
綜合以上兩部分參數討論，本研究建立之物種分布模型可視為物種分布界線，
預測結果與實際觀察的物種分布範圍吻合，已足以進行推論，作者初步比較與廣 義線性模型（generalized linear model, GLM）的結果，兩者相似，但 GLM 過於保 守而高估預測分布範圍，以本研究中所採用的 MAXENT 有較好的預測表現，故
25
並未整合多個模型的預測結果。
然而相較前人研究（Marske et al., 2009）的空間尺度，以 30 秒解析度的資料 在台灣尺度中建立預測模型，相對空間解析度略顯不足，以物種分布概念中所提 到的三項影響因子：環境因子、種間交互作用及可達性，受限於資料涵蓋的範圍 較大，本研究僅考慮環境因子中的氣候因素，若能考量特定的物種之其他環境因 子，如土壤、季風，或種間關係，如植群型、土壤微生物，將可增加物種預測分 布的準確度（Sinclair et al., 2010）。另一方面，尺度問題可能受限於實際環境因子 取樣密度，後續經空間內插補足未調查區域，特別是山地區域，以台灣的狀況而 言，高山區域的氣候資料多倚賴玉山氣象站的資訊，因此在山地的環境因子自相 關情形較為嚴重。
III. 多樣性評估分析
由物種預測分布結果顯示，遭受氣候變遷衝擊下，台灣山薰香屬的分布面積 的改變趨勢一致，皆有明顯面積縮減的現象，然而在遺傳多樣性的喪失中，假定 在相較物種演化的速率而言極短暫的百年內，物種的傳播能力既使適合長距離傳 播，但仍然侷限於小尺度的播遷，因此結合物種預測分布與遺傳多樣性的地理座 標資訊，進行遺傳多樣性喪失分析（Balint et al., 2011; Pfenninger et al., 2012）。相 較於分布面積改變的程度，遺傳多樣性喪失幅度較小，原因可能為棲地縮小，部 分稀有的基因單型並不能適應加劇的暖化而消失，或者具有未取樣的族群並無納 入分析所致。另外，由親緣網狀圖破裂情形顯示（圖二十），喪失的基因單型偏 向分支末端，靠近祖先型的基因單型保留較多，此有趣的現象可能與祖先型多分 布於分布核心區域有關，也暗示喪失的遺傳多樣性可能為稀有的基因單型。
藉由切分至物種以下的具遺傳訊息的操作單位，解析度提高至族群層級，對 於評估物種受到氣候變遷衝擊下的影響，明顯提高對於種內變異與氣候變遷之間
26
的關連，除了物種分布模型提供的族群整體的概況，更能評估細部各地理分區的 族群遺傳多樣性的影響程度，對於保育工作而言，提供更多資訊規劃各區域的保 育策略或評定優先順序。
然而以上問題可藉由提高分子標記的解析度並縮小研究區域，以釐清遺傳多 樣性消失的成因以及與微棲地的關係，達到有效率監控實際上的氣候變遷衝擊下，
物種真正受到的危害。
27
結論
本研究提供一個結合地理及遺傳資訊的實例，整合物種分布模型及族群遺傳 分析，評估高山植物受氣候變遷衝擊下的棲地變動及遺傳多樣性喪失，在台灣山 薰香屬的例子中，推論過去的演化歷史與間冰期避難所模型吻合，於本世紀的氣 候變遷中，其預測分布面積將受到氣候變遷的衝擊大幅縮小，並且最適分布海拔 提高 400 公尺以上，遺傳多樣性喪失的空間分布北、中、南部族群喪失比率最為 嚴重，加上東北往南方向遞減的脆弱度分布，顯示北部的遺傳多樣性核心區域受 到顯著的影響，與過去研究相符，本篇建議於北、中、南部各設置暖化監測樣區，
進一步驗證氣候變遷對於高山物種在演化適應及族群遷移，以及最重要的保種和 維護物種的遺傳多樣性。
28
參考文獻
吳明洲。2004。臺灣產碎雪草屬植物之系統生物研究。博士論文，臺灣大學，植 物科學研究所。
李彥希、顏宏旭。2011。全球暖化對國家公園環境變遷影響先驅計畫，第 1-132 頁。內政部營建署。
孫世鐸。2011。應用物種分布模式於氣候變遷下監測樣區之篩選。碩士論文，臺 灣大學，森林環境暨資源學系。
許晃雄、周佳、吳宜昭、盧孟明、陳正達、陳永明。2012。台灣氣候變遷的關鍵 議題。臺灣醫學 16: 459-470。
黃啟俊。2006。台灣產玉山杜鵑複合群之親緣地理學研究。碩士論文，成功大學，
生命科學系。
Balint, M., S. Domisch, C. H. M. Engelhardt, P. Haase, S. Lehrian, J. Sauer, K.
Theissinger, S. U. Pauls, and C. Nowak. 2011. Cryptic biodiversity loss linked to global climate change. Nature Climate Change 1: 313-318.
Bogaert, J., P. Van Hecke, and I. Impens. 1999. A reference value for the interior-to- edge ratio of isolated habitats. Acta Biotheoretica 47: 67-77.
Casalegno, S., G. Amatulli, A. Camia, A. Nelson, and A. Pekkarinen. 2010.
Vulnerability of Pinus cembra L. in the Alps and the Carpathian mountains under present and future climates. Forest Ecology and Management 259: 750-761.
Chou, C.-H., T.-J. Huang, Y.-P. Lee, C.-Y. Chen, T.-W. Hsu, and C.-H. Chen. 2011.
Diversity of the alpine vegetation in central Taiwan is affected by climate change based on a century of floristic inventories. Botanical Studies 52: 503-516.
Chung, J. D., T. P. Lin, Y. L. Chen, Y. P. Cheng, and S. Y. Hwang. 2007.
Phylogeographic study reveals the origin and evolutionary history of a
29
Rhododendron species complex in Taiwan. Molecular Phylogenetics and Evolution
42: 14-24.Chung, K.-F. 2006. Phylogenetics and phylogeography of the South Pacific alpine plant
Oreomyrrhis: Insights into global alpine biogeography. Doctoral Dissertation,
Washington University, Division of Biology and Biomedical Sciences, Program in Evolution, Ecology and Population Biology, St. Louis.Chung, K.-F. 2007. Inclusion of the south pacific alpine genus Oreomyrrhis (Apiaceae) in Chaerophyllum based on nuclear and chloroplast DNA sequences. Systematic
Botany 32: 671-681.
Chung, K.-F., C.-I. Peng, S. R. Downie, K. Spalik, and B. A. Schaal. 2005. Molecular systematics of the trans-Pacific alpine genus Oreomyrrhis (Apiaceae):
phylogenetic affinities and biogeographic implications. American Journal of
Botany 92: 2054-2071.
Clement, M., D. Posada, and K. A. Crandall. 2000. TCS: a computer program to estimate gene genealogies. Molecular Ecology 9: 1657-1659.
DeChaine, E. G. and A. P. Martin. 2004. Historic cycles of fragmentation and expansion in Parnassius smintheus (Papilionidae) inferred using mitochondrial DNA.
Evolution 58: 113-127.
Dormann, C. F., J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carré, J. R. G. Marquéz, B. Gruber, B. Lafourcade, P. J. Leitão, T. Münkemüller, C. McClean, P. E. Osborne, B. Reineking, B. Schröder, A. K. Skidmore, D. Zurell, and S. Lautenbach. 2013.
Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36: 027-046.
Elith, J., S. J. Phillips, T. Hastie, M. Dudik, Y. E. Chee, and C. J. Yates. 2011. A statistical explanation of MaxEnt for ecologists. Diversity and Distributions 17:
30
43-57.
Elith, J., C. H. Graham, R. P. Anderson, M. Dudik, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. Overton, A. T. Peterson, S. J. Phillips, K. Richardson, R. Scachetti-Pereira, R. E. Schapire, J. Soberon, S.
Williams, M. S. Wisz, and N. E. Zimmermann. 2006. Novel methods improve prediction of species' distributions from occurrence data. Ecography 29: 129-151.
ESRI. 2011. ArcGIS Desktop, version 10.
Excoffier, L., G. Laval, and S. Schneider. 2005. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary
Bioinformatics 1: 47-50.
Franklin, J. and J. A. Miller. 2009. Mapping Species Distributions: Spatial Inference and Prediction New York: Cambridge University Press.
Guisan, A., O. Broennimann, R. Engler, M. Vust, N. G. Yoccoz, A. Lehmann, and N. E.
Zimmermann. 2006. Using niche-based models to improve the sampling of rare species. Conservation Biology 20: 501-511.
Hijmans, R. J., S. E. Cameron, J. L. Parra, P. G. Jones, and A. Jarvis. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal
of Climatology 25: 1965-1978.
Hsieh, C.-F. 2002. Composition, endemism and phytogeographical affinities of the Taiwan flora. Taiwania 47: 298-310.
Hutchinson, G. E. 1957. Population studies - animal ecology and demography - concluding remarks. Cold Spring Harbor Symposia on Quantitative Biology 22:
415-427.
Kropf, M., J. W. Kadereit, and H. P. Comes. 2002. Late Quaternary distributional stasis
31
in the submediterranean mountain plant Anthyllis montana L. (Fabaceae) inferred from ITS sequences and amplified fragment length polymorphism markers.
Molecular Ecology 11: 447-463.
Lenoir, J., J. C. Gegout, P. A. Marquet, P. de Ruffray, and H. Brisse. 2008. A significant upward shift in plant species optimum elevation during the 20th century. Science 320: 1768-1771.
Liu, C. R., P. M. Berry, T. P. Dawson, and R. G. Pearson. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28: 385-393.
Marmion, M., M. Parviainen, M. Luoto, R. K. Heikkinen, and W. Thuiller. 2009.
Evaluation of consensus methods in predictive species distribution modelling.
Diversity and Distributions 15: 59-69.
Marske, K. A., R. A. B. Leschen, G. M. Barker, and T. R. Buckley. 2009.
Phylogeography and ecological niche modelling implicate coastal refugia and trans-alpine dispersal of a New Zealand fungus beetle. Molecular Ecology 18:
5126-5142.
Mathias, M. E. and L. Constance 1955. The genus Oreomyrrhis (Umbelliferae), a problem in south Pacific distribution. Berkeley: University of California
Publications in Botany 27 347-416.
Matsui, T., T. Yagihashi, T. Nakaya, H. Taoda, S. Yoshinaga, H. Daimaru, and N. Tanaka.
2004. Probability distributions, vulnerability and sensitivity in Fagus crenata forests following predicted climate changes in Japan. Journal of Vegetation Science 15: 605-614.
McClish, D. K. 1989. Analyzing a portion of the ROC curve. Medical Decision Making 9: 190-195.
Murray, M. G. and W. F. Thompson. 1980. Rapid isolation of high molecular weight
32
plant DNA. Nucleic Acids Research 8: 4321-4326.
Ohsawa, T. and Y. Ide. 2008. Global patterns of genetic variation in plant species along vertical and horizontal gradients on mountains. Global Ecology and Biogeography 17: 152-163.
Parry, M. L., O. F. Canziani, J. P. Palutikof, P. J. v. d. Linden, and C. E. Hanson [eds.].
2007. Climate Change 2007: Impacts, Adaptation and Vulnerability : Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, New York.
Pauls, S. U., C. Nowak, M. Bálint, and M. Pfenninger. 2013. The impact of global climate change on genetic diversity within populations and species. Molecular
Ecology 22: 925-946.
Pfenninger, M. and D. Posada. 2002. Phylogeographic history of the land snail
Candidula unifasciata (Helicellinae, Stylommatophora): fragmentation, corridor
migration, and secondary contact. Evolution 56: 1776-1788.Pfenninger, M., M. Balint, and S. Pauls. 2012. Methodological framework for projecting the potential loss of intraspecific genetic diversity due to global climate change.
BMC Evolutionary Biology 12: 224.
Phillips, S. J., R. P. Anderson, and R. E. Schapire. 2006. Maximum entropy modeling of species geographic distributions. Ecological Modelling 190: 231-259.
R Core Team. 2013. R: A language and environment for statistical computing. website:
http://www.R-project.org/.
Ramirez, J. and A. Jarvis. 2008. High Resolution Statistically Downscaled Future Climate Surfaces. Cali, Columbia: International Center for Tropical Agriculture
(CIAT); CGIAR Research Program on Climate Change, Agriculture and Food
Security (CCAFS).
33
Sala, O. E., F. S. Chapin, J. J. Armesto, E. Berlow, J. Bloomfield, R. Dirzo, E. Huber- Sanwald, L. F. Huenneke, R. B. Jackson, A. Kinzig, R. Leemans, D. M. Lodge, H.
A. Mooney, M. Oesterheld, N. L. Poff, M. T. Sykes, B. H. Walker, M. Walker, and D. H. Wall. 2000. Global biodiversity scenarios for the year 2100. Science 287:
1770-1774.
Shih, F.-L., S.-Y. Hwang, Y.-P. Cheng, P.-F. Lee, and T.-P. Lin. 2007. Uniform genetic diversity, low differentiation, and neutral evolution characterize contemporary refuge populations of Taiwan fir (Abies kawakamii, Pinaceae). American Journal
of Botany 94: 194-202.
Sinclair, S. J., M. D. White, and G. R. Newell. 2010. How useful are species distribution models for managing biodiversity under future climates? Ecology and Society 15: 8.
Soberon, J. and M. Nakamura. 2009. Niches and distributional areas: concepts, methods, and assumptions. Proceedings of the National Academy of Sciences 106: 19644- 19650.
Tamura, K., D. Peterson, N. Peterson, G. Stecher, M. Nei, and S. Kumar. 2011. MEGA5:
Molecular Evolutionary Genetics Analysis Using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Molecular Biology and
Evolution 28: 2731-2739.
Torres, R., J. Pablo Jayat, and S. Pacheco. 2013. Modelling potential impacts of climate change on the bioclimatic envelope and conservation of the Maned Wolf (Chrysocyon brachyurus). Mammalian Biology - Zeitschrift für Säugetierkunde 78:
41-49.
Tsukada, M. 1966. Late Pleistocene vegetation and climate in Taiwan (Formosa).
Proceedings of the National Academy of Sciences of the United States of America
55: 543-548.34
Waltari, E., R. J. Hijmans, A. T. Peterson, A. S. Nyari, S. L. Perkins, and R. P. Guralnick.
2007. Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions. PLoS ONE 2: e563.
Williams, S. E., L. P. Shoo, J. L. Isaac, A. A. Hoffmann, and G. Langham. 2008.
Towards an integrated framework for assessing the vulnerability of species to climate change. PLoS Biology 6: 2621-2626.
Wu, M.-J., S.-F. Huang, T.-C. Huang, P.-F. Lee, and T.-P. Lin. 2005. Evolution of the
Euphrasia transmorrisonensis complex (Orobanchaceae) in alpine areas of Taiwan.
Journal of Biogeography 32: 1921-1929.
35
圖一、 氣候變遷衝擊下的物種脆弱度評估架構（Williams et al., 2008）
