不同光周期條件下谷子株高的全基因組關聯分析

賈小平,張 博，全建章，李劍峰，王永芳，袁璽壘

(1.河南科技大學 農學院，河南 洛陽 471023；2.河北省農林科學院 谷子研究所，國家谷子改良中心，河北 石家莊 050035)

株高是農作物的重要性狀，對作物的抗倒伏性、產量有重要影響。著名的“綠色革命”就是將半矮稈性狀引入水稻育種，選育出矮化耐密植的新品種，使產量得到提升[1]。明確株高的分子遺傳基礎有利于利用分子育種手段選育理想株型的新品種，使作物的生態適應能力、水分利用效率、機械收割能力、產量潛能等獲得綜合改善。

株高受遺傳網絡調控，目前發現的矮稈表型多是由赤霉素、油菜素生物合成或者信號識別途徑產生缺陷造成的[2-4]。此外，一些和細胞分裂、伸長相關的基因突變也會導致植株矮化，如水稻中編碼F-box域基因D3發生突變會導致矮化表型，編碼微管相關蛋白的dgl1基因缺陷也會導致植株矮化[5-6]。

谷子(Setariaitalica)作為我國的傳統作物，營養組成均衡，其較小的基因組(約510 Mb)、嚴格自花授粉、高光合等特點可以作為一個理想C4作物以及能源作物的模型[7-8]。株高與谷子倒數第一、第二、第三節倒伏指數存在一定負相關，而與穗長、穗粗、穗粒質量、穗碼數4個產量相關性狀呈顯著正相關[9]。近10年基于雙親雜交分離群體的數量性狀位點(QTL)定位方法已經用于谷子株高遺傳基礎的解析[10-13]。雖然通過上述研究定位了一些株高相關QTL位點，但是由于雙親的等位變異有限，檢測到的QTL位點也存在一定局限性，而且構建分離群體也較為耗時耗力。

全基因組關聯分析以自然群體為研究對象，利用分布于整個基因組的高密度DNA標記來發掘與目標性狀關聯的遺傳位點，相比傳統的基于雙親雜交分離群體的定位方法省去了構建分離群體的過程，而且可以檢測到更多的目標性狀關聯位點。目前全基因組關聯分析已經廣泛應用于植物數量性狀定位研究[14-17]。谷子基因組測序已經完成[18-19]，但是有關全基因組關聯分析的研究較少。Gupta等[20]調查了184份谷子材料的產量相關性狀，用50對簡單序列重復(SSR)標記進行了性狀與標記的關聯分析，發現有2個SSR標記與旗葉寬、籽粒產量有較強關聯，且位于5號染色體。Jia等[21]對916種不同的谷子資源進行重測序，鑒定了258萬個單核苷酸多態性(SNP)位點，并使用80萬個高質量的SNP標記來構建谷子基因組的單倍型圖譜，與47個性狀進行了關聯分析，獲得了多個地理環境穩定表達的QTL位點，這是首次報道的基于重測序的谷子全基因組關聯分析方面的研究。雖然目前已有較多關于谷子株高QTL定位的報道，但是針對不同光周期條件下谷子株高定位方面的研究還未開展。本研究在海南種植區短日照、洛陽種植區中日照、吉林種植區長日照3個不同光周期條件下，調查了98份來自國內以及國外的不同生態區谷子材料株高，并對這98份谷子材料進行重測序，利用測序數據開發海量SNP位點，開展位點與株高的全基因組關聯分析，發掘與株高性狀緊密連鎖的分子標記，并且對關聯位點區域存在的候選基因進行預測，以期為進一步克隆控制株高的基因奠定基礎。

1 材料和方法

1.1 試驗材料及表型鑒定

98份谷子材料包括了89份國內材料，來自河南、河北、山東、山西、陜西、甘肅、寧夏、內蒙古、吉林、遼寧、黑龍江、新疆、湖南、四川等地區；9份國外材料，來自美國、法國、日本、朝鮮、南非及國際半干旱研究所，材料信息見表1。

2015年5月中旬至10月中旬將98份谷子材料分別種植于河南科技大學開元校區試驗田和吉林市農業科學院試驗田，2015年11月初至2016年2月末將98份谷子材料種植于海南省樂東縣九所鎮河北農林科學院谷子研究所試驗田。每份谷子材料種植2行，行長2 m、行距50 cm、株距3～5 cm。按當地常規管理，每份材料于成熟期選取行中部谷子(共10株)測量株高，取平均值作為最終測量值。用SPSS 19.0軟件繪制直方圖及計算株高的廣義遺傳力。

