地理與生態因子在島嶼上黃芩屬物種的族群分化與快速種化上的效應

全文

(1)國立臺灣師範大學理學院生命科學系博士論文 Department of Life Science College of Science. National Taiwan Normal University Doctoral Dissertation. 地理與生態因子在島嶼上黃芩屬物種的族群分化與快速種化上的效應 The geographic and ecological effects on the population differentiation and rapid divergence of island Scutellaria (Lamiaceae). 黃秉宏 Huang, Bing-Hong 指導教授：廖培鈞博士 Advisor: Liao, Pei-Chun, Ph.D. 中華民國 109 年 6 月 June 2020.

(2) 致謝博士班的時光雖不算長，卻也深深的影響了我非常多，不管是知識拓展還是辯證思考等等。我的研究都是由我的指導教授廖培鈞教授主持的科技部計畫所支持的(MOST 102-2621-B-003-005-MY3)，亦師亦友的廖老師也在平時予我啟發亦多，這份論文沒有廖老師的支持，是完全無法完成的，裡面的論文若有任何帶來的成就或喜悅，廖老師真的是功不可沒。此外，多項分析方法例如 ordination analysis 或 regression based analysis 在黃老師的指導下才能運用在論文中，其中在口試前後，黃老師也不斷提醒我在使用這些方法的迷思，獲益良多。而江友中老師是這個類群的啟蒙老師，從我在當助理的時候就將許多材料的訊息給我，而我也有一部分的樣本蒙江老師提供，從 purposal defense 至口試期間對於這類群的知識跟討論都不吝提供，也啟發非常多關於論文及材料相關的切入面相。而李勇毅老師一路互相討論學習，也常常談及一些喜愛的蘭科科學議題，多次讓我的議題也有新的思考方向，有時候我都會想，我是否能鍾情的把一類型的植物做到極致，並能像李老師與他們的團隊一樣，了解到如此透徹呢？而指導老師廖老師的介紹之下，趙怡姍老師與我的合作緊鑼密鼓的進行，期間相互討論的過程中，我也在此就彼此想法激盪出對我論文討論及未來進展有關的方向，另外，趙老師在蕨類雜交議題之精，相比之下黃芩的雜交真是太單純了，有時候在合作過程我也接觸不少在對雜交的科學議題想像，我的界線還是太小了。研究過程中，協助最大的，莫過於一起參與進來的怡雯、蕭錦隆與桓翊，怡雯與我一起經營實驗室黃芩的維持，這些材料都藉她的巧手才能活下來，換我獨自用我的 black finger 大概樣本數會瞬間少一半吧？蕭錦隆先生協助台北黃芩研究的鑑定與地點指出幫助特別多，而桓翊在大學期間協助補足黃芩的 genotyping 才能讓台北黃芩的故事更加完整。實驗室的睿澤(樂高)協助種子 SEM 拍攝，以及與敏心升成博士班與準博士班，雖說是後進，但是以很快的速度共同接手我放手的行政事務，讓我能專心在論文上，也提供很多討論，這是我望塵莫及的速度，接下來進步的棒子就交給你們了，期待早日能看你們一起升成博士候選人甚至畢業到各自經營相互合作的日子。實驗室經過的成員們，已畢業的紹瑋、鈞凱、怡婷在各項實驗互相精進的日子，現在看你們各有出路，非常為你們開心。致維學長在期間的學術各項討論，也受益良多，雖然跟學長現在在不同單位，還希望能有再交流的機會。然後在撰寫這本論文的時間，時任博後的建棣一直對我的各種分類議題有非常多的協助，甚至延伸台北黃芩的研究建立的了一篇正式分類處理，將那份研究最後一塊拼圖補齊。遠在北京的高健也在這些議題上交換意見許多，倍感溫馨，欣蓓在教學後也在我實驗上有相當多的協助及合作，一直是教學相長啊。明威、祖恩、沛煒跟豐愷也是實驗室的後起之秀，平常討論真的少不了。世穎老師實驗室的翊韶跟偉銘在最後口試及文章的協助指教也讓文章流暢些，在這麼趕的狀況下，真的非常謝謝上面各位的大力幫忙。最後，這篇論文並非我的終點，我依然有許多以台灣產植物為材料，探討台灣植物物種、棲地及遺傳多樣性的議題想繼續探討，對於看著這本論文的讀者，如果願意給我任何指教，我會衷心的感謝。. I.

(3) 目次致謝................................................................................................................................ I 摘要................................................................................................................................ 1 Abstract .......................................................................................................................... 3 第一章、概論................................................................................................................ 5 1.1 物種分化的力量............................................................................................. 5 1.2 高特有種比例的島嶼例子............................................................................. 6 1.3 研究材料介紹................................................................................................. 6 圖 1-1 台灣產黃芩物種地理分布簡圖 ....................................................... 8 1.4 本論文研究議題............................................................................................. 9 表 1-1 台灣之黃芩屬植物分布區域。 ....................................................... 9 第二章台灣產近緣黃芩屬植物的分歧及種間基因交流測定................................ 10 2.1 前言............................................................................................................... 10 2.2 材料與方法................................................................................................... 12 2.2.1 材料收集............................................................................................. 12 表 2-1、採樣地訊息及各物種及族群代號 ............................................... 12 2.2.2 基因多樣性整理................................................................................ 13 2.2.3 遺傳組成及基因交流檢測................................................................ 13 2.3 結果............................................................................................................... 14 2.3.1 遺傳多樣性........................................................................................ 14 表 2-2 五種黃芩及其族群的遺傳多樣性指數計算平均值及其 95%信賴區間。.......................................................................................................... 15 2.3.2 物種遺傳結構分析............................................................................ 15 圖 2-1 以 R 套件 vegan 內的 find.cluster 指令算出各個以五種黃芩的 SSR 遺傳資料計算出來 cluster 數對應的 BIC 值，BIC 變化轉折處約莫在 10 左右。.......................................................................................................... 16 圖 2-2 五種黃芩遺傳組成及群聚分析。(A) 利用 DAPC(上)及 STRUCTURE(中、下)進行的 individual assignment test 的結果，族群及物種代號見表 5-1。(B) 利用 DAPC 的結果，以 BIC 值的變化約略在 10 個 cluster 為最好的分群數，以此繪製的散佈圖。 ................................. 17 圖 2-3 檢視 STRUCUTURE K=10 的結果發現該群聚的 log likelihood 約略分成兩群，重新以 log likelihood 較高的獨立運算重新取得共識，相較於圖 2-2A, K=10 的結果，未發現明顯的種間遺傳混雜，但有更明顯的族群間遺傳分化情形。.............................................................................. 17 2.3.3 基因交流計算.................................................................................... 17 圖 2-4 以 IMa3 計算出各個分支間基因交流的程度。樹型及分化時間由於無確切參考的微衛星體突變率，故參照 Chiang et al. (2012a)之樹型及分化時間(kilo years ago, kya)。箭號及其大小代表基因交流的相對強度及方向，並以紅色代表 indica group 內的基因交流，而藍色是至少有一端為非 indica group 成員的基因交流。 ................................................... 18 2.4 討論................................................................................................................ 18 II.

(4) 2.4.1 更新世末期的台灣產黃芩具有旺盛的種間基因交流..................... 18 2.4.2 族群間具有明顯結構......................................................................... 19 第三章島嶼上的空間環境異質性趨使島嶼物種持續分化以致生態種化............ 22 3.1 前言................................................................................................................ 22 圖 3-1 布烈氏黃芩(S. playfairi，淺藍色)及田代氏黃芩(S. tashiroi，淺黃色)於台灣島的分佈圖。兩者實際的標本紀錄詳見圖 3-5。圖中可見兩者於中央山脈(底圖深色處為海拔高度，愈深者海拔愈高)東南處有明顯的重疊。本圖利用 R 套件 maptools (Lewin-Koh et al., 2011)繪製而成。 24 3.2 材料與方法.................................................................................................... 25 3.2.1 採樣地點與遺傳分析......................................................................... 25 表 3-1 田代氏黃芩與布烈氏黃芩於各個採樣點的物候特性 ................. 25 3.2.2 遺傳分析............................................................................................. 25 3.2.3 利用 ABC 評估物種種化模型 .......................................................... 26 圖 3-2 藉 ABC 演算評估田代氏黃芩與布烈氏黃芩的演化情境：(A) 持續隔離模型(CI)；(B) 祖先基因交流模型(AM)；(C) 二次接觸模型(SC)；以及(D) 持續基因交流模型(CM)。其中 AM 有最高的 posterior probability，各個以 AM 情境底下計算出來的參數及其 marginal density distribution 也羅列如下：(E) 基因交流 (F) 有效族群大小(effective population sizes) (G) 分歧與發生基因交流的時間，及 (H) 中性基因座的平均突變率。MD, marginal density; RD, relative density. ................... 27 3.2.4 以降維的環境因子計算物種間的生態區位差異............................. 27 3.2.5 利用生態區位模擬來預測物種分佈趨勢......................................... 28 3.2.6 利用冗餘分析(redundancy analysis)評估環境變量對遺傳變異的影響.................................................................................................................. 29 3.3 結果................................................................................................................ 29 3.3.1 遺傳多樣性與族群結構..................................................................... 29 圖 3-3 以 BayeScan，依照 Bayes Factor 大於 10 作為標準。找出 3 個 positive-outlier 基因座 (aus9-2, aus9-3, aus9-5) 以及一個 negative-outlier 基因座 (ST5-179) ....................................................................................... 30 表 3-2 利用 19 組微衛星基因座計算出之各個族群的遺傳多樣性指數31 表 3-3 不同基因座類型進行 AMOVA(analysis of molecular variance)結果整理.............................................................................................................. 31 圖 3-4 以 STRUCTURE 軟體進行 BCA 分析的結果(A) ΔK 以及 lnP(K) 在不同的歸群數(K)的變化圖。紅色線為其十次獨立模擬計算出來 lnP(K)，也就是各個 K 值計算出的 posterior probability 及其標準差，而黑色線則呈現 K 值變化 ΔK 變化的趨勢。結果顯示 ΔK 在 K = 2 來到最大值。(B) 當 K=2 或 3 時，每個個體計算出來的遺傳群聚結果(Q-matrix)，這裡同時呈現中性基因座及適應性基因座的結果。.......................................... 32 表 3-4 微衛星體各基因座的遺傳多樣性 ................................................. 33 圖 3-5 藉 DAPC 對樣本進行遺傳結構分析。分別以 density plot(A)以及 bar plot(B)呈現物種的階層各個樣本的遺傳組成，僅有極少量的遺傳混雜被發現。另也以分 scatter plot(C)以及 bar plot(D)以族群的階層分析各個樣本的遺傳組成。可發現除田代氏黃芩的霧鹿族群外，其餘族群或多或少都有些遺傳混雜及結構。.................................................................. 34 III.

