基于重測(cè)序的短節(jié)間番茄品系變異檢測(cè)分析

李毅豐， 姜悅暢， 孫中鋒， 王 帥，2， 唐貝貝， 師海林， 張 寧，2， 王玉斌， 毛秀杰，2

(1.河北科技師范學(xué)院園藝科技學(xué)院，河北秦皇島 066004；2.河北省特色園藝種質(zhì)挖掘與創(chuàng)新利用重點(diǎn)實(shí)驗(yàn)室，河北秦皇島 066004； 3.河北省承德市農(nóng)業(yè)農(nóng)村局，河北承德 067000)

番茄(L.)是重要的蔬菜作物和模式植物。實(shí)際生產(chǎn)應(yīng)用中，集約化的育苗方式增加了番茄育苗密度，植株間相互遮陰，促使幼苗發(fā)生徒長(zhǎng)，因此番茄節(jié)間長(zhǎng)度的研究受到廣泛關(guān)注。番茄節(jié)間長(zhǎng)度是影響番茄產(chǎn)量的重要農(nóng)藝性狀之一，不僅會(huì)對(duì)番茄栽培密度產(chǎn)生影響，還對(duì)番茄的葉面積指數(shù)、光能利用率、干物質(zhì)積累以及果實(shí)品質(zhì)產(chǎn)生影響。番茄節(jié)間長(zhǎng)度縮短不僅能夠有效降低植株高度，抑制徒長(zhǎng)，增加土地利用率，減少土壤中營(yíng)養(yǎng)物質(zhì)的消耗，還能降低苗期管理和果實(shí)采收的人工成本，因此，對(duì)短節(jié)間番茄進(jìn)行研究具有重要意義。

隨著測(cè)序技術(shù)的迅速發(fā)展，通過(guò)高通量測(cè)序技術(shù)能夠?qū)υS多作物基因進(jìn)行研究。周世奇等對(duì)選育的航天突變體煙草NC89-M與野生型NC89進(jìn)行全基因組重測(cè)序，在NC89中檢測(cè)獲得了1 848 013個(gè)單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)、398 922個(gè)小片段插入和缺失(insertion-deletion，InDel)、41 969個(gè)染色體結(jié)構(gòu)變異(structural variation，SV)，在NC89-M中檢測(cè)到 1 876 219 個(gè)SNP，402 011個(gè)InDel，42 699個(gè)SV，采用生物信息學(xué)分析方法得出代謝通路和次生代謝產(chǎn)物合成相關(guān)基因突變數(shù)目最多，通過(guò)變異基因功能注釋獲得-基因和基因，分別調(diào)控開(kāi)花時(shí)間和側(cè)生器官發(fā)育與葉緣形狀。杜海東等對(duì)3個(gè)自封頂番茄品系進(jìn)行基因組重測(cè)序，通過(guò)變異檢測(cè)分析3個(gè)樣本檢測(cè)到 5 968 501 個(gè)SNP和485 114個(gè)InDel，與參考基因組比對(duì)后共發(fā)生33 473個(gè)變異基因。對(duì)CDS區(qū)域的變異基因進(jìn)行KEGG數(shù)據(jù)庫(kù)比對(duì)，發(fā)現(xiàn)主要集中在基礎(chǔ)代謝和玉米素的生物合成，通過(guò)多序列比對(duì)獲得了16個(gè)調(diào)控封頂花序數(shù)的關(guān)鍵基因。

目前，關(guān)于番茄節(jié)間長(zhǎng)度的研究，多集中于對(duì)長(zhǎng)節(jié)間番茄的研究，而對(duì)短節(jié)間番茄突變體的研究鮮見(jiàn)報(bào)道。因此，本研究對(duì)2個(gè)節(jié)間長(zhǎng)度不同的番茄品系進(jìn)行全基因組重測(cè)序，并對(duì)其SNP、InDel、拷貝數(shù)變異(copy number variation，CNV)、SV等4種變異類(lèi)型進(jìn)行深度挖掘，以及基因變異對(duì)代謝通路的影響，從而為分子水平改良番茄節(jié)間長(zhǎng)度的研究奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料