1.2 谷子材料重測序

98份谷子材料長到三葉期時采集幼嫩葉片，利用十六烷基三甲基溴化銨(CTAB)法提取基因組DNA，提取的DNA滿足OD260/280=1.8～2.0，且沒有明顯降解，送上海美吉生物技術公司進行重測序。

表1 98份谷子材料信息Tab.1 Information of 98 foxtail millet varieties

1.3 SNP標記的開發

采用GATK軟件進行SNP、InDel(插入、缺失)信息的鑒定，并用Samtools提供的vcfutils工具對檢測到的SNP、InDel位點進行過濾，再利用ANNOVAR軟件以及參考基因組的基因組特征描述(gff)信息對SNP位點注釋。

1.4 群體結構分析

采用主成分分析法確定98份谷子材料的群體結構，選擇均勻分布谷子9條染色體的3 600個SNP位點，用GCTA軟件進行主成分分析，得到樣品的群體聚類情況。

1.5 連鎖不平衡(LD)分析

用連鎖不平衡系數(r2)衡量染色體上SNP位點間的連鎖不平衡水平，r2=0時則說明該群體處于連鎖平衡或者無LD，r2越大，SNP位點間距離越接近；一般r2大于或等于0.8，則視為兩位點屬于連鎖不平衡。用r2衰退到一半時對應的位點間的距離作為連鎖不平衡衰減距離，用VCFtools軟件來計算連鎖不平衡的衰減的數值，并以非線性回歸模型繪制出r2隨著SNP位點間物理距離的增加而衰減的趨勢圖。

1.6 株高與SNP標記的關聯分析及候選基因的篩選

以群體結構作為固定效應，個體親緣關系作為隨機效應，校正群體結構和個體間親緣關系的影響，采用Tassel Version 5軟件中的混合線性模型(Mixed Linear Model, MLM)進行性狀與標記間的關聯分析，當SNP位點的P<0.000 1時視該SNP位點與性狀關聯。把與株高關聯的SNP標記定位到谷子的基因組上，并確定相對應的置信區域，在關聯SNP標記連鎖不平衡衰減距離的區域范圍內，搜尋與性狀相關的候選基因。對篩選出的基因與基因本體聯合會(GO)、非冗余(NR)、京都基因與基因組百科全書(KEGG)等功能數據庫比對，實現對候選基因的注釋。

2 結果與分析

2.1 株高性狀的表型分析

在3個不同光周期條件下，谷子株高的表型變異和頻率分布分別見表2和圖1。谷子株高性狀均表現較大的變異范圍，其中海南種植區短日照條件下株高變異在58.5～152.4 cm，平均株高為108.1 cm；洛陽種植區中日照條件下株高變異在70.7～155.4 cm，平均株高為110.5 cm；吉林種植區長日照條件下變異在70.7～169.3 cm，平均株高127.6 cm。利用約束最大似然法估算出株高的廣義遺傳力為0.501。隨著日照時間的增長，谷子的株高也呈現增高趨勢，海南、洛陽2個種植區光周期間株高差異不顯著，但是海南、洛陽和吉林3個種植區光周期間株高存在顯著差異(P<0.01)。

表2 不同種植區光周期條件下株高的表型變異Tab.2 Phenotypic variation of plant height under photoperiod conditions of different cultivation regions

注：不同大寫字母代表極顯著(P<0.01)，不同小寫字母代表顯著(P<0.05)，相同小寫字母代表不顯著(P>0.05)。

Note: The different capital letters represent extremely significant (P<0.01), the different lower-case letters represent significant (P<0.05), the same lower-case letter represents no significant (P>0.05).

從圖1可以看出，海南種植區短日照條件下株高集中分布在100～120 cm，這個區間往兩側移動呈現遞減趨勢，基本符合正態分布規律；洛陽種植區中日照條件下株高同樣集中分布在100～120 cm，這個區間往兩側移動也表現出一定的遞減趨勢，總體上接近正態分布規律；吉林種植區長日照條件下株高分布表現一定雙峰模式，在120 cm附近有一個峰，另一個峰在140 cm附近。株高符合多基因控制的數量性狀遺傳模式。

2.2 98份谷子材料基因組重測序結果