(5) 3.3.2 種化模型............................................................................................. 34 表 3-5 利用 ABC 計算出 AM model 的各項參數.................................... 35 3.3.3 布烈氏黃芩與田代氏黃芩的棲位分化............................................. 35 圖 3-6 以五個環境因子用 Principal component analysis (PCA) 繪製而成的二維 (A) 及三維 (B)的散佈圖。PCA 的結果顯示兩物種的棲位有很大的重疊。.................................................................................................. 36 圖 3-7 以 Kruskal–Wallis 檢測環境因子轉換而成的前三 PC 軸是否在田代氏黃芩及布烈氏黃芩間有顯著分化...................................................... 36 表 3-6 利用 Generalized linear model (GLM) 檢測五項參數是否能顯著解釋環境參數轉換而成的 PC1 .................................................................. 37 表 3-7 利用多變數邏輯迴歸檢測環境變量與物種分化之關係 ............. 38 3.3.4 生態棲位模擬..................................................................................... 38 3.3.5 環境變量能否解釋兩種黃芩的遺傳變異......................................... 40 表 3-8 The 利用 generalized linear model (GLM)來檢視環境因子與 dbRDA 前兩軸之相關性 ............................................................................ 40 圖 3-8 以 dbRDA 繪製而成的 ordisurf 散布圖，並檢視四個環境因子在 ordisurf plot 上面的方向。 ......................................................................... 41 3.4 討論................................................................................................................ 41 3.4.1 非異域種化的遺傳證據..................................................................... 41 3.4.2 生態因子趨使物種持續維持生殖隔離............................................. 43 表 3-9 distance-based redundancy analysis (dbRDA)的結果統整。我們以遺傳資料轉換成的 PC1 及 PC2 作為遺傳組成(SSR)，而 dbRDA 計算公式為 SSR ~ alt + AET + bio3 + bio14 ........................................................ 43 表 3-10 以 dbRDA 檢測各項限制因子能否解釋遺傳多樣性 ................. 44 3.5 結論................................................................................................................ 45 第四章區域性環境異質性形塑兩種未分化完全的島嶼物種的遺傳組成............ 46 4.1 前言................................................................................................................ 46 圖 4-1 藉 19 個 bioclim 因子將(A) 向天盞與(B)台北黃芩的標本記錄點位(黑色十字標記)進行物種空間模擬分佈預測。 ................................... 48 圖 4-2 台北黃芩與向天盞的各種演化情境。並以 ABC 進行情境比較及參數演算。(A)兩物種分化後歷經瓶頸效應後再次產生近期基因交流(後稱二次接觸情境)；(B)台北黃芩為向天盞近期拓殖的族群，並持續與向天盞有基因交流(後稱創始者情境) (C)台北黃芩為多次由向天盞拓殖出來的族群；(D)台北黃芩為單次由向天盞拓殖出來的族群。遷移率在此由 ABC 計算，虛線外框代表物種分歧的情形。 ................................... 50 4.2 材料與方法.................................................................................................... 50 4.2.1 分佈現況與採集................................................................................. 50 圖 4-3 本研究採樣的各個族群小堅果紋路。由 SEM 的圖顯見個體見的種子大小有相當大的差異，但種間並無明顯差異。最右下角的圖是以 TC 族群為例拍攝未成熟的小堅果照片，箭頭所指處即為傘狀紋路所在。代號為 MK 的為台北黃芩的模式產地的樣本。 ..................................... 51 4.2.2 定序與基因型讀取............................................................................. 51 4.2.3 中性檢測與遺傳多樣性..................................................................... 52 表 4-1 採樣族群資訊及當地的 bioclimatic variable. ................................ 53 IV.

(6) 4.2.4 利用貝氏群聚分析來檢視族群的遺傳組成..................................... 54 4.2.5 利用 Discriminant analysis of principal components 法區分物種及族群.................................................................................................................. 54 4.2.6 利用 approximate Bayesian computation 評估種化情境 .................. 54 4.2.7 藉生態棲位模擬預測物種潛在分佈地點......................................... 55 4.2.8 利用 Mantel test 檢測 isolation-by-distance 及 isolation-by-environment ............................................................................. 56 4.2.9 檢測氣候對族群遺傳組成的影響..................................................... 56 4.3 結果................................................................................................................ 57 4.3.1 成熟種子的小堅果外衣紋路形態差異............................................. 57 4.3.2 中性檢測及遺傳多樣性..................................................................... 57 圖 4-4 以 FST based analysis 對 11 組微衛星基因座進行中性檢測(a) Fdist 法 (b) BayeScan 法。兩種方法都未找出 FST 值過高的(positive outlier)或過低者(negative outlier)，顯示這些基因座都未違背選汰中性。 .......... 58 表 4-2 各採樣族群的遺傳多樣性參數列表 ............................................. 59 圖 4-5 向天盞(XH, YH, WL,CL, TC 及 YL 族群)和台北黃芩(DA, MK 及 EG)採樣點對應的葉綠體基因座單套型，以及其重建的 minimum spanning network。 ..................................................................................... 60 表 4-3 利用 Analysis of molecular variance (AMOVA) 分析兩種黃芩的遺傳變異.......................................................................................................... 60 4.3.3 向天盞與台北黃芩的族群結果......................................................... 60 圖 4-6 利用 STRUCTURE 進行 BCA 分析之結果(A) 以 Evanno et al. (2005)法利用 ΔK 評估最歸群數為 3，其次為 2，顯示這兩物種大致上可分為(B)兩到(C)三群。 ............................................................................... 61 圖 4-7 由 DAPC 的結果顯示台北黃芩與向天盞的族群結構。(A)利用單一 discriminant function 區別向天盞(bar，藍)及台北黃芩(tpe，紅)。 (B)DAPC 以兩個 discriminant function 做的散佈圖，藍色系者為向天盞，而紅色系者為台北黃芩。.......................................................................... 62 4.3.4 台北黃芩的多次起源與近期族群分化可用以解釋為何兩種間無明顯遺傳分化.................................................................................................. 62 圖 4-8 利用 Kruskal–Wallis test 來檢視微衛星基因座轉換過的 PC/DA 軸是否能顯現物種階層或族群階層的遺傳分化。...................................... 63 圖 4-9 ABC 運算結果彙整 (A-D)呈現圖 4-2B 內各個參數的 density plot： (A)估算的祖先有效族群量 NA 及 NB；(B)現今族群的 NB 及 NT 的有效族群量；(C)各個遷移率 (D)向天盞的溯祖時間(T1)及台北黃芩的拓殖時間 (T2)。(E-H)呈現圖 4-2C 內各個參數的 density plot：(E)向天盞的有效族群量(NB)及三個台北黃芩各自的有效族群量(NT1-NT3)；(F)三次台北黃芩拓殖的時間；(G)遷移率；(H)圖 4-2B 及圖 4-2C 各自估算出來的突變率 (μB 及 μC)。 .................................................................................................. 64 4.3.5 台北黃芩與向天盞的 Grinnellian niche 明顯交疊 ........................... 65 表 4-4 多項式邏輯迴歸分析檢測各環境因子對物種出現能否解釋，並利用 type-II ANOVA 檢測各別因子的顯著性。 .......................................... 65 圖 4-10 將向天盞及台北黃芩族群的各氣候因子調閱出來進行 PCA 分析，(a)僅用溫度類的 bioclimatic variable(bio2, bio8,及 bio9) (b) 僅用降 V.