本試驗(yàn)以短節(jié)間番茄品系CH和普通節(jié)間番茄品系DH為試驗(yàn)材料，均為無(wú)限生長(zhǎng)型。試驗(yàn)在河北科技師范學(xué)院園藝園林試驗(yàn)站7號(hào)溫室中進(jìn)行，于2021年1月16日進(jìn)行播種育苗，待幼苗長(zhǎng)至5～6張真葉時(shí)進(jìn)行移栽。采用日光溫室栽培，大行距為70 cm，小行距為40 cm，株距為30 cm，應(yīng)用滴灌方式進(jìn)行水肥管理。在開(kāi)花結(jié)果期取幼嫩組織用于基因組重測(cè)序分析。

1.2 試驗(yàn)方法

1.2.1 全基因組DNA提取 通過(guò)十六烷基三甲基溴化銨(CTAB)法對(duì)樣本莖段進(jìn)行DNA的提??；用0.8%瓊脂糖凝膠電泳對(duì)DNA質(zhì)量進(jìn)行檢測(cè)；采用核酸蛋白分析儀對(duì)提取的DNA進(jìn)行定量。

1.2.2 全基因組重測(cè)序 對(duì)質(zhì)檢合格的DNA采用超聲進(jìn)行隨機(jī)打斷。采用TruSeq DNA PCR-Free Prep kit的標(biāo)準(zhǔn)建庫(kù)流程制備測(cè)序文庫(kù)，對(duì)序列末端進(jìn)行修復(fù)、3′端添加多聚腺苷酸(PloyA)、5′端添加含有文庫(kù)特異性標(biāo)簽、純化、測(cè)序文庫(kù)模板富集等步驟完成測(cè)序文庫(kù)的制備。對(duì)完成的測(cè)序文庫(kù)在Agilent Bioanalyzer上進(jìn)行質(zhì)檢；采用Quant-iT PicoGreen dsDNA Assay Kit對(duì)文庫(kù)進(jìn)行定量(文庫(kù)濃度大于2 nmol質(zhì)量合格)；根據(jù)所需測(cè)序量對(duì)合格的各上機(jī)測(cè)序文庫(kù)梯度稀釋后按相應(yīng)比例混合使用；應(yīng)用NovaSeq測(cè)序儀進(jìn)行雙端測(cè)序。

1.2.3 數(shù)據(jù)統(tǒng)計(jì) 采用BWA(0.7.12-r1039)mem程序?qū)⒔?jīng)過(guò)過(guò)濾后得到的高質(zhì)量數(shù)據(jù)比對(duì)到參考基因組上，比對(duì)的參數(shù)均按照bwamem的默認(rèn)參數(shù)。采用ANNOVAR軟件對(duì)SNP位點(diǎn)和InDel位點(diǎn)進(jìn)行注釋。采用GATK(https：//www.broadinstitute.org/gatk/)(RealnTimes = 1)進(jìn)行CNV檢測(cè)，使用BreakDancer進(jìn)行結(jié)構(gòu)變異(SV)檢測(cè)。

2 結(jié)果與分析

2.1 CH與DH農(nóng)藝性狀分析

在開(kāi)花結(jié)果期對(duì)2個(gè)番茄品系的株高、總節(jié)間長(zhǎng)度、平均節(jié)間長(zhǎng)度進(jìn)行分析(表1、圖1)。結(jié)果表明，CH株高均極顯著高于DH；CH與DH相比總節(jié)間長(zhǎng)度不存在顯著差異；CH與DH之間的平均節(jié)間長(zhǎng)度均存在極顯著差異。

表1 CH與DH開(kāi)花結(jié)果期農(nóng)藝性狀分析

2.2 變異檢測(cè)分析

2.2.1 CH與DH基因組重測(cè)序結(jié)果與分析 由表2可見(jiàn)，CH和DH這2個(gè)樣本的Reads總數(shù)分別為 168 148 334、152 257 248個(gè)，過(guò)濾后Reads數(shù)量及過(guò)濾后Reads數(shù)量占原始Reads數(shù)量的百分比分別為160 322 604個(gè)(95.35%)、145 443 866個(gè)(95.53%)；堿基總數(shù)分別為25 222 250 100、22 838 587 200 bp，過(guò)濾后Reads堿基數(shù)及高質(zhì)量Reads堿基數(shù)占原始?jí)A基總數(shù)的百分比分別為 23 522 020 839 bp(93.26%)、21 329 453 341 bp(93.39%)。GC含量在36.83～36.98之間；Q20≥96.82%、Q30≥91.85%。