98份谷子材料重測序獲得序列與豫谷1號參照基因組比對，基因組覆蓋率在86.10%～97.55%，測序深度在6.407 9～13.327 5倍，平均測序深度為6.8倍，表明基因組被覆蓋地較均勻，測序的隨機性較好。利用基因組重測序檢測到900多萬SNP位點，經過濾獲得4 482 208個高質量的SNP標記，這些SNP標記有4 037 579個在基因間隔區，占所有標記總數的90%。剩余的SNP標記主要分布在內含子區域、轉錄起始位點的上游區域、基因下游區域，其相應數目為122 698，102 355，81 070。在外顯子區域內的標記可以分為5種情況：非同義變異、終止密碼子獲得變異、終止密碼子丟失變異、同義變異、未知變異，其數目分別為31 477，459，66，25 381，332個，共計57 715個。在非編碼RNA區域內的標記也可根據位置可分為5種：非編碼RNA外顯子區域、非編碼RNA內含子區域、非編碼RNA 3′UTR、非編碼RNA 5′UTR、非編碼RNA剪接位點區域，其具體的SNP標記數目分別為14 371，28 855，39，21，101，共計為43 387(表3)。

圖1 株高在不同種植區光周期條件下的頻率分布Fig.1 Plant height frequency under photoperiod conditions of different cultivation regions

表3 SNP類型統計Tab.3 SNP types statistics

2.3 群體結構分析

利用所開發的3 000個SNP標記，通過 GCTA 軟件進行主成分分析，得到98個谷子材料的主要成分聚類情況。如圖2所示，98個谷子材料聚為三維，pca1、pca2 、pca3 分別為3個主成分，說明最佳分群數為3個。

2.4 LD分析

從圖3可以看出，當r2值降低到最高值的一半時，SNP標記間距離為47.5 kb，因此，本研究中谷子基因組LD衰減距離為47.5 kb。較小的LD衰退距離有利于縮小后續關聯分析候選區域，提高關聯結果的精確度。

2.5 谷子株高的全基因組關聯分析

將重測序開發的SNP標記與3個不同光周期條件下的株高數據進行關聯分析，共檢測到10 703個SNP位點與株高性狀顯著關聯(P<0.000 1)，其中海南種植區短日照條件與株高關聯的SNP位點數量最多，為10 608個，且多數集中在1號染色體上(10 458個)，其次是洛陽種植區中日照條件，關聯的SNP位點數量為78個，吉林種植區長日照條件關聯的SNP位點數量最少，為17個(圖4)。在海南種植區短日照條件下關聯的SNP位點在谷子的9條染色體上均有分布，在洛陽種植區中日照條件下關聯的SNP位點分布在除3號、8號以外的其他染色體上，在吉林種植區長日照條件下關聯的SNP位點分布在除2號、9號以外的其他染色體上。其中1號染色體上的3個標記(SNP13861443、SNP14872616、SNP18601830)在海南、洛陽2個種植區光周期條件下同時檢測到，而在吉林種植區長日照條件下未檢測到，這3個株高關聯SNP位點分別位于1號染色體的13 861 443，14 872 616，18 601 830 bp處。其余SNP標記不能在2個或者3個光周期條件下同時檢測到。因此，推定1號染色體可能存在中、短日照條件表達的控制株高的QTL位點，而在吉林種植區長日照條件下也在1號染色體檢測到3個與株高關聯的SNP位點(SNP8797106、SNP8797129、SNP9880380)(表4)，說明谷子1號染色體可能存在對光周期具有不同反應的控制株高的QTL位點。

圖2 98個谷子材料PCA聚類Fig.2 PCA cluster diagram of 98 foxtail millet varieties

圖3 LD衰減曲線Fig.3 The attenuation curve of LD

在1號染色體與株高關聯的SNP位點兩側100 kb范圍內搜索候選基因，在海南種植區搜索到10個候選基因，這些基因功能主要包括維生素E生物合成、細胞膜組成、代謝過程、DNA結合以及一些功能未知基因。在洛陽種植區區搜索到4個候選基因，其中2個基因功能為代謝、DNA結合，還有2個功能未知基因。在吉林種植區搜索到1個候選基因，功能為線粒體組成。在海南、洛陽種植區共同搜索到的候選基因有3個，包括1個功能未知基因、1個成蛋白基因(DRF)、一個DNA結合蛋白基因(ESCAROLA)，只有成蛋白基因外顯子區檢測到SNP位點(SNP14876527)，其余2個候選基因未發現SNP位點(表5)。因此，推測成蛋白基因突變可能導致了海南、洛陽種植區谷子株高的變異，該基因為株高主要候選基因。