(7) 雨類的 bioclimatic variable(bio13, bio18,及 bio19) 。氣候因子參照 WorldClim (http://www.worldclim.org/bioclim)資料 ................................. 66 4.3.6 藉生態棲位模擬預測物種的分佈..................................................... 66 表 4-5 統整 Mantel test 與 multiple matrix regression with randomization (MMRR)藉以了解族群遺傳組成(gen), 地理距離(geo)及環境差異(env) 之間的關聯性.............................................................................................. 67 4.3.7Isolation-by-distance 及 isolation-by-environment 的模型比較 ........ 67 4.3.8 以環境因子解釋族群遺傳結構......................................................... 68 4.4 討論................................................................................................................ 68 4.4.1 向天盞與台北黃芩的形態特徵難以區別......................................... 68 4.4.2 遺傳及棲地分化不明顯暗示種化不完全......................................... 69 表 4-6 dbRDA 分析結果統整，顯示氣候因子(限制因子)顯著解釋遺傳組成的變異...................................................................................................... 70 4.4.3 氣候因子並不與物種區別有關......................................................... 71 4.4.4 氣候異質性致使區域性適應............................................................. 71 4.4.5 夏末初秋的降雨是主要解釋遺傳分佈的環境因子......................... 72 圖 4-11 利用 partial dbRDA 針對六個環境因子進行散佈圖及 ordisurf plot ............................................................................................................... 73 圖 4-12 向天盞與台北黃芩的各個記錄點統出而成的(a)每季(三個月)的平均溫度及(b)每雨的降雨量。 ................................................................. 74 4.5 結論................................................................................................................ 75 第五章結語................................................................................................................ 76 附表.............................................................................................................................. 78 附表 1 各個 bioclimatic variable 代號對照表 ........................................... 78 附表 2 向天盞(S. barbata)與台北黃芩(S. taipeiensis)的葉形差異 .......... 79 附表 3 作者於博士班期間協助實驗室進行的研究發表......................... 80 參考文獻...................................................................................................................... 83. VI.

(8) 摘要. 研究背景：島嶼物種的快速種化可以由隨機的力量、遷徙拓殖、基因交流以及天擇等力量驅動。然而島嶼面積不大，頻繁的物種接觸導致的基因交流理應減少物種或族群分化，現今我們常以生態力量趨動的生態種化解釋之，因為生態因子如環境氣候異質性等使族群間有區域性適應阻擋部分基因組的基因交流，區域性適應影響的基因組週邊區域因連鎖而產生搭便車效應加強了該區段的基因交流阻隔，而漸漸的這些阻隔區域增加以到一定程度並完全種化，這些力量可以共同作用以加快島嶼物種的分化速度。研究具高特有種比例的物種有助了解這些力量在祖先物種拓殖入島嶼後加快分化速度的相對重要性。台灣產的黃芩屬特有種比例比起其他東亞及東南亞的鄰近島嶼或陸塊都高上許多，全部共 8 種，即有 6 種為特有種。過去的研究指出台灣產黃芩屬的多樣性來自至少三次的拓殖事件以及島嶼在地的快速種化。而台灣產的黃芩利用分子定年計算出的最早分歧時間不過 50 萬年內，其中被命名作 indica group 的印度黃芩、布烈氏黃芩、田代氏黃芩及南台灣黃芩形成單一分支且分歧時間短於 20 萬年內，被認為可能是其中一次拓殖事件後的在地快速種化。我們想以台灣產黃芩屬成員作材料切入生態種化如何與台灣高特有性之關聯性。方法我們使用了微衛星體基因型檢測作為快速變異的遺傳資料進行分析。並輔以溯祖理論分析法及演化模型選擇找出最適合解釋台灣產黃芩是否一如生態種化之預測分化早期有基因交流。並以環境因子與遺傳分化進行關聯性分析，了解遺傳變異或族群分化是否與環境因子有關。結果我們優先了解這群物種來檢測 indica group 是否主要來自單次的快速種化，以利後續設計實驗了解台灣產黃芩的快速種化之因。然而，依據溯祖理論 IMa3 分析，我們發現 indica group 在分歧時伴隨旺盛的基因交流，然而在使用貝氏歸群分析時卻顯示這些近期分歧且鄰域物佈的物種甚少種間遺傳混雜，顯示物種明顯分化， 1.

(9) 考量過去旺盛的基因交流及台灣島面積不大，物種易接觸，這麼短的分歧時間實難分化出這麼多特有種，我們藉分歧時有種間基因交流的情形，再依生態種化的概念，推論生態因子可能參與在這些快速種化的黃芩的物種維持與分化。我們挑選兩組共域/鄰域分佈的特有種，分別是田代氏與布烈氏黃芩組，以及向天盞與台北黃芩組來問答這個議題。我們利用以 approximate Bayesian computation 的方法進行溯祖模型選擇，發現田代氏黃芩與布烈氏黃芩有強烈的祖先性基因交流，然而現今已有明顯的分化及基因交流的阻斷。兩個物種的分化都可以由在地環境異質性及多環境因子組合的棲位分化來解釋。而向天盞及台北黃芩則展現不同的結果，不論是溯祖模擬模型選擇、群聚分析及單套形的樹形/網狀分析都無法清楚的區分特有的台北黃芩及廣布於東亞的向天盞。而可供區別的形態特徵如小堅果外衣紋路及葉型，在觀測後也證實僅是採集植物個體發育時期不同所致，並非可供物種區分的特徵。組合台北黃芩及向天盞的族群後發現族群間有明顯的遺傳結構，這些遺傳結構以廣義線性模型及座標軸分析等方法發現夏秋降雨量在族群的差異與族群遺傳分化有關。結論因此，無論是向天盞/台北黃芩組，或是布烈氏黃芩/田代氏黃芩組都發現明顯的族群尺度分化，這些分化都可以由氣候異質性及區域性適應解釋。換言之，氣候異質性及區域性適應均能促進台灣產黃芩屬成員讓族群分化，更可能加速種化以讓多樣性速率提升。. 關鍵字：黃芩、快速種化、棲位模擬、ABC 溯祖理論分析. 2.

(10) Abstract. Background: Species diversification within an island can be invoked by stochastic process (eg. drift or geographic isolation), migration (or gene flow) and ecological factors (eg. divergent selection lead population differentiation). However, gene flow is often considered to homogenize the genetic component between populations or species as area of island may be so small that species may contact to each other. The ecological factors such as environment heterogeneity may lead local adaptation which can block gene flow in parts of the genome, and the hitchhiking effect may expand the isolation area of genome until the speciation complete. These processes are not mutually exclusive and may accelerate diversification rate of colonized species. Genus with high endemism in an island is therefore suitable for accessing the relative importance of both neutral and adaptive process on diversification of island species. The genus Scutellaria L. comprised eight species with six of them are endemic species, which is suitable for ecological speciation study. The endemism is almost the highest in East Asia. Previous studies revealed very recent divergence of these species (within 500 kilo years ago, kya) and at least three times colonization with subsequently radiation may invoke the high endemism of Taiwanese Scutellaria. Amongst them, S. indica, S. tashiroi, S. playfairii and S. austrotaiwanensis (hereafter: indica group) are the most monophyletic group in which relatively short divergent time (less than 200 kya). We aimed to ask whether the ecological species can be associated with rapid species diversification in Taiwanense Scutellaria. Methods We adopt microsatellite as our molecular marker due to its high mutation rate and variation. We also used coalescence analyses and evolutionary model choice to evaluate the magnitude of historical gene flow, so that we can check whether species divergence was accomplished with gene flow as expectation of ecological speciation. Results According to results of coalescence simulation using IMa3, highly historical gene flow can be found. among members of indica group during divergence, but Bayesian. clustering analysis revealed very few interspecific admixture, implying clear species 3.

(11) divergence. Homogenizing genetic components between species can be often accived by historical gene flow. Therefore it is hard to believe that geographic factors along can lead to complete speciation in these rapid divergent Scutellaria. The historical gene flow often correlated to ecological speciation. Hence the ecological factors may also involve in these rapid divergent Scutelliara. We focus whether ecological factors can be associated with the divergence of the endemic Scutellaria in Taiwan. We adopt two sympatric or parapatric and morphological similar pair species as case, S. tashiroi/S. playfairii and S. barbata/S. taipeiensis, respectively, for understanding whether differentiation can be associated with neutral geographic distance or climatic driven local adaptation. We found that the clear divergence, local adaptation and niche differentiation can be found in former group even with consideration of geographic factors. Though gene flow can be detected with approximate Bayesian computation coalescence simulation analysis but differentiation in multidimensional hypervolume niches blocks the gene flow nowadays, indicating heterogeneity in environment drive divergence of S. tashiroi and S. playfairii. On the other hand, coalescence simulation, clustering analysis and tree/network analysis do not found the divergence between widespread species, S. barbata, and endemic S. taipeiensis. The distinct characters such as nutlet coat and leaf shape merely reflect different developmental stage or phenotypic plasticity. However, population level differentiation can be found to be associated with climatic factors such as Summer-Fall precipitation. Conclusions These two paired studies revealed that local adaptation driven by environmental factors may lead population differentiation and may subsequently lead to rapid speciation.. Keywords: Scutellaria, rapid speciation, niche modeling, ABC coalescence simulation. 4.

(12) 第一章、概論 1.1 物種分化的力量促進物種分化的力量相當多，例如異域分佈、天擇、環境因子及區域性適應等，然而基因交流在此扮演的角色則相當多元。Wright (1931)就嘗試討論過 gene swamping 一類的議題-當基因交流夠強，而當地特有的環境壓力不夠大時，基因交流可以將不適應的等位基因(maladaptive allele)帶入當地，相反的，如果當地特有的環境壓力夠大時，族群依然可以藉由這種天擇力道完成環境互斥天擇 (environmental antagonistic selection)保持住族群內獨有的適應性等位基因，並藉以完成區域性適應(local adaptation)並維持族群分化(Smith, 1966; Endler, 1977)，甚至進一步促成種化。這樣的論點雖然發展久遠，但仍然需要一定程度的地理隔離或是鄰域分佈等地理障礙在。一直到 Wu (2001)提出藉分歧性天擇與基因組結構(genomic architecture)共同作用，受分歧性天擇的位點及周邊連鎖的區域抵抗基因組部份區域不會有因基因交流帶來的外來等位基因，而搭便車效應(genome hitchhiking)逐漸擴大這樣區域，基因組的其他基因座也因應生成新的受分歧性天擇基因座減少物種雜交產生的基因交流，受天擇影響的基因座與週邊區域連鎖，產生搭便車效應進而持續擴大隔離的區域，最終生殖隔離產生，物種便可在有基因交流的情況分化而不需有地理隔離的時期。直到基因定序普及化，以及溯祖理論大量應用在種化議題上，在 isolation-with-migration 的概念就此形成(Wakeley & Hey, 1998; Hey, 2006; Hey, 2010; Pinho & Hey, 2010) 並開始被大量應用在各式各樣的物種中，。並開始應用在種化、天擇與基因交流上(cf. Niemiller et al., 2008)。然而，實際上現生共域的物種並一定是共域種化所致，例如過去分歧歷史曾有異域分佈或鄰域分佈的時期(cf. Albert B. Phillimore et al., 2008)，又或者在一些基因交流受限但分佈廣的物種呈現 isolation-by-distance(IBD)來解釋族群的分歧 (Hoelzer et al., 2008)。除此之外，像是連續性的環境變異(例如溫度)也可能限制了不適應(maladaptative)的遷入者(immigrantion)進而限制了遷移方向(bias dispersal)造成族群分化，這種稱之為 isolation-by-environment(IBE)，一樣能促進種化(Wang & Bradburd, 2014)。然而，相鄰的區域常有相似的環境條件，致使環境因子常與地理因子呈現一定程度的相關性，導致要區別兩者並不容易。因此，利用現今常利用 regression based analysis 或 ordination analysis 來區別地理和環境 5.

(13) 因子對於族群間或近緣種間分歧的相對重要性來討論物種間的分歧(Lee & Mitchell-Olds, 2011; Sexton et al., 2014; Huang et al., 2016c)。 1.2 高特有種比例的島嶼例子地理距離及隔離常是促使島嶼物種異域種化的關鍵，若島嶼上有破碎崎嶇的地形，也可能賦予島嶼有多樣氣候異質性，促成分歧性天擇(divergent selection)、區域性適應(local adaptation)或 IBE。而兩者也可能同時促進島嶼的族群分化或加速物種的種化(Kubota et al., 2014; Kubota et al., 2015)。而且相較大陸而言，小尺度的地理距離的島嶼因板塊擠壓常產生破碎地形或是海拔高低落差，破碎與海拔落差大的地形更能加強島上地理隔離及氣候區域性適應的強度，甚至可能加速種化(cf. Wanderley et al., 2018)。其中，找出可能解釋物種分化的環境因子，可以進一步了解生態因子如何趨動族群分化甚至到種化-生態種化(ecological speciation)(Nosil et al., 2009a; Nosil et al., 2009b; Nosil & Feder, 2012)。除此之外，利用環境因子進行物種生態棲位模擬(ecological niche modeling)可以找出與物種出現相關的環境因子以外及種間棲位分化外(Gorel et al., 2019; Hu et al., 2019)，也可以了解種以內不同地區的族群分化由哪些氣候因子主導(Gorel et al., 2019)。要了解這些議題，大量的遺傳資料、合適的島嶼物種以及齊全的環境資料是必要的。台灣是相當適合研究島嶼生物地理、生態種化、IBD/IBE 等相關議題的好地點。台灣位處熱帶與亞熱帶交界，為大陸型島嶼，在冰河期海平面下降後，大陸棚形成陸橋，讓台灣周邊陸塊的物種可以遷移至台灣，而冰河期結束後海平面上升導致陸橋消失，到達台灣的族群旋即與其他族群形成隔離，使得台灣物種來源多元(Hsieh, 2002)，加以台灣海拔落差非常大，落差可將近 4,000 公尺，小小一片島嶼溫度落差極大，相當於可以一路由熱帶雨林看到類似寒帶的氣候，因此有著足夠的氣候異質性，而東亞島弧中又以台灣面積最大，加以中央山脈鎮於島中，使得東西部台灣有一定程度的地理阻隔，這些指標都有機會加速島內物種跟族群的分化，這也可能反應台灣特有種比例高的原因(Hsieh, 2002; Hsu & Wolf, 2009; Jump et al., 2012)。 1.3 研究材料介紹 6.

(14) 其中，黃芩屬是吻合這條件的材料之一。黃芩(Genus Scutellaria)為唇形科黃芩亞科(Scutellarioideae)植物，全世界約有 350 種，亞洲大陸約百餘種(Li and Hedge 1994)，而台灣產共有八種，其中六種為特有種(台北黃芩[S. taipeiensis]、台灣黃芩[S. taiwanensis]、田代氏黃芩[S. tashiroi]、布烈氏黃芩[S.playfairi]、南台灣黃芩 [S. austrotaiwanensis]、長葉黃芩[Hsieh (2013)])及兩種非特有種(印度黃芩[S. indica]、向天盞[S. barbata])。台灣產特有種黃芩屬成員除了長葉黃芩僅見於南投人倫林道以外，其餘彼此間呈現共域(sympatric)或鄰域(parapatric)分布(表 1-1、圖 1-1)(Hsieh & Huang, 1995, 1997; Hsieh, 2013)，台灣產黃芩均屬黃芩亞屬 (Scutellarioideae)，過去研究了解黃芩屬底下的分類、染色體數與其地理分佈三者之間有相當明顯的趨勢，例如 2n=26 均屬黃芩亞屬，而 n=13(2n=26)在東亞及東南亞也最為常見，並被認為是東亞這帶最基礎的染色體數(basic chromosome number)(Ranjbar & Mahmoudi, 2013)，而台灣產的黃芩不僅形態相似(表 1-1)，染色體數也均為 2n=26，除暗示這些黃芩親緣關係之近，也排除這類黃芩屬成員以基因組尺度的重組造成種化的可能性(eg. Polyploidy)，其中，布烈氏黃芩、南台灣黃芩及田代氏黃芩關係又更為密切，這三者在過去曾被處理為呂宋黃芩(S. luzonica)底下的變種，外形也極為相似，僅可由田代氏黃芩的花冠筒基部無垂直折曲而另外兩者有，以及葉柄跟葉片刻痕比例區別布烈氏黃芩及南台灣黃芩 (Yamazaki, 1992; Hsieh & Huang, 1995, 1997)。另外一組則為向天盞與台北黃芩，區別特徵僅有葉長寬比以及小堅果外衣紋路可供鑑定參照(Huang et al., 2003)。兩者的整合葉綠體與核基因的分析也極少變異，定年也短於 2 萬年內(Chiang et al., 2012a)，表示他們分歧時間約莫在最末盛冰期時分歧(Last glacial maximum, LGM) 至今，更需謹慎參照整合以了解各項因子貢獻在其快速種化的比例如何。根據我們過去對黃芩植物的採樣及觀察發現部份分布重疊之種類，在共域棲地可找到中間型或連續性變異之小型植株(如台東新武呂溪一帶布烈氏黃芩與田代氏黃芩)，台灣產的黃芩屬雖未有實質報導授粉者是誰，然而已有大量的同屬其他成員被報導為小型蜂或蠅類(Zhang et al., 2011; Wiesenborn, 2019)，與野外觀察一致，傳粉者分化的可能性因此不高。暗示生殖隔離不完全或有基因漸滲 (introgression)之可能性，共域或鄰域分布的近緣種中可能因為接觸而產生種間基因交流，而基因交流容易帶入不適應當地的等位基因，間接削弱分歧性天擇的力 7.

(15) 量，或是藉由帶入其他物種的適應性等位基因加速適應性演化，進而影響種化的速度(Servedio et al., 2013)，了解這些力量彼此的相關性有助以黃芩作為起點，了解台灣植物的高特有性與環境之間的關係。除此之外，由於黃芩屬成員散佈能力有限，種子多以果實開裂時的小堅果彈出的重力傳播，僅偶發性的洪水可能協助長距離傳播(Williams, 1992)。而傳粉者多為小型蜂類，飛行能力有限且不具備專一性(Zhang et al., 2011)。再加以黃芩屬成員多為自交親合且許多成員可以走莖或地下部作無性繁殖，族群內遺傳多樣性因此大多偏低，最後，黃芩屬成員常被報導有種類/族群有長期使用閉鎖花 (cleistogamy)作為繁殖策略，降低族群內異型合子率，上述這些情形又使得局部族群常為單一遺傳組成(Sun, 1999)，再加上長期氣候變遷導致地理切割事件 (vicariance)對於黃芩屬的族群/物種分化有一定程度的貢獻，在過去的研究亦發現台灣產的黃芩分化及偶發性族群擴遷均能對應特定的地質事件(eg. Marine isotope stage, MIS)(Chiang et al., 2012a)。如此說明破碎的地理結構及長遠的氣候變遷同時作用於台灣產黃芩之分歧。以及黃芩的拓遷對長期氣候變遷之敏感。. 圖 1-1 台灣產黃芩物種地理分布簡圖 8.

(16) 1.4 本論文研究議題然而台灣島形成不過兩百萬年左右，而過去曾以核基因座與葉綠體基因座的序列分析對台灣產之黃芩屬(Scutellaria, Labiatae)植物之來源及種化時間作過推論，提出台灣原生之黃芩為多次起源，且來台後又經歷快速種化而形成高特有種比例(6/8)的假說，綜合上述，台灣產的黃芩屬成員的快速種化，單用台灣的氣候變遷造成的黃芩長距離遷移、地理切割及地理隔離來解釋其能在 60 萬年內種化至如今六種特有種之多，其中 indica group 則可在 20-30 萬年間即有三個特有種在台灣分化(Chiang et al., 2012a)，在有物種間分佈範圍大幅重疊的情況下，地理因子很難單獨解釋物種的快速分化，再加以台灣氣候異質性高，生態因子作用在族群分化甚至加速種化的生態分歧或生態種化可能參與於快速種化中。生態種化最主要的兩個切入層面是：1)分歧過程不必然地理障礙，僅需少數基因受分歧性天擇影響，故大部分基因座的基因交流依舊順暢，因此使用多基因座分析時，我們預期會看到一定程度的 historical gene flow。2)生態因子諸如氣候等等可以解釋遺傳變異或物種分化，這可以用生態棲位模擬、迴歸及相關性分析了解遺傳變異與氣候因子是否有關聯性，以間接支持生態因子作用的可能性。因此，我們擬以上面的策略，分別以三個獨立研究闡述生態因子如何維持物種和族群的分化，甚至可能加速台灣產黃芩的快速種化。表 1-1 台灣之黃芩屬植物分布區域。物種向天盞 S. barbata Don 台北黃芩 S. taipeiensis Huang, Hsiao, et Wu 印度黃芩 S. indica L.. 特有種與否分佈否台灣中、北部(少數分布於南部)、中國大陸、喜喜拉雅山區是台灣北部(台北山區海拔約 200m). 田代氏黃芩 S. tashiroi Hayata 布烈氏黃芩 S. playfairi Kudo 南台灣黃芩 S. austrotaiwanensis Hsieh & Huang 台灣黃芩 S. taiwanensis Wu 長葉黃芩 S. hsiehii Hsieh. 否是是是. 台灣中、北部 2400m 以下山區、東亞、西伯利亞、東南亞台灣東部(花蓮、台東) 台灣中、南部(嘉義、屏東、高雄) 恒春半島，陡坡或林緣、海拔 400~500m. 是是. 台灣東、南部(台東、屏東) 台灣中部 (南投). 粗體標示的物種為 indica group 的成員。. 9.