表2 CH與DH的堿基數(shù)據(jù)統(tǒng)計(jì)

2.2.2 CH和DH與參考基因組對(duì)比情況 以Heinz 1706為參考基因組。2個(gè)樣本的比對(duì)率在99.61%～99.89%之間(表3)。平均測(cè)序深度為23X，1X覆蓋度在96.63%～99.74%之間，4X 覆蓋度在95.62%～99.52%之間，10X 覆蓋度在93.19%～96.97%之間，20X 覆蓋度在46.47%～52.31%之間(表4)。綜上，比對(duì)結(jié)果正常，可用于后續(xù)相關(guān)分析。

表3 CH與DH序列比對(duì)結(jié)果統(tǒng)計(jì)

表4 CH與DH比對(duì)測(cè)序深度和覆蓋度結(jié)果統(tǒng)計(jì)

2.2.3 SNP檢測(cè)與注釋 依據(jù)與參考基因組的比對(duì)結(jié)果(表5)可知，CH和DH的總SNP變異個(gè)數(shù)相同，為3 893 609個(gè)，CH和DH的純合基因型數(shù)分別為186 147、3 506 620個(gè)，CH中純合基因型低于DH；雜合基因型數(shù)分別為1 106 264、251 093個(gè)；未知基因型數(shù)分別為36 007、2 994個(gè)；與參考基因組不一致的純合基因型數(shù)分別為 2 565 191、132 902個(gè)。由SNP在染色體上的分布(圖2)可知，SNP在9號(hào)染色體的24.89～27.79 Mb、35.32～38.21 Mb區(qū)域內(nèi)富集，初步判斷在此區(qū)域內(nèi)發(fā)生基因突變，其余SNP均勻分布在6號(hào)染色體的1.74～8.11 Mb、12.16～31.26 Mb和9號(hào)染色體的7.53～18.53 Mb、21.42～62.53 Mb之間。

表5 CH與DH SNP檢測(cè)結(jié)果統(tǒng)計(jì)

依據(jù)參考基因組比對(duì)結(jié)果，對(duì)CH和DH中的SNP進(jìn)行比較(表6)，結(jié)果表明，CH和DH中共檢測(cè)到1 086 531個(gè)SNP。對(duì)2個(gè)樣本之間的全基因組SNP變異進(jìn)行注釋?zhuān)l(fā)現(xiàn)SNP變異主要集中在基因間區(qū)，其比例約占總數(shù)的78.03%；其次是發(fā)生于內(nèi)含子區(qū)域內(nèi)，其比例約占總數(shù)的8.48%；而發(fā)生于CDS區(qū)域的非同義突變比例約占CDS區(qū)域總數(shù)的60.85%。

表6 CH與DH SNP注釋結(jié)果

基因組SNP突變可以分成6類(lèi)，分別為T(mén)：A>C：G、T：AA：T、C：G>G：C、C：G>T：A。SNP頻譜分析如圖3所示，由圖3可知，T：A>C：G和C：G>T：A為主要的SNP突變型。

2.2.4 InDel檢測(cè)及注釋 為了定位目標(biāo)性狀，每組性狀相關(guān)樣本一起call群體InDel，經(jīng)過(guò)濾對(duì)最終得到的群體SNP在每個(gè)樣本中的數(shù)目做統(tǒng)計(jì)見(jiàn)表7。由表7可知，CH和DH的純合基因型數(shù)分別為15 327、259 512個(gè)，CH中純合基因型數(shù)小于DH；雜合基因型數(shù)分別為76 637、29 143個(gè)，CH中雜合基因型數(shù)多于DH；未知基因型數(shù)分別為3 180、733個(gè)；與參考基因組不一致的純合基因型數(shù)分別為265 384、71 140個(gè)。

