第一節 全基因組關聯性分析 (GWAS)
1. 外表型調查
本次 GWAS 的 240 個品種隨機分為 4 個批次種植，1 批次內含 60 品種。但 從冷處理組的原始資料的分佈情形可見 (附圖 1)，批次間存在效應，導致各批 次的平均值有所差異：第二批表現優於其它三批次，第三批的新根表現偏差，
第四批的老根表現也偏差。各批次皆含60 種原，如此的數量已足夠龐大，因此 我們假定各批次間的材料不應有偏向，意即若存在批次效應，不會是取樣偏差 所造成的。面對批次效應，我們採取的策略是以平均值對批次間進行校正，將 各批次的平均調整至同一水平 (附圖 1)。
冷處理組的原始資料經校正後，以直方圖觀察 TRls、TRn、NARls、NARn、 OARls及OARn等六個性狀的分佈狀況 (圖 5)。六個性狀皆呈連續分佈，並無明 顯的左偏或右偏。如此結果支持我們以關聯性定位為主要方法，尋找與目標性 狀具高度關聯之遺傳區域。
六個性狀彼此也存在一定的相關性，相關係數介於 0.50 ~ 0.96 (表 2)。唯新 根性狀 (NARls、NARn) 與老根性狀 (OARls、OARn) 間相關係數皆不顯著，這 說明新根的出生與老根的繼續生長，很可能是由不同的生理機制所控制。另 外，無論是總根、新根或老根，其總長性狀皆與總數性狀呈高相關 (0.82、
0.83、0.96)，總長性狀本已就包含總數性狀的意義在內，高相關的結果說明，
我們可能可單獨以總長性狀作為三種根的耐冷指標。
在此我們關心的是地下部的總伸展面積，也因為 TRls與其它五個性狀皆有 一定的相關性，往下的分析以TRls為優先檢視的性狀，其它性狀為輔。
圖 5 六個性狀的直方分佈圖
圖中六個性狀的排列，由左上至右下分別是TRls、NARls、OARls、TRn、NARn
及OARn。
表 1 6 個性狀間的相關係數
TRls TRn NARls NARn OARls OARn
TRls 1.00 0.82*** 0.79*** 0.64*** 0.58*** 0.50***
TRn 1.00 0.58*** 0.75*** 0.63*** 0.66***
NARls 1.00 0.83*** -0.03 -0.06
NARn 1.00 -0.01 -0.03
OARls 1.00 0.96***
OARn 1.00
** p-value < 0.001
2. 次族群結構分析
將台灣品種的資料扣除後，繼續分析本次 GWAS 中屬於 44Kpanel 的 237 個種原，其中一個種原 (NSFTV.ID = 2) 因無法查得其基因型資料，故此後的討 論亦將其捨棄。水稻具有強烈的族群結構，前人研究中提到主要可分為5 種次 族群 (subpopulation)，分別是：芳香稻 (Aromatic, Aro)、aus (AUS)、秈稻 (indica, IND)、溫帶型稉稻 (temperate japonica, TEJ)、熱帶型稉稻 (tropical japonica, TRJ) (Garris et al., 2005)。以軟體 STRUCTURE 將所有種原分群，預設 群數設為5 (K = 5)，若無法與任一次族群具 80% 相似者，則歸類為 Admix
個PC 總共解釋 55.78% 的總變異量，這樣結果大致同於 Zhao et al. (2011) 所 述，因此本次GWAS 採用前 5 個 PC 作為族群結構的估計值。
表 2 次族群變因對各性狀的 VE
Trait TRls TRn NARls NARn OARls OARn
VE (%)* 24.93 31.02 19.89 22.23 14.63 15.65
*VE (variance explained)：變方分析中將次族群平方和除以總平方和
圖 6 各次族群種原 TRls
盒鬚圖中，每盒鬚中的粗線代表中位數 (Q2)；盒的上下兩線分別是第三四分位 數 (Q3)、第一四分位數 (Q1)，兩者差值為四分位距 (interquartile range, IQR)；
盒外的上下兩線則代表Q3 + 1.5 IQR、Q1 – 1.5 IQR，作為是否離群的界線，於 兩者以外的數值視為離群值 (outlier)，以圓圈表示。以灰色代表 ADM、紫色代 表Aro、黃色代表 AUS、紅色代表 IND、藍色代表 TEJ、淺藍色代表 TRJ，橫 虛線為所有種原之平均。
圖 7 236 個種原於前 5 個 PC 的投影狀況
括號內的百分比為該PC 對整體資料的解釋變異程度。灰色代表 ADM、紫色代 表Aro、黃色代表 AUS、紅色代表 IND、藍色代表 TEJ、淺藍色代表 TRJ。
3. 關聯性定位
模式建構的好壞，可藉由百分位圖 (quantile-quantile plot, qqplot) 檢視。假 定控制一性狀的基因座不應過多，SNP 之回歸係數的 P-value 應服從均一分佈 (uniform distribution)。那麼以均一分佈的變數之百分位數為期望值，眾 P-value 之百分位數為觀測值，兩者所繪製的關聯散佈圖，即百分位圖，其不應偏離斜 域。曼哈頓圖為SNP 與其回歸係數的 P-value 的散佈圖，為方便觀察，P-value 取以10 為底的負對數值 ( –log10(P-value) )，而不同染色體的 SNP 以不同顏色 表示。顯著標準的訂定則參考錯誤發現率 (false discovery rate, FDR)，同樣使用 44K panel 作定位族群 (mapping population) 的數篇前人研究都指出，當顯著標 準訂於 –log10(P-value) > 4 時，可使 FDR 低於 25% (Famoso et al., 2011; Zhao et al., 2011; Ueda et al., 2015)，即發現顯著的 SNP 裡，其真正為不顯著的數量少於
25%。本次建立的模式，TRls之FDR 計算的結果大致同於上述等研究 (附表 個SNP (id6010626, 20.37 Mb / id6010718, 20.49 Mb / id6010812, 20.70 Mb) 在 TRn的P 模式達顯著水準， id6010718 在 TRn的P+K 模式、id6010812 在
(PC1) 作為共變項 (Famoso et al., 2011; Hsu and Tung, 2015)。定位結果如圖
TRls的P 模式下，達顯著水準 –log10(P-value) > 4 的 SNP 共有 24 個。以最 SNP (id2007347) 的差異。依上述觀察可將各種 haplotype 分為兩種集團：與 Hap 2 相近的溫帶型稉稻集團，有 Hap 1 及 Hap 2；與 Hap 2 有所差異的非溫帶
indica 亞種的種原，24 個秈稻、11 個 aus，說明 indica 亞種內可能也有具耐冷 特性的遺傳區域，只是不如溫帶型稉稻那樣普遍，效應也不如溫帶型稉稻強 烈：因為溫帶型稉稻經其生長環境的淘汰，普遍具有耐冷特性；秈稻的生長環 境並沒有這樣的自然選拔，所以耐冷特性不如溫帶型稉稻普遍，但也不代表完 全沒有。此耐冷特性很可能來自秈稻次族群而非aus 次族群，因為秈稻的 TRls
同樣表現不錯，aus 表現不佳 (圖 6)。從這樣的結果也可說明，不同次族群間 (溫帶型稉稻 vs. 秈稻)，其控制耐冷特性的 QTL 應不相同。而也因秈稻次族群 的耐冷特性不夠普遍、效應不夠強，我們無法在全族群的P 模式下找到秈稻次 族群與耐冷特性有關的QTL。
(a) P 模式
(b) P+K 模式
圖 8 兩種模式下，TRls的曼哈頓圖及百分位圖 (a) P 模式、 (b) P+K 模式
圖 9 175 種原間，Hap 2 種原與其餘種原的表現 紅色代表Hap 2，藍色代表其餘非 Hap 2 的種原。
(a) 全種原，P 模式
(b) 溫帶型稉稻 (TEJ)，P 模式 (c) 熱帶型稉稻 (TRJ)，P 模式
(d) 秈稻 (IND)，P 模式 (e) aus (AUS)，P 模式
圖 10 P 模式下，全種原的與各次族群的曼哈頓圖
表 3 GWAS 達顯著水準之 SNP 所涵蓋範圍有與前人文獻中 QTL 相重疊者 Traits Chr. SNP range
(Mb)
Reference QTL QTL pos.
(Mb)
Associated Traits Stress All 236 accessions (P model)
NARls, NARn 1 8.59 ~ 8.72 Horii et al. (2006) - 3.43 ~ 10.75 RAL (the length of the main root axis) 30 days TRls, TRn,
NARls, NARn
2 18.10 ~ 18.44 Courtois et al. (2003) - 10.26 ~ 20.73 DRS
(deep root weight per unit shoot weight)
45 days / Drought
TRls, TRn 4 6.54 ~ 7.51 Andaya and Mackill (2003) qCTS4-1 0.68 ~ 6.57 CT (cold tolerance，leaf yellowing) 28 days / Cold TRls, TRn 4 6.54 ~ 11.87 Courtois et al. (2003) - 6.57 ~ 11.29 DRW (deep root weight)
DRT (deep root weight per tiller)
45 days / Drought
TEJ accessions (P model)
TRls 2 18.10 ~ 19.44 Courtois et al. (2003) - 10.26 ~ 20.73 DRS
(deep root weight per unit shoot weight)
45 days / Drought
IND accessions (P model)
TRls 1 41.69 ~ 41.70 Zheng et al. (2006) - 40.29 ~ 41.92 ARN (adventitious root number) 14 days / Drought
表 4 第 2 條染色體上 18.20 ~ 18.60 Mb 間的 haplotype 分析
約400 kb 的區間內，14 個達顯著水準的 SNP，能將各種原分為 7 種 haplotype。任何種原在此 14 個 SNP 若有缺值者，則視為未歸類 (undefined)。
SNP
1 2 3 4 5 6 7 8 9 10 11 12 13 14
ud2001000 ud2001009 ud2001011 id2007257 id2007258 id2007268 id2007273 id2007278 id2007296 id2007301 id2007312 id2007317 id2007323 id2007347
Subpopulation Phenotype
Physical position
(bp)
18222134 18252287 18305871 18307289 18308009 18344580 18346536 18347547 18357037 18359019 18383771 18385072 18407944 18442319 ADMIX
Japonica Indica Total Mean of TRls
(S.D.)
Hap # Aro TEJ TRJ AUS IND
1 G T T T A G C C T G C G T A 1 0 0 0 0 0 1 22.51 (NA) 2 G T T T A G C C T G C G T C 5 5 50 1 0 0 61 31.29 (7.64) 3 T C C A G A G T C T T A A A 18 0 5 46 0 1 70 24.25 (6.16) 4 T C C A G A G T C T T A A C 0 0 0 0 0 1 1 41.90 (NA) 5 T C C T G A G T C G T G A A 2 0 0 0 11 24 37 28.20 (10.79) 6 T C C T G G G T C G C G T C 0 0 0 0 2 0 2 20.84 (3.26) 7 T C C T G G G T C G T G T C 1 0 0 0 2 0 3 28.07 (3.69)
8 Undefined 8 0 7 7 23 16 61 26.25 (8.43)
第二節 親本挑選與外表型驗證
第 2 條染色體上 18.20 ~ 18.60 Mb 候選區域間的 haplotype 分析，找到 Hap 2 此表現良好的 haplotype，其有利於接下來雜交親本的選擇。在此階段對候選 親本再進行一次外表型的驗證，增加各種原的重複數，以確認其耐冷特性的穩 定程度。驗證階段有16 個候選親本 (表 5)，選擇的依據分述如下：表現良好的 種原，其TRls排於236 個種原的最前 40 名，且屬於溫帶型稉稻、Hap 2，並來 自高緯度的溫帶地區 ( > 40˚ N /S)，共 9 個；但秈稻也有發現幾個表現不錯的 種原，故同樣保留2 個名額；表現不良的種原，其 TRls排於最末40 名，不屬 於溫帶型稉稻、Hap 2，共 5 個。候選各種原的詳細資料則列於附錄表 2。
16 個候選種原的驗證結果如圖 11。批次間的差異僅在第一批的對照組 (B1-01 ~ B1-03)，與其它二批的對照組相比稍為浮動，其餘表現皆相似；盆裝 間的差異也不大。
GWAS 時採用冷處理下的觀測值作遺傳定位，在此也就冷處理組進行討 論。本次外表型驗證時每盆內的試驗單位完全相同，皆為16 種原、3 重複，因 此可偵測是否存在盆裝效應，甚至是批次間的效應。先去除具缺值的種原 (GN633、GN678、GN875) 後，對剩餘 13 個種原的冷處理組觀測值作變方分析 (表 6)。批次 (batch) 間的效應不顯著，盆裝 (box) 間的效應也不顯著，這說明 候選眾種原的冷處理表現，無論優劣都算穩定，受其它外在因素的干擾有限；
種原 (variety) 間的差異達顯著水準，P-value < 0.001，代表候選種原間確實存 在耐冷特性上的差異，挑選方具意義。
去除 GWAS 時未種植的種原 (GN693)，剩餘 15 個種原再以相關分析觀察 GWAS 與本次驗證的耐冷表現 (圖 12)，兩次試驗的相關係數為 0.6286，可發現
大多數的種原都很貼近斜率為1 的斜直線，即兩次試驗的表現十分接近，但有
Class Number Subpopulation Hap # Country
R* 9 TEJ Hap 2 > 40˚ N /S
圖 11 16 個候選種原的外表型驗證
x 軸所示編號 B1-01 代表第 1 批的第 1 盆裝，其餘以此類推，不同盆裝以虛線 相隔，不同批次以實線相隔；y 軸為各單株之 TRls。圖中以藍色表示DAS 15 時 移至15/13˚C (日/夜溫) 自然光室的冷處理組 (Trt)、紅色表示 DAS 15 時仍生長 於30/25˚C (日/夜溫) 自然光室的對照組 (Con)。
圖 12 外表型試驗與 GWAS 試驗中外表型的相關分析
左圖為15 個種原的結果；右圖為去除 GN619，14 個種原的結果。圖中紅點代 表R (Resistance) 的種原，藍點代表 S (Sensitive) 的種原；星型點則代表預備試 驗時也有種植的種原。
第三節 F
2族群對關聯性定位的驗證
1. 外表型分佈
F2族群的TRls分佈如圖13、圖 14。冷處理組與正常處理組相比整體的 TRls有所下降 (圖 13)，代表冷處理確實會抑制根系的生長；而不同盆裝間也沒 有效應導致盆裝間的TRls有表現上的差異。正常處理組有稍微右偏的情形 (圖 14 a)，常態性檢定 Shapiro-test 的結果 P-value 為 0.00259，無法認定其為常態分 佈，而兩親本Kon Suito (GN721) 及 TCS10 的平均值分別為 33.19 及 25.43，兩 親本的TRls表現並沒有特別的差距，不過整體的變異範圍較大。冷處理組的分 佈並無明顯偏歪 (圖 14 b)，常態性檢定 Shapiro-test 的結果 P-value 為 0.3798，
無法拒絕其為常態分佈，而兩親本Kon Suito 及 TCS10 的平均值分別為 16.13 及6.49，由圖中可見其在整個族群的分佈中比起正常處理組更為偏向兩端，但 整體的變異範圍卻不如正常處理組。先前試驗中親本的挑選指標為冷處理下的 TRls，表現優良的候選種原與不具耐冷特性的品種雜交，其後裔於冷處理下比 起正常處理更符合常態分佈，這樣的結果符合我們的設想。
另外我們也觀察到超越親本分離的現象，無論哪種處理下 Kon Suito 皆非表 現最好者，TCS10 也非最劣者，F2族群中仍有比親本表現更好或更差的後裔，
這也說明控制根系生長及耐冷特性都可能由超過一組的基因共同控制，並非由 單一基因所決定。
圖 13 F2族群於正常處理與冷處理下各單株的TRls
x 軸由左至右依盆號排列，兩種處理各 4 盆。以紅色點代表正常處理組，藍色 點代表冷處理組。
(a) (b)
圖 14 F2族群的TRls分佈
左圖為正常處理組，右圖為冷處理組；圖下方橫線為兩親本之平均±標準差，
左圖為正常處理組，右圖為冷處理組；圖下方橫線為兩親本之平均±標準差，