(17) 第二章台灣產近緣黃芩屬植物的分歧及種間基因交流測定 2.1 前言島嶼物種輻射種化及多次拓殖事件常造成島嶼高特有性，舉例來說 Gillespie et al. (1994)在研究夏威夷島上 Tetragnatha 屬蜘蛛時，合併來自可能來源地的其他同屬成員重建親緣關係樹來闡述這些物種單次拓殖至島上並快速輻射種化。 Papadopulos et al. (2011)選用 Lord Howe Island 作為研究樣地，以島上特有的 Coprosma 屬植物屬單一分支及相對短的分化時間推論其為島嶼快速種化的例子。而以上兩例均為海洋性島嶼，島嶼距離大陸相當遙遠，拓殖不易減少物種拓殖成功的機率，再現拓殖成功便可能性較低，因而相對容易找出單次拓殖的例子，再依地形複雜程度導致島嶼區域氣候異質性高增高區域性適應(local adaptation)，該次拓殖便可能在地大量輻射種化(Paulay, 1994; Mallet et al., 2014; Papadopulos et al., 2014)，以上的例子均以觀察到該島特有種組成單一分支便能以簡約性法則推測為單次拓殖後大量種化，且以 Lord Howe Island 的例子來看，該島面積小，物種分歧時伴隨暫時異域分歧的可能性大幅降低，但區域性特性依然存在，因此 Papadopulos et al. (2011)能更有信心推論出島內輻射種化，以及島內的輻射種化多為同域的生態種化的結論。然而大陸性島嶼，特別是在冰河時間海平面下降形成陸橋者(例如台灣對應東亞陸塊、馬達加斯加對應東非陸塊)，在冰河時期可允許大陸與島嶼物種相互遷移，姐妹種異域種化並平行拓殖至台灣(parallel colonization)以及多次拓殖時常發生，要驗證小尺度的環境異質性導致種化加速便不易回答(Coyne & Orr, 2004; Coyne, 2007)。即使台灣因地形崎嶇造成的海拔落差(~4000 公尺)而因此有了多樣的氣候因子及特有種比例(Hsieh, 2002)也不易找到這樣的例子。要了解這些島嶼近緣種間以單次拓殖並因生態因子輻射種化的例子，第一步必需重建過去這些近緣種與鄰近地區的同屬成員作材料進行生物地理分析重建祖先分佈點，例如馬達加斯加島的 Canarium 便被推論為單次由東南亞長距離拓殖而適應當地，且該屬在馬達加斯加均為特有種(Federman et al., 2015)。台灣雖鮮有單次拓殖的例子，但也有以生物地理分析發現平行拓殖，但至少一分枝具有高特有性的例子，例如佛甲草(Ito et al., 2017)或小蘗(Yu & Chung, 2014)，同時也有相反情形的例子，台灣特有的脊椎動物大多是由鄰近區域多次拓殖後冰河期隔離產生分化而形成今 10.