圖4 海南、洛陽、吉林3個種植區株高的曼哈頓散點圖Fig.4 Manhattan plot of plant height in Hainan, Luoyang and Jilin cultivation regions

表4 谷子材料1號染色體與株高關聯的SNPs位點Tab.4 The SNPs associated with plant height on Chr.1 of foxtail millet

表5 海南、洛陽種植區共同檢測到的株高關聯候選基因Tab.5 The candidate genes associated with plant height in Hainan and Luoyang cultivation regions

3 結論與討論

本研究在海南種植區短日照、洛陽種植區中日照2個光周期條件下獲得3個與谷子株高關聯的SNP位點，這3個位點都位于1號染色體上。王曉宇等[10]在谷子2號染色體定位到2個控制株高的QTL位點；趙美丞[11]利用圖位克隆法在谷子9號染色體獲得一個半顯性矮稈基因SiDw1；Wang等[12]、Zhang等[13]利用高密度的SNP標記將谷子株高控制位點定位于1、4、5、6、7、9號染色體上。前人在1號染色體上將株高QTL位點定位在184.84，144.32 cM 2個峰值點附近區域[12-13]，本研究在1號染色體獲得的3個株高關聯SNP位點位于13 861 443，14 872 616，18 601 830 bp，盡管本研究關聯位點位置使用的是物理距離，而前人報道的關聯位點使用的是遺傳距離，Wang等[12]報道的控制株高的QTL位點位于1號染色體的末端，本研究3個關聯位點位于1號染色體的中部靠前，初步判斷本研究獲得的QTL位點與前人報道的不在同一個區間內。Zhang等[13]在宣化種植區長日照條件下檢測到1號染色體144.32 cM位置存在株高QTL位點，而在海南三亞種植區短日照條件下并未在1號染色體檢測到控制株高的QTL位點，且長日照檢測到的QTL位點位于1號染色體的后部，據此推斷本研究在1號染色體關聯到的株高控制位點應該為在短日照、中日照條件特定表達的新的QTL位點。

本研究在海南、洛陽種植區2個光周期條件關聯到3個候選基因，但是2個候選基因在基因內部未檢測到SNP位點，只有一個成蛋白類基因(LOC101783280)在外顯子區檢測到一個SNP位點(SNP14876527)，推測該基因突變可能和株高相關。目前，在水稻等作物報道的株高相關基因主要包括參與赤霉素生物合成及信號傳遞、油菜內酯生物合成、獨角金內酯生物合成3類基因[22]。成蛋白或者形成素(Formin)是一類多結構域蛋白，具有保守的串聯FH1、FH2結構域，對微絲的聚合有重要作用。植物中的成蛋白分為2個亞家族，第一亞家族在氨基端具有跨膜域和信號肽，第二亞家族具有PTEN結構域。植物成蛋白基因也具有微絲成核、成束、聚合等作用，在細胞分裂、極性生長中具有重要作用[23-24]。已有研究表明，水稻成蛋白基因OsFH5突變會導致植株矮化、頂部節尖彎曲等表型[25]，本研究在海南、洛陽2個種植區光周期條件下都檢測到與株高關聯的成蛋白類基因，該基因突變是否是導致谷子植株矮化的主要原因仍需進一步做基因表達分析及轉基因驗證。

在海南、洛陽、吉林3個種植區不同光周期條件下調查了98份谷子材料的株高性狀，3個種植區株高變異幅度在58.5～169.3 cm，廣義遺傳力為0.501。對98份谷子材料進行重測序開展SNP標記與株高的全基因組關聯分析，獲得10 703個與株高關聯的SNP位點，其中只有1號染色體的3個SNP(SNP13861443、SNP14872616、SNP18601830)在海南、洛陽2個種植區光周期條件下能被穩定檢測到，在3個關聯SNP候選區域檢測到3個候選基因，只有推定的成蛋白基因(LOC101783280)在外顯子區檢測到一個SNP位點(SNP14876527)，因此，該成蛋白基因可能是海南短日照、洛陽中日照光周期條件下導致谷子株高變異的主要候選基因。