表7 CH與DH InDel檢測(cè)結(jié)果統(tǒng)計(jì)

為鑒定CH與DH之間的InDel變異，分別進(jìn)行比較(表8)。結(jié)果表明，CH和DH之間共檢測(cè)到277 973個(gè)InDel。對(duì)2個(gè)樣本之間的全基因組InDel變異進(jìn)行注釋?zhuān)l(fā)現(xiàn)InDel變異主要集中在基因間區(qū)域，其比例約占總數(shù)的76.53%；其次是發(fā)生于內(nèi)含子區(qū)域內(nèi)，其比例約占總數(shù)的11.67%；而發(fā)生于CDS區(qū)域的移碼突變比例約占總數(shù)的0.55%。

表8 CH與DH InDel注釋結(jié)果統(tǒng)計(jì)

2.2.5 CH與DH CNV變異分析 采用CNVnator v0.2.7檢測(cè)3個(gè)樣本的全基因組中存在的CNV，通過(guò)CNV統(tǒng)計(jì)(表9)可知，CH中檢測(cè)出的CNV為 12 756 個(gè)，拷貝數(shù)缺失10 970(85.9%)，拷貝數(shù)增加1 786(14.1%)；在DH中檢測(cè)出的CNV為 13 041 個(gè)，拷貝數(shù)缺失11 576(88.8%)，拷貝數(shù)增加1 465(11.2%)。因此得出，CH和DH的CNV變異主要是發(fā)生拷貝數(shù)缺失。

表9 CH與DH CNV變異分析

2.2.6 CH與DH SV變異分析 染色體的結(jié)構(gòu)變異主要包括缺失(deletion，DEL)、插入(insertion，INS)、倒位(inversion，INV)、染色體內(nèi)易位(intra-chromosomal translocation，ITX)、染色體間易位(inter-chromosomal translocation，CTX)。采用Breakdancer1.3.7版本檢測(cè)染色體結(jié)構(gòu)變異，5種結(jié)構(gòu)變異的數(shù)量統(tǒng)計(jì)見(jiàn)表10。CH的SV數(shù)量最多為20 501個(gè)，其中染色體缺失的數(shù)量最多為8 176個(gè)，占總變異的39.9%，染色體間易位和染色體內(nèi)易位分別占總變異的24.6%、20.0%，插入和倒位分別占總變異的5.5%、10.0%。DH的SV數(shù)量最少為13 015個(gè)，其中變異最多的2種類(lèi)型為染色體間易位和缺失分別為4 668個(gè)、4 468個(gè)，占總變異的35.9%、34.3%，插入、倒位、染色體內(nèi)易位分別占總變異的3.7%、6.5%、19.6%。CH和DH的SV變異主要發(fā)生的變異類(lèi)型為缺失。

表10 CH與DH SV變異分析

2.2.7 CH與DH代謝水平的差異分析 為了解基因變異引起代謝物質(zhì)產(chǎn)生差異，因此對(duì)CH與DH節(jié)間中的代謝物質(zhì)進(jìn)行了分析。根據(jù)差異代謝物結(jié)果，對(duì)CH與DH的KEGG通路進(jìn)行比較，獲得KEGG通路富集圖(圖4)。由圖4可知，在CH和DH中共注釋到20個(gè)代謝通路，其中有3個(gè)通路被顯著富集，包含吲哚生物堿生物合成、二萜生物合成、草莽酸途徑生物堿的生物合成。

2.2.8 番茄短節(jié)間形成相關(guān)基因挖掘 為挖掘參與調(diào)控番茄短節(jié)間形成的相關(guān)基因，以SNP以及InDel中的基因作為切入點(diǎn)，進(jìn)行變異基因挖掘。通過(guò)代謝通路富集情況，篩選出3個(gè)與節(jié)間長(zhǎng)度相關(guān)聯(lián)的變異基因，各基因功能注釋、變異類(lèi)型見(jiàn)表11，根據(jù)功能注釋獲得了控制赤霉素氧化酶、赤霉素受體、生長(zhǎng)素響應(yīng)因子的基因。