(18) 天的台灣動物的特有性(He et al., 2018)。但是僅靠祖先地理重建也僅能排除/考量島嶼物種分化過程並無異域時期或離開角島嶼的可能性，要證明其為生態種化主導島嶼生物的多樣性，物種間分化過程中伴隨著祖先基因交流及生態因子造成的天擇效應，以至完成生殖隔離的過程亦需檢視(Wu, 2001; Coyne, 2007; Nosil et al., 2009b; Coyne, 2011; Nosil & Feder, 2012)。在同域種化(sympatric speciation)的過程中，兩族群藉由區域性適應或分歧性天擇(divergent selection)的作用使基因組少部分區域並不允許基因交流，形成小規模的生殖隔離，而大多數區域則容許基因交流，然而基因組則以搭便車效應(hitchhiking)擴增基因組內生殖隔離的成份。漸漸減少物種間基因交流 (interspecific gene flow)的機會，最終生殖隔離完成(Wu, 2001; Rundle & Nosil, 2005; Nosil et al., 2009a; Nosil & Feder, 2012)。基因交流在遞減的過程中仍扮演著多樣的角色，多數案例認為頻繁的基因交流可以讓遺傳分化完全的物種再次產生混雜。舉例來說：外來的近緣種與原生種產生二次接觸(secondary contact)後常會產生遺傳混雜(genetic admixture)，可能致使分化中斷，產生逆向種化(reverse speciation)，並可能因此減少當地的物種多樣性(Seehausen et al., 2008)。或者可能藉二次接觸(secondary contact)產生的基因漸滲使適應性的等位基因流至另外一方(Hedrick, 2013)。利用溯祖分析 Isolation-with-migration(IM)模型可以計算過去物種分歧之間是否產生基因交流(Wakeley & Hey, 1998; Pinho & Hey, 2010)。除此之外，IM 後續的演算法 IMa2 及 IMa3 提供我們可以進一步檢視多個已知親緣關係的近緣種，在分歧過程中其基因交流旺盛的程度(Hey, 2009; Hey et al., 2018)。綜合上述，本研究選用台灣產黃芩屬其中一個特有性高的分枝作為研究材料，該分枝於 Chiang et al. (2012a)中發現的三次拓殖入台灣的分枝中，含有最高特有種比例的分枝，在該研究被暫時命名為 indica group，由廣佈於東南亞、南亞及鄰近島嶼的印度黃芩(Scutellaria indica)，同為特有種的田代氏黃芩(S. tashiroi)、布烈氏黃芩(S. playfairi)及南台灣黃芩(S. austrotaiwanensis)所組成。該類群分歧時期經估算在印度黃芩拓殖入台灣後在約莫 27 萬年內便形成現在三個特有種的情形，且這三個特有種與印度黃芩在現今分佈棲地均有或多或少的重疊，在台灣東部如花蓮太魯閣或東南部霧鹿均可見到兩個或以上的物種共域(表 2-1；圖 1-1)、這四種黃芩在過去分類上處理均極為相似，繁殖特徵僅田代氏黃芩沒有另外三種黃芩 11.

(19) 在花冠基部有 S 型折曲(Sigmoid corolla)(Hsieh & Huang, 1995, 1997)，傳粉者應為蜂類及蠅類，偶有蝶類及日行性蛾類訪花，缺乏專一性的傳粉者(Zhang et al., 2011; Wiesenborn, 2019)，及共域地常有難以辨識的小苗出現，使得這四種黃芩在遺傳上是否如其分類所示完全分化？或是大量的遺傳混雜？抑或是由生態種化所述，早期有基因交流但現今已遞減甚至中止？了解這些議題是了解黃芩屬於台灣產的特有性的其中一個關鍵之一。因此，本研究的目的羅列如下：(1) 了解台灣產黃芩的遺傳組成，並檢視種間分化情形。 (2) 檢視台灣產黃芩分歧過程是否有基因交流？. 2.2 材料與方法 2.2.1 材料收集收集的植物葉片以矽膠乾膠後保存(表 2-1)。其 DNA 參考 Doyle and Doyle (1987)的 CTAB(cetyl trimethylammonium bromide)萃取之，以 1%瓊脂電泳 (argarose)確認完整性後，以 1X TE buffer 並保存在-20 ˚C 冰箱。表 2-1、採樣地訊息及各物種及族群代號物種/代號. 地點 (代號). 座標 (經/緯度). N. 南台灣黃芩(Aus) 屏東恒春 (HC). 120.1753. 23.95333. 19. 屏東里龍山 (LLS). 120.7197. 22.17639. 31. 台南南化 (NH). 120.3611. 23.01861. 14. 南投奧萬大 (AWD). 121.1768. 23.9534. 9. 新北烏來 (WL). 121.5539. 24.93382. 16. 日本大阪生駒山 (JP). 135.6806. 34.67957. 5. 屏東浸水營 (JSY). 120.7525. 22.42263. 31. 台東霧鹿 (WU). 121.0375. 23.17164. 1. 台東蘭嶼 (LY). 121.5102. 22.07737. 2. 花蓮太魯閣 (TRK). 121.5116. 24.17946. 4. 屏東霧台 (WT). 120.7219. 22.74091. 10. 台東霧鹿 (WU). 121.1132. 23.13664. 19. 印度黃芩(Ind). 台灣黃芩(Tai) 田代氏黃芩(Tas). 布烈氏黃芩(Pla). 12.

