基于重測序的短節間番茄品系變異檢測分析

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

(1.河北科技師范學院園藝科技學院，河北秦皇島 066004；2.河北省特色園藝種質挖掘與創新利用重點實驗室，河北秦皇島 066004； 3.河北省承德市農業農村局，河北承德 067000)

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

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

目前，關于番茄節間長度的研究，多集中于對長節間番茄的研究，而對短節間番茄突變體的研究鮮見報道。因此，本研究對2個節間長度不同的番茄品系進行全基因組重測序，并對其SNP、InDel、拷貝數變異(copy number variation，CNV)、SV等4種變異類型進行深度挖掘，以及基因變異對代謝通路的影響，從而為分子水平改良番茄節間長度的研究奠定基礎。

1 材料與方法

1.1 試驗材料

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

1.2 試驗方法

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

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

1.2.3 數據統計 采用BWA(0.7.12-r1039)mem程序將經過過濾后得到的高質量數據比對到參考基因組上，比對的參數均按照bwamem的默認參數。采用ANNOVAR軟件對SNP位點和InDel位點進行注釋。采用GATK(https：//www.broadinstitute.org/gatk/)(RealnTimes = 1)進行CNV檢測，使用BreakDancer進行結構變異(SV)檢測。

2 結果與分析

2.1 CH與DH農藝性狀分析

在開花結果期對2個番茄品系的株高、總節間長度、平均節間長度進行分析(表1、圖1)。結果表明，CH株高均極顯著高于DH；CH與DH相比總節間長度不存在顯著差異；CH與DH之間的平均節間長度均存在極顯著差異。

表1 CH與DH開花結果期農藝性狀分析

2.2 變異檢測分析

2.2.1 CH與DH基因組重測序結果與分析 由表2可見，CH和DH這2個樣本的Reads總數分別為 168 148 334、152 257 248個，過濾后Reads數量及過濾后Reads數量占原始Reads數量的百分比分別為160 322 604個(95.35%)、145 443 866個(95.53%)；堿基總數分別為25 222 250 100、22 838 587 200 bp，過濾后Reads堿基數及高質量Reads堿基數占原始堿基總數的百分比分別為 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的堿基數據統計

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

表3 CH與DH序列比對結果統計

表4 CH與DH比對測序深度和覆蓋度結果統計

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

表5 CH與DH SNP檢測結果統計

依據參考基因組比對結果，對CH和DH中的SNP進行比較(表6)，結果表明，CH和DH中共檢測到1 086 531個SNP。對2個樣本之間的全基因組SNP變異進行注釋，發現SNP變異主要集中在基因間區，其比例約占總數的78.03%；其次是發生于內含子區域內，其比例約占總數的8.48%；而發生于CDS區域的非同義突變比例約占CDS區域總數的60.85%。

表6 CH與DH SNP注釋結果

基因組SNP突變可以分成6類，分別為T：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檢測及注釋 為了定位目標性狀，每組性狀相關樣本一起call群體InDel，經過濾對最終得到的群體SNP在每個樣本中的數目做統計見表7。由表7可知，CH和DH的純合基因型數分別為15 327、259 512個，CH中純合基因型數小于DH；雜合基因型數分別為76 637、29 143個，CH中雜合基因型數多于DH；未知基因型數分別為3 180、733個；與參考基因組不一致的純合基因型數分別為265 384、71 140個。

表7 CH與DH InDel檢測結果統計

為鑒定CH與DH之間的InDel變異，分別進行比較(表8)。結果表明，CH和DH之間共檢測到277 973個InDel。對2個樣本之間的全基因組InDel變異進行注釋，發現InDel變異主要集中在基因間區域，其比例約占總數的76.53%；其次是發生于內含子區域內，其比例約占總數的11.67%；而發生于CDS區域的移碼突變比例約占總數的0.55%。

表8 CH與DH InDel注釋結果統計

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

表9 CH與DH CNV變異分析

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

表10 CH與DH SV變異分析

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

2.2.8 番茄短節間形成相關基因挖掘 為挖掘參與調控番茄短節間形成的相關基因，以SNP以及InDel中的基因作為切入點，進行變異基因挖掘。通過代謝通路富集情況，篩選出3個與節間長度相關聯的變異基因，各基因功能注釋、變異類型見表11，根據功能注釋獲得了控制赤霉素氧化酶、赤霉素受體、生長素響應因子的基因。

表11 變異基因功能注釋

2.2.9 候選基因表達分析 基于qRT-PCR 技術檢測開花結果期中與節間長度相關聯的變異基因，并進行相對表達量分析。由圖5可得，3個變異基因在2個試驗材料中均表現出差異。基因在DH中的表達量是CH中的2.85倍；基因在DH中的表達量是CH中的5.87倍；基因在DH中的表達量是CH中的1.58倍。在DH中3個基因的表達量明顯高于CH，以此推斷、、在番茄節間長度的調控中發揮著至關重要的作用。

3 討論與結論

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

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