表11 變異基因功能注釋

2.2.9 候選基因表達(dá)分析 基于qRT-PCR 技術(shù)檢測(cè)開(kāi)花結(jié)果期中與節(jié)間長(zhǎng)度相關(guān)聯(lián)的變異基因，并進(jìn)行相對(duì)表達(dá)量分析。由圖5可得，3個(gè)變異基因在2個(gè)試驗(yàn)材料中均表現(xiàn)出差異?；蛟贒H中的表達(dá)量是CH中的2.85倍；基因在DH中的表達(dá)量是CH中的5.87倍；基因在DH中的表達(dá)量是CH中的1.58倍。在DH中3個(gè)基因的表達(dá)量明顯高于CH，以此推斷、、在番茄節(jié)間長(zhǎng)度的調(diào)控中發(fā)揮著至關(guān)重要的作用。

3 討論與結(jié)論

番茄是重要的蔬菜作物，其節(jié)間長(zhǎng)度受多種因素影響，包括環(huán)境因子、植物激素、遺傳因子等。目前已有植物節(jié)間長(zhǎng)度分子研究的報(bào)道，劉根忠等通過(guò)對(duì)226份番茄核心種質(zhì)資源進(jìn)行全基因組關(guān)聯(lián)分析，鑒定出7個(gè)與節(jié)間長(zhǎng)度相關(guān)的數(shù)量性狀基因座(QTL)，通過(guò)對(duì)顯著的SNP位點(diǎn)進(jìn)行候選基因分析，挖掘出控制番茄節(jié)間長(zhǎng)度的主效基因和；劉忠祥等通過(guò)對(duì)玉米進(jìn)行外源赤霉素(GA)的噴施，確定玉米對(duì)GA的響應(yīng)正常，在利用分子標(biāo)記技術(shù)將控制株高的主效QTL qPH3.2共定位在第3染色體上。研究發(fā)現(xiàn)，不同的土壤環(huán)境會(huì)對(duì)番茄內(nèi)源激素含量產(chǎn)生影響，從而對(duì)植株的生長(zhǎng)發(fā)育產(chǎn)生影響；控制節(jié)間發(fā)育相關(guān)的基因通過(guò)影響代謝通路來(lái)控制植物節(jié)間發(fā)育。前人研究發(fā)現(xiàn)，分子標(biāo)記輔助育種能夠提高植物優(yōu)良性狀的選擇效率，對(duì)植物遺傳育種具有重要意義。隨著番茄全基因組測(cè)序的完成以及二代測(cè)序的迅速發(fā)展，使番茄基因組的變異檢測(cè)分析成為可能。

本研究通過(guò)對(duì)短節(jié)間番茄品系CH和普通節(jié)間番茄品系DH進(jìn)行全基因組重測(cè)序，與參考基因組Heinz 1706番茄相比，在DH中檢測(cè)到 3 893 609 個(gè)SNP、360 528個(gè)InDel、13 041個(gè)CNV、13 015個(gè)SV；在CH中檢測(cè)到3 893 609個(gè)SNP、360 528個(gè)InDel、12 756個(gè)CNV、20 501個(gè)SV。結(jié)果表明，2個(gè)番茄材料主要發(fā)生的變異類(lèi)型為單核苷酸多態(tài)性(SNP)和小片段插入和缺失(InDel)，與油菜、煙草和水稻的研究結(jié)果相吻合。對(duì)2個(gè)樣本的KEGG通路進(jìn)行比較，富集到的代謝通路包括吲哚生物堿生物合成、二萜生物合成、莽草酸途徑生物堿的生物合成。為后續(xù)對(duì)變異基因分子標(biāo)記的開(kāi)發(fā)、分子標(biāo)記輔助育種、基因克隆以及基因功能驗(yàn)證的研究奠定重要基礎(chǔ)。