(20) 為了評估各族群的遺傳多樣性，我們選用共顯性(codominant)的微衛星標記 (microsatellite)進行基因型讀取(genotyping)。我們在過去開發出共 19 組於台灣產黃芩屬物種間通用的微衛星引子(Chiang et al., 2012b)，我們參考其測試條件進行 PCR，並以 10% polyacrylamide gel electrophoresis (acrylamide: bisacrylamide 29:1) 進行膠體電泳，膠體電泳的結果以軟體 Quantity One v4.6.2 區分不同長度的 PCR 產物進行基因型讀取。 2.2.2 基因多樣性整理我們以 GenAlEx 6.5(Peakall & Smouse, 2012)計算各個族群與基因座的遺傳多樣性指數。並為確認是否有基因座受到天擇影響而造成後續模型演算評估有偏差，我們以 Bayescan v2.1(Foll & Gaggiotti, 2008)確認哪些基因座偏離中性， Bayescan 在計算每個基因座的 Fst 指數時，區分為族群間的 alpha 及族群內的 beta，若該基因座受到分歧性天擇，則族群間的 alpha 便大於零，藉由檢測該基因座 alpha 為零及大於零，該軟體採用 logarithm of posterior odds[log(PO)]評估兩者的模型，該值大於 1 者則代表該基因座可能為受到分歧性天擇作用，而 Bayescan 設定如下，一開始共計 20 次 pilot runs，每個 pilots runs 計有 5,000 筆資料。先移除前 50,000MCMC steps 作為 burn-in，接著每 10 步取樣一次，記錄 10,000 筆資料，共計 150,000MCMC steps。 2.2.3 遺傳組成及基因交流檢測過濾掉適應性基因座，其餘基因座分別進行 STRUCTURE(Pritchard et al., 2000)與 R 套件的 adegenet 模組 DAPC(Jombart, 2008) (discriminant analysis of principal components)來了解各物種內每個個體的遺傳組成及混雜程度。其中 STRUCTURE 假定族群為理想族群(ideal population)並以各基因座的等位基因頻率評估個體組成機率，而 DAPC 則是併用 discriminant analysis(DA)及 PCA 兩種作法進行 ordination analysis，以最能解釋變異的前兩組成份(principal components, PCs)來區分各別物種或族群。我們同時使用兩者可以更區分台灣產黃芩間的遺傳組成及混雜程度(Manel et al., 2005; Jombart et al., 2010)。而 STRUCTURE 則是以 admixture model，並不納地理位置作為考量(without LOC prior)進行對每個群聚值 13.

(21) (K)進行 10 次獨立的模擬，每次模擬包含 100,000 次的 burn-in 及 1,000,000 次的 MCMC 模擬。每個群聚值的十次獨立模擬的共識以 CLUMPP (Jakobsson & Rosenberg, 2007)進行統整，確認共識值 H’>0.9 後輸出獨立模擬的最後共識值，若 H’值低於 0.9，則將各個 independent runs 分別檢視，並移除 log liklihood 過低的 indepdent runs (remove outlier)再進行 CLUMPP 分析。而 DAPC 則以 BIC 值找出最好的 cluster 數，並利用 find.clusters 指令將基因型轉成 PCs，再以 dapc 指令將上述的 PCs 進行 DA，將結果以 scatter plot 及 individual assignment plot 的方式呈現個體混雜程度。除此之外，為了計算黃芩屬過去的基因交流程度，我們使用 isolation-with-migration (IM)演算法計算種與種之間的 historical gene flow。其中 IMa3 運用 Bayesian MCMC 與 likelihood 合併計算參數(m, theta, and t)，原理大致與 IMa2 相同，然而，IMa3 與 IMa2 不同的地方在前者可設立一個 ghost species 作為外群(Hey et al., 2018)，主因為本研究使用台灣黃芩(Tai)作為外群，然而台灣黃芩遷入台灣應為獨立於 indica group 的事件(Chiang et al., 2012a)，使得我們無法確認是否有滅絕或異地的近緣種也曾與我們的採樣族群有基因交流，並為此錯誤計算基因交流及共祖的有效族群大小(Slatkin, 2005; Hey et al., 2018; Lynch & Sethuraman, 2019)。因此我們也依照 IMa3 指示加入 ghost species 作為外群以減少 gene flow 錯誤估算。 2.3 結果 2.3.1 遺傳多樣性經 Bayescan 檢測後發現所有基因座的[log(PO)]均小於 0，顯示這 19 個基因座並無偏離選汰中性。故 19 組基因座均適合放入探討純粹的族群動態 (demography)分析中而暫時不需考慮天擇對各別基因座的影響。經計算各族群內的等位基因數平均值可發現，族群內帶有的等位基因數不多(1.05±0.05 2.79±0.44，表 2-1)，同時也發現族群的異型合子觀測值都低於預期值(表 2-2)。這些情況都反應在偏高的近交係數上(inbreeding coefficient)。與其他黃芩屬成員自交親合及小型蜂類授粉的情況一致。除此之外，Shannon index 在物種的階層普遍偏高(0.22±0.07 - 0.85±0.20，表 2-2)，卻在族群的階層相對低很多(0.04±0.04 0.60±0.13，表 2-2)，顯示黃芩屬底下各種內族群間的α-diversity 異質性很高。這些數據也與 AMOVA 的結果一致，但是大部分族群內的遺傳多樣性指數常常遠 14.

(22) 低於族群合併後的值，這暗示族群內的變異相對小且異質性低。表 2-2 五種黃芩及其族群的遺傳多樣性指數計算平均值及其 95%信賴區間。 Species/populations N. Na. Ne. I. Ho. He. uHe. F. S. austrotaiwanensis 64 屏東里龍山 31 台南新化 14. 2.05±0.22. 1.53±0.12. 0.44±0.09. 0.05±0.03. 0.27±0.06. 0.27±0.06. 0.84±0.05. 1.53±0.16. 1.28±0.09. 0.24±0.07. 0.05±0.03. 0.15±0.05. 0.16±0.05. 0.71±0.10. 1.16±0.12. 1.13±0.10. 0.09±0.06. 0.03±0.02. 0.06±0.04. 0.06±0.04. 0.43±0.04. 19. 1.16±0.09. 1.09±0.06. 0.08±0.05. 0.06±0.05. 0.05±0.03. 0.05±0.03. 0.01±0.21. 30. 2.53±0.23. 1.80±0.15. 0.64±0.09. 0.03±0.01. 0.38±0.05. 0.38±0.05. 0.93±0.03. 台北烏來. 16. 1.74±0.23. 1.35±0.11. 0.31±0.09. 0.04±0.02. 0.18±0.05. 0.19±0.05. 0.83±0.06. 南投奧萬大. 9. 1.37±0.14. 1.33±0.12. 0.23±0.08. 0.00. 0.16±0.06. 0.17±0.06. 1.00. 日本大阪生駒山. 5. 1.05±0.05. 1.05±0.05. 0.04±0.04. 0.05±0.05. 0.03±0.03. 0.03±0.03. 1±0.03. 29. 2.95±0.34. 1.86±0.19. 0.65±0.11. 0.11±0.03. 0.37±0.06. 0.37±0.06. 0.71±0.05. 屏東霧台. 10. 2.53±0.32. 1.87±0.20. 0.60±0.12. 0.18±0.05. 0.34±0.07. 0.36±0.07. 0.43±0.09. 台東霧鹿*. 19. 2.00±0.28. 1.55±0.19. 0.40±0.11. 0.07±0.03. 0.23±0.06. 0.24±0.06. 0.73±0.05. 屏東佳樂水 S. indica. S. playfairi. S. taiwanensis 屏東浸水營 S. tashiroi. 31. 1.37±0.11. 1.31±0.11. 0.22±0.07. 0.04±0.02. 0.16±0.05. 0.16±0.05. 0.78±0.05. 68. 5.21±1.14. 2.64±0.45. 0.85±0.20. 0.15±0.04. 0.38±0.08. 0.38±0.08. 0.56±0.05. 台東霧鹿*. 22. 3.16±0.73. 1.92±0.33. 0.54±0.17. 0.10±0.03. 0.26±0.08. 0.26±0.08. 0.66±0.04. 台東蘭嶼. 15. 2.16±0.43. 1.74±0.28. 0.45±0.13. 0.13±0.06. 0.25±0.07. 0.26±0.07. 0.6±0.09. 花蓮慕谷慕魚. 16. 2.79±0.44. 1.86±0.20. 0.60±0.13. 0.20±0.05. 0.33±0.07. 0.34±0.07. 0.36±0.05. 花蓮太魯閣. 15. 2.63±0.47. 1.88±0.25. 0.58±0.14. 0.20±0.05. 0.32±0.07. 0.33±0.07. 0.39±0.06. *共域族群 N: 樣本數; Na: 等位基因數; Ne: 有效等位基因數; I: Shannon 指數; Ho: 異型合子觀測值; He = 異型合子期望值; uHe = 校正後之異型合子期望值; F: 族群固定指數，又名近交係數. 2.3.2 物種遺傳結構分析兩種群聚分析在分群數上並不相同，利用貝氏群聚分析計算的△K 值顯示分群數 K 為 5 最佳，STRUCTURE 在 K=5 時，台灣產的特有種的遺傳混雜非常少，均小於 10%，而台灣產的印度黃芩有一小部分與布烈氏共享的遺傳組成，但這組成卻不見於日本的族群(JP)。對比之下，依照 BIC 的結果，DAPC 在 8-10 個 clusters 時 BIC 值己無明顯變化(圖 2-1)。兩種分析均顯示物種間遺傳混雜的個體並不多，僅有一個位於屏東霧台的布烈氏黃芩個體在 DAPC 分析上被歸為與印度黃芩相似的群聚中, 然而該個體在 STRUCTURE 卻與其他同為布烈氏黃芩的個體歸在一起(圖 2-2A)。經檢視 DAPC scatter plot，發現兩物種的分群非常接近不易區別 (圖 2-2B)，加以屏東霧台(WT)當地並無印度黃芩的目擊及標本記錄，初步排除該個體為近代雜交或回交個體的可能性，應為錯誤 assignment，或是印度黃芩與布烈氏黃芩有共享較多的等位基因。另外，相對 STRUCTURE，DAPC 顯示了南 15.

