不同光周期條件下谷子株高的全基因組關(guān)聯(lián)分析

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

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

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

株高受遺傳網(wǎng)絡(luò)調(diào)控，目前發(fā)現(xiàn)的矮稈表型多是由赤霉素、油菜素生物合成或者信號(hào)識(shí)別途徑產(chǎn)生缺陷造成的[2-4]。此外，一些和細(xì)胞分裂、伸長相關(guān)的基因突變也會(huì)導(dǎo)致植株矮化，如水稻中編碼F-box域基因D3發(fā)生突變會(huì)導(dǎo)致矮化表型，編碼微管相關(guān)蛋白的dgl1基因缺陷也會(huì)導(dǎo)致植株矮化[5-6]。

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

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

1 材料和方法

1.1 試驗(yàn)材料及表型鑒定

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

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

1.2 谷子材料重測序

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

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

1.3 SNP標(biāo)記的開發(fā)

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

1.4 群體結(jié)構(gòu)分析

采用主成分分析法確定98份谷子材料的群體結(jié)構(gòu)，選擇均勻分布谷子9條染色體的3 600個(gè)SNP位點(diǎn)，用GCTA軟件進(jìn)行主成分分析，得到樣品的群體聚類情況。

1.5 連鎖不平衡(LD)分析

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

1.6 株高與SNP標(biāo)記的關(guān)聯(lián)分析及候選基因的篩選

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

2 結(jié)果與分析

2.1 株高性狀的表型分析

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

表2 不同種植區(qū)光周期條件下株高的表型變異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可以看出，海南種植區(qū)短日照條件下株高集中分布在100～120 cm，這個(gè)區(qū)間往兩側(cè)移動(dòng)呈現(xiàn)遞減趨勢，基本符合正態(tài)分布規(guī)律；洛陽種植區(qū)中日照條件下株高同樣集中分布在100～120 cm，這個(gè)區(qū)間往兩側(cè)移動(dòng)也表現(xiàn)出一定的遞減趨勢，總體上接近正態(tài)分布規(guī)律；吉林種植區(qū)長日照條件下株高分布表現(xiàn)一定雙峰模式，在120 cm附近有一個(gè)峰，另一個(gè)峰在140 cm附近。株高符合多基因控制的數(shù)量性狀遺傳模式。

2.2 98份谷子材料基因組重測序結(jié)果

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

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

表3 SNP類型統(tǒng)計(jì)Tab.3 SNP types statistics

2.3 群體結(jié)構(gòu)分析

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

2.4 LD分析

從圖3可以看出，當(dāng)r2值降低到最高值的一半時(shí)，SNP標(biāo)記間距離為47.5 kb，因此，本研究中谷子基因組LD衰減距離為47.5 kb。較小的LD衰退距離有利于縮小后續(xù)關(guān)聯(lián)分析候選區(qū)域，提高關(guān)聯(lián)結(jié)果的精確度。

2.5 谷子株高的全基因組關(guān)聯(lián)分析

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

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

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

在1號(hào)染色體與株高關(guān)聯(lián)的SNP位點(diǎn)兩側(cè)100 kb范圍內(nèi)搜索候選基因，在海南種植區(qū)搜索到10個(gè)候選基因，這些基因功能主要包括維生素E生物合成、細(xì)胞膜組成、代謝過程、DNA結(jié)合以及一些功能未知基因。在洛陽種植區(qū)區(qū)搜索到4個(gè)候選基因，其中2個(gè)基因功能為代謝、DNA結(jié)合，還有2個(gè)功能未知基因。在吉林種植區(qū)搜索到1個(gè)候選基因，功能為線粒體組成。在海南、洛陽種植區(qū)共同搜索到的候選基因有3個(gè)，包括1個(gè)功能未知基因、1個(gè)成蛋白基因(DRF)、一個(gè)DNA結(jié)合蛋白基因(ESCAROLA)，只有成蛋白基因外顯子區(qū)檢測到SNP位點(diǎn)(SNP14876527)，其余2個(gè)候選基因未發(fā)現(xiàn)SNP位點(diǎn)(表5)。因此，推測成蛋白基因突變可能導(dǎo)致了海南、洛陽種植區(qū)谷子株高的變異，該基因?yàn)橹旮咧饕蜻x基因。

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

表4 谷子材料1號(hào)染色體與株高關(guān)聯(lián)的SNPs位點(diǎn)Tab.4 The SNPs associated with plant height on Chr.1 of foxtail millet

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

3 結(jié)論與討論

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

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

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