(23) 台灣黃芩(Aus)及田代氏黃芩(Tas)有較明顯的族群分化但帶有部分的遺傳混雜，三個南台灣族群均有自己優勢的遺傳組成成分，而霧鹿(WU)的田代氏黃芩遺傳組成則明顯異於其餘族群，離島的蘭嶼(LY)則相較混雜。相對檢視貝氏群聚法 K=10 的結果，其共識值並不高(H’ = 0.76)，因此，我們分別檢視各個 K=10 時各個獨立運算的結果，多數結果確認南台灣黃芩的新化族群 XH 與其餘南台灣黃芩有明顯分化，以及田代氏黃芩的霧鹿(WU)族群的遺傳組成的確異於其他族群，但僅有兩個獨立運算支持田代氏黃芩與南台灣黃芩有遺傳混雜，將 log likelihood 最高最接近的前四名作共識後得到的圖也未見種間遺傳混雜(圖 2-3)，因此排除 WU 族群的異質性是肇因於物種雜交而成(圖 2-2、圖 2-3)。總體而言，不管是 DAPC 或是 STRUCTURE 均指出這五種黃芩屬成員在遺傳上已有相當明顯的分化，且沒有證據顯示這五種黃芩有近代的雜交，暗示這五種黃芩之間已有相當的生殖隔離。除此之外，部分台灣產黃芩屬成員(如 Aus 及 Pla)在族群間有明顯的族群結構，顯示族群間的連接性(connectivity)並不明顯。. 400 250. 300. 350. BIC. 450. 500. 550. Value of BIC versus number of clusters. 0. 20. 40. 60. 80. 100. Number of clusters. 圖 2-1 以 R 套件 vegan 內的 find.cluster 指令算出各個以五種黃芩的 SSR 遺傳資料計算出來 cluster 數對應的 BIC 值，BIC 變化轉折處約莫在 10 左右。. 16.

(24) 圖 2-2 五種黃芩遺傳組成及群聚分析。(A) 利用 DAPC(上)及 STRUCTURE(中、下)進行的 individual assignment test 的結果，族群及物種代號見表 5-1。(B) 利用 DAPC 的結果，以 BIC 值的變化約略在 10 個 cluster 為最好的分群數，以此繪製的散佈圖。. 圖 2-3 檢視 STRUCUTURE K=10 的結果發現該群聚的 log likelihood 約略分成兩群，重新以 log likelihood 較高的獨立運算重新取得共識，相較於圖 2-2A, K=10 的結果，未發現明顯的種間遺傳混雜，但有更明顯的族群間遺傳分化情形。 2.3.3 基因交流計算在考慮來有未採樣物種的情況，IMa3 偵測到祖先的基因交流(historical gene flow)(圖 2-4)，雖然 indica group 在 STRUCTURE 及 DAPC 分析中已顯示現今的 indica group 的四個物種是生殖隔離相當完整的物種，我們仍在共域分佈的田代氏黃芩與布烈氏黃芩發現有相當明顯的祖先基因交流(historical)，布烈氏黃芩與南台灣黃芩的共祖也與田代氏黃芩有相當明顯的基因交流。即便是現今物種分佈重疊最小的南台灣黃芩(見於恒春半島，偶見於台灣西南部阿里山山脈以南)與印度黃芩(主要為台灣西部與東北部)亦有相當明顯的單向基因交流。唯 indica group 17.

(25) 三個特有種的共祖與印度黃芩沒有偵測到明顯的基因交流(圖 2-5)。總結而言，不論是 indica group 或是台灣黃芩，彼此分歧之時均伴有相當程度基因交流，這些基因交流發生的時間，若以微衛星體保守估計為 1 x 10-6 per generation per locus 計算，indica group 的分歧伴隨著基因交流的時間可以從 19.5kya 一直到 1.5kya(圖 2-5)均有明顯的祖先基因交流。這結果展現不同於 STRUCTURE 與 DAPC 無種間混雜的結果。. 圖 2-4 以 IMa3 計算出各個分支間基因交流的程度。樹型及分化時間由於無確切參考的微衛星體突變率，故參照 Chiang et al. (2012a)之樹型及分化時間(kilo years ago, kya)。箭號及其大小代表基因交流的相對強度及方向，並以紅色代表 indica group 內的基因交流，而藍色是至少有一端為非 indica group 成員的基因交流。 2.4 討論 2.4.1 更新世末期的台灣產黃芩具有旺盛的種間基因交流本研究以台灣產黃芩 indica group 現今遺傳多樣性組成、溯祖計算其基因交流程度，闡明 indica group 在輻射快速種化中依然有旺盛的基因交流，在本篇研究也發現雖然現今 indica group 已形成強烈的生殖隔離，甚少種間混雜(圖 2-2)， 18.

(26) 然而溯祖類分析(IMa3)在更新世末期間(19.5kya, 圖 2-4)，indica group 也可能因為分化或遷移的過程中產生或中斷基因交流，舉例來說，現今常有共域族群或鄰域族群(圖 1-1)indica group 三個特有種之間便有旺盛的基因交流，然而，在初拓殖入台灣時，indica group 的共祖，可能藉由黃芩種子可藉偶發性的洪水產生長距離遷移的能力(Williams, 1992; Cruzan, 2001; Nelson & Goetze, 2010)，從古閩江 (Paleo-Mingjing)拓殖入台灣北部，其中印度黃芩往台灣西方擴散，而其餘三種特有種的共祖則往台灣東部擴散(Chiang et al., 2012a)，而台灣中間夾著雪山山脈及中央山脈，台灣缺乏可橫貫東西部的河流系統，黃芩授粉能力及種子傳播能力非常受限(Williams, 1992; Cruzan, 2001; Nelson & Goetze, 2010; Zhang et al., 2011)，不容易跨越南北走向的山脈系統，也因此能解釋為何在圖 2-4 中 t1 至 t2 間無觀察到明顯的印度黃芩與三特有種共祖的基因交流。但三特有種及其共祖彼此之間一直偵測到頻繁的基因交流。 2.4.2 族群間具有明顯結構在 DAPC 與 STRUCTURE K=10 的結果中，三個特有種都展現了有族群分化的趨勢，南台灣黃芩的新化族群(XH)明顯不同於地緣性較強的恆春半島(LLS & HC)族群，布烈氏黃芩在東部(WU)跟西部(WT)兩族群主要組成也有異，田代氏黃芩雖然族群間混雜程度較前兩者略高些，但仍可見到霧鹿(WU)與蘭嶼(LY)有自己主要的組成，這也反應在族群的多樣性指數不高，但族群固定指數卻非常高 (F 值, 表 2-2)。暗示族群間的基因交流並不旺盛。要解釋為何族群結構強，可以由地理因子、黃芩本身的遷移歷史及環境因子來推論，如同上述，黃芩的傳粉機制主要由小型蜂類與蠅類，有效傳播傳播距離短(Zhang et al., 2011; Wiesenborn, 2019)，易產生近交，其中這三種黃芩都可藉走莖進行無性繁殖，以及黃芩有閉鎖花(cleistogamy)現象會促進自交(Sun, 1999; Denda, 2002)，又可能僵化族群內的遺傳組成，造成高比例的同型合子(表 2-2)。台灣島南北狹長，然而主要的溪流多為東西走向，要藉溪流進行長距離遷移需靠河流或洪水才能產生跨平原的移動，然而這三種特有黃芩主要都居於中低海拔山區(Hsieh & Huang, 1995, 1997)，高地致使難洪水現象不易波及而難以造成南北遷移，東西橫貫又受山脈阻絕，綜合這些繁殖系統(mating system)以及台灣的地理，. 19.

(27) 便可對台灣黃芩屬成員產生一定的阻絕性，使得族群間雖然具有遺傳混雜但還是能觀察到明顯的遺傳結構。但是，若完全以地理阻隔解釋台灣產黃芩屬的快速種化仍有疑問。首先，我們的研究指出 indica group 在分化時有大量的基因交流，大量的基因交流理應促使物種之間的差異性下降，減緩分化的速度甚至倒轉整個種化的過程(reverse speciation)(Seehausen, 2006)。若黃芩屬間地理隔離是主要的因子，這些祖先基因交流的結果似乎與其衝突。可能的解釋為：1)在 Chiang et al. (2012a)利用 DIVA(Dispersal-Vicariance Analysis)-based 分析重建 indica group 拓殖台灣的路徑之時，發現冰期-間冰期的氣候變遷促成 indica group 的祖先產生偶發性的長距離遷移並伴隨著地理切割事件(vicariance)。我們推論這些長距離遷移及地理切割事件改變了 indica group 祖先的分佈範圍，使現今相互交疊的 indica group 物種在過程曾有隔離，在這些黃芩屬花部特徵及傳粉者近乎無差異(Hsieh & Huang, 1995, 1997; Zhang et al., 2011; Hsieh, 2013; Wiesenborn, 2019)的情況下，加以部分區域如霧鹿(WU)有觀察到部分無法區分為布烈氏黃芩或田代氏黃芩的小苗(疑似雜交個體，兩者成體葉長寬比相差甚多，小苗介於中間，個人觀察)，這些情況可能說明 indica group 的合子前隔離機制並不完整(prezygotic isolation)，不僅如此，有部分近緣種黃芩屬成員人工雜交後種子仍有萌發能力，雜交的 F1 產生的花粉仍有部分活性(Roland Von, 1991)，這也說明黃芩屬成員經常會有合子後隔離 (postzygotic isolation)尚未完全的情況，因此，與親本回交產生基因漸滲便是可能的。這時候如果有一段時期是異域分佈的，由於傳播能力不強，種間基因交流的頻率便可能因此下降。2)即便沒有因氣候變遷產生的黃芩屬南渡北遷造成分佈重疊現象，indica group 的成員在分化時也可能因台灣島破碎崎嶇的地形導致強烈的氣候及棲地異質性，這些異質性可能對族群內產生分歧性天擇(divergent selection)，或對族群間產生區域性適應(local adaptation)，使得不同個體/族群基因組部分區域產生隔離，即便其他區域仍有基因交流(porous genome)(Wu, 2001)，直到這些隔離區域因搭便車效應(hitchhiking)在基因組擴張，這樣的生態種化 (ecological speciation)便能解釋為何現今 indica group 成員無論在形態及 STRUCTURE 分群結果等遺傳資料是可清楚區別的群體(物種)(表 1-1；圖 2-3)，卻又能在溯祖分析(IMa3)中偵測到大量祖先的基因交流(圖 2-4)。 20.

(28) 因此，我們推論生態因子諸如生物氣候(bioclim)等，可以減少物種間的基因交流，並維持物種的分化，甚至與地理隔離一併加快這些物種的分化。要證明這假說，我們必須取得族群間的地理資訊、遺傳資料以及環境因子等，進行溯祖分析確認近緣種/姐妹種間的在分歧時基因交流改變的趨勢，並將環境因子在控制地理因子的情況下與遺傳因子進行迴歸性分析(regression analysis)確認其關聯性。然而，台灣產黃芩共計八種，至少形成三個獨立的分支(Chiang et al., 2012a)，在溯祖分析上要 IMa3 的 tutorials(https://github.com/jodyhey/IMa3)即不建議在三個種以上的情況進行過度複雜的探討，物種過多在進行迴歸性分析要控制的變量及不確定性也會隨之增多，因此，我們在下兩個章節，會分別選擇兩組近緣並有明顯共域/鄰域族群的台灣黃芩屬成員，分別是田代氏黃芩與布烈氏黃芩組，以及台北黃芩和向天盞組，闡述環境因子如何影響物種及族群的遺傳分化，來確認生態種化或生態分歧(ecological divergence)是否能解釋台灣產黃芩的多樣性。. 21.

(29) 第三章島嶼上的空間環境異質性趨使島嶼物種持續分化以致生態種化 3.1 前言分歧性天擇不僅可以加速種化，也能持續的維物種的分化。早期談到物種分化大多著眼於地理障礙或地理距離(isolation-by-distance, IBD)，然而近年愈來愈多的研究指出環境變異比地理距離更能解釋物種或族群的遺傳分化(Coyne, 2011; Shafer & Wolf, 2013; Sexton et al., 2014)，這些研究都指出周遭環境變化對物種演化的重要性。Via (2009)提出天擇在生態種化上主要參與於兩個階段，一個是藉分歧性天擇阻擋基因交流，此階段可以區分異域種化(allopatric speciation)及同域種化(sympatric speciation) (Wu, 2001; Via, 2009)。另一階段則是在物種分化後，各別物種因面臨的生境特性而受到獨有的天擇篩選作用，這些獨有的天擇作用間接降低物種雜交的可能性，進而維持分化狀態。(Palumbi, 2009; Via, 2009; Zitari et al., 2012)。然而第一階段並不容易驗證，因為種化當下的生態因子難以重建，因此要證明生態種化(ecological speciation)相當困難。但今時的環境因子可藉地理資訊系統所提供的圖層取得，因而仍能驗證 Via (2009)所提的種化第二階段，也就是生態因子對種化後分歧的影響。另外，由於適應性分歧的種化過程，種間基因交流會依分化時間隨之遞減的特性 (Wu, 2001; Scopece et al., 2007)，我們可以藉由溯祖(coalescence theory)探測過去種間基因交流的變化來檢測兩物種分歧是否吻合生態種化的假說。再者，藉由天擇的影響，可以比單純藉遺傳漂變(stochastic drift)更快加速物種(或族群)的分歧(Slatkin, 1987; Ellstrand & Elam, 1993)，只要這些天擇持續維持甚至加強物種間的分歧，物種間基因交流會更快速的遞減。目前以建立模型為基礎的 approximate Bayesian computation(ABC)來評估物種間分化與基因交流已臻成熟。我們可以利用 ABC 來評估各種演化情境，因此我們檢測完全無基因交流、基因交流遞減或持續性的基因交流等情境(scenario)哪一個最能解釋我們的遺傳資料，以了解種間基因交流的變化。另外，配合環境因子與遺傳變異進行相關性分析，我們可以知道環境異質性及區域性適應如何驅動物種分化。台灣產黃芩中，布烈氏黃芩(S. playfairi)與田代氏黃芩(S. tashiroi)相當適合作此議題的材料。兩者均為台灣特有種。布烈氏黃芩主要分佈於台灣西南部，並一路延伸至中央山脈東南側，而田代氏黃芩主要分佈於中央山脈東側、東南側以及 22.

(30) 海岸山脈，兩者有部分族群於中央山脈東南側共域(圖 3-1)，加以兩者外觀相當相似，主要區別的地方有 1)布烈氏黃芩花冠筒基部呈下折後 90 度斜上生長 (sigmoid corolla)，而田代氏黃芩僅為單純斜上生長。2) 布烈氏黃芩上下唇瓣均為淺白至淺藍色，偶有淺粉紅色個體，然而上下唇瓣的顏色並無明顯區別。相反的，田代氏黃芩下唇瓣為深藍色，其餘花冠筒部位則為白色到淺藍色，3)布烈氏黃芩為頂生總狀花序，田代氏黃芩鮮為頂生總狀花序，常為腋生總狀花序 4)布烈氏黃芩好生於森林底層或略有遮蔭的林緣岩坡下，而田代氏黃芩則喜歡生於開闊垂直聳立的峭壁縫中，也因為有此一微棲地的差異，田代氏黃芩常懸垂於懸崖上受陽光直晒，而布烈氏黃芩植株則相對挺直於林緣有蔭蔽的岩坡邊(Hsieh & Huang, 1995) ，然而在共域處可見到難以鑑定的小型植株，暗示兩物種關係之近。兩物種的分化時間相當短，僅約有 0.2 Mya，而且台灣產黃芩屬成員的花色調節基因如 myb gene family 及 anthocyanin biosynthesis pathway 也受到的正向天擇作用 (Huang et al., 2015; Huang et al., 2016b, a)，加以這兩種所屬的 indica group 在更新世後期拓殖到台灣後，各自物種的遷移-地理切割 (dispersal-vicariance)事件都與氣候變遷息息相關(Chiang et al., 2012a)，這些跡象都暗示這兩個物種的 post-speciation divergence 與生物地理事件、氣候變遷及環境異質性導致的天擇有關。生態因子例如地理-環境異質性，產生的分歧性天擇或區域性適應能有效的解釋物種分化及族群結構，因此，我們可以藉評估環境因子在影響物種分佈及遺傳結構的重要性，以進一步了解生態因子如何驅動並保持物種分化。舉例來說，線性相似矩陣(linear dependent similarity matrix)一類的多變量分析可以評估控制因子(constraint factors)如何影響應變因子(response factors)，藉此了解地理距離或環境距離如何解釋遺傳距離(Balkenhol et al., 2009; He et al., 2013)。除此之外，利用氣候因子進行物種生態棲位模擬(ecological niche modeling, ENM)也可以探索哪些環境因子與該物種的分佈有關，哪些關鍵因子在物種間有分化，這些都常被用來檢測生態種化與分歧(Glor & Warren, 2011; Blair et al., 2013)。再加上 ABC 的演化假說推演分析(hypothetico-deductive evolutionary analyses)更針對各個假說取得其描述式統計(summary statistics)，在不需 exact likelihood calculations 的狀況下，確立最好的種化模型。(Bertorelle et al., 2010; Csillery et al., 2010)。這樣的流程已 23.