新疆核桃‘新露’與‘新萃豐’全基因組重測序分析

武鵬雨，李冬，吳翠云，包建平，金強，虎海防 ，張銳*

（1塔里木盆地生物資源保護利用兵團重點實驗室，新疆 阿拉爾 843300）

（2南疆特色果樹高效優質栽培與深加工技術國家地方聯合工程實驗室，新疆 阿拉爾 843300）

(3塔里木大學園藝與林學學院，新疆 阿拉爾 843300)

（4新疆佳木果樹學國家長期科研基地，新疆 溫宿 843100）

核桃（Juglans regia L.）隸屬于胡桃科（Juglandaceae）、核桃屬（Juglans），是世界4大堅果（核桃、扁桃、板栗、腰果）之一，是不可或缺的堅果和木本油料樹種［1］。我國核桃栽培歷史悠久，是核桃屬（Juglans）植物的起源和分布中心之一［2］，新疆作為我國核桃重要產區之一，核桃栽培面積達3.5×105hm2，其中塔里木盆地周圍綠洲產量最高，位居全國前列［3］。新疆栽培核桃資源豐富，擁有野核桃和紙皮核桃等種質資源，也有薄皮核桃如‘扎343’‘新新2號’‘新萃豐’等品種，殼薄如紙，透光可見殼內果仁形狀，薄皮核桃硬殼的平均厚度一般小于1.5 mm［4］，在薄皮不露仁栽培種中殼厚最厚的是‘新萃豐’，厚度為1.45 mm［5］；而‘新露’是天然內果皮發育不全的露仁核桃品種，其內果皮的殼面在生長發育過程中發生自動破裂使核桃仁露出，果仁裸露雖方便食用，但不耐運輸儲藏。對擁有早實、薄殼、抗逆性強的‘新萃豐’和天然果殼裸露的‘新露’進行研究，可以使新疆核桃內果皮發育的遺傳背景清晰化，為新疆核桃的保護與科學利用提供理論依據。

隨著測序成本的降低和已知基因組序列物種的增多，基因組測序已成為研究植物分子育種、群體進化中最為迅速有效的方法之一［6］。基因組測序技術為不同物種提供了高質量的參考基因組序列，基因組重測序技術根據供參考的高質量序列進行整個基因組、區域或基因的重新測序，每個個體的測序數據與高質量的參考數據進行比對，不同樣本之間的遺傳變異可以被檢測為高度可信的序列差異，并在此基礎上對個體或群體進行比較基因組學分析，可以從基因組水平對重要性狀的候選功能基因或QTL進行定位與分析［7］，通過與參考基因組比較獲得多態性插入與缺失（InDel）標記，使用這些分析標記用于基因鑒定［8］，構建雜交群體使用基因組重測序單核苷酸多態性位點（SNP）數據庫篩選與目的形狀相關基因［9］。本研究通過Illumina NovaSeq 6000測序平臺對新疆栽培核桃中的‘新露’和‘新萃豐’進行全基因組重測序分析，對其SNP和InDel等變異位點進行深度挖掘，從基因組層面分析新疆核桃內果皮發育過程中出現的果仁裸露這一性狀的機理，通過檢測變異位點所在區域，為后續內果皮發育等性狀的相關基因的挖掘、品種選育等研究提供一定的基礎。

1 材料與方法

1.1 試驗材料

試驗所選‘新露’與‘新萃豐’品種材料種植于新疆佳木果樹學國家長期科研基地，2019年選取生長情況一致的4年生核桃植株各4株，在5月份新出嫩葉時采摘葉片，放入液氮中保存帶回實驗室。

1.2 試驗方法

‘新露’與‘新萃豐’的DNA提取方法按照DNA提取試劑盒的使用說明操作（TianGen）。DNA樣品的質量控制通過0.8%瓊脂糖凝膠電泳檢測，同時采用紫外分光光度計對DNA進行定量分析。

1.2.1 測序

對提取的DNA進行超聲隨機打斷，隨后進行末端損傷修復及連接接頭，在與測序接頭連接后進行PCR擴增，對測序文庫模板進行富集并純化文庫產物，連接產物通過2%瓊脂糖凝膠電泳進行最終片段選擇與純化，選擇插入片段大約為400 bp回收后根據制造商建庫流程構建測序文庫。在Illumina Nova-Seq 6000測序平臺上對構建的文庫進行雙端（PE 150 bp）高通量測序。

1.2.2 數據分析

獲得原始數據后采用FastQC將數據進行質量控制，進而獲得高質量數據。使用Burrows-Wheeler-Alignment Tool軟件將高質量數據比對到已測序完成的核桃參考基因組，參考基因組數據為：JrSerr v1.0(GCA_004785585.1)［10］，利用 picard 軟件將比對得到的sam結果文件轉化為bam文件，使用GATK軟件對SNP和InDel進行檢測，并用Unified Genotyper對變異位點位置進行獲取，使用ANNOVAR軟件對獲取的SNP和InDel進行功能注釋。獲取‘新露’與‘新萃豐’相同位置的SNP及InDel后，再根據非同義SNP（nonsynonymous，nsSNP）/InDel獲取相關基因。候選基因的富集分析遞交于AgriGO軟件用于富集Gene ontology terms。

2 結果與分析

2.1 測序數據預處理

以Illumina Nova Seq測序平臺提供的初始測序數據為原始數據（Raw Data），即本次測序中獲得的各樣本得到的短序列數（Reads）和堿基總數（Bases），對2個核桃樣本產出的數據進行匯總（見表1）。2個品種共得到232 715 892個短序列，其中‘新露’核桃短序列數目127 572 756個，Q30比為91.30%；‘新萃豐’核桃中Reads數目105 143 136個，Q30比為91.40%。為剔除Illumina平臺錯誤率對結果的影響，需對原始數據進行質量控制，包括去除低質量序列，去除接頭，質量過濾和長度過濾等獲得純化數據（Clean Data），以進行后續工作。質量控制后‘新露’與‘新萃豐’分別獲得121 633 052個和100 378 460個短序列，測序總有效短序列與初始下機測序的比例分別為95.34%，95.46%。

表1 測序數據匯總

2.2 短序列匹配統計

將過濾后獲得的高質量數據比對到核桃參考基因組，統計結果表明：‘新露’中總匹配的短序列數為121 332 322個，占所有短序列數的98.1%，重復序列占比15.15%，覆蓋全基因的深度為26.76 x，當覆蓋深度≥20時，覆蓋全基因組的百分比為75.32%；‘新萃豐’中總匹配的短序列數為99 895 457個，占所有短序列數的97.96%，重復序列占比14.02%，覆蓋全基因的深度為22.45 x，覆蓋深度≥20時，覆蓋全基因組的百分比為59.92%。

2.3 SNP統計分析

‘新露’與‘新萃豐’核桃樣品檢測后獲得SNP位點統計信息（見表2），得到變異位點所在區域及變異影響。‘新露’樣本的總SNP數量為4 700 220個，雜合SNPs數量為2 645 454個，純合SNPs數量為2 054 766個，編碼區非同義突變數量（nsSNP）為81 580個，密度占比1.74%，SNP位于內含子區域內（Intronic）總數為535 307個，在基因間區域內總數（Intergenic）為3 469 438個。

表2 ‘新露’與‘新萃豐’SNPs數目與類型檢測結果

‘新萃豐’樣本經質量過濾后獲得4 675 346個高質量的SNP位點，編碼區內非同義突變總數為83 811個，SNP位于內含子區域內總數為541 875個，位于基因間區域的SNP占比73.2%。

2.4 InDel統計分析

在基因組中檢測長度小于50 bp的小片段插入或缺失（InDel）并進行變異位置及類型的注釋，對‘新露’與‘新萃豐’的InDels變異位點信息進行統計（見表3）。‘新露’的總InDel數量為648 603個，編碼區內移碼插入總數（Frame shift insertion）和移碼缺失（Frameshift deletion）數量分別為3 560個、2090個，InDel在基因5’UTR內總數（5’UTR），3’UTR內總數（3’UTR）和InDel位于不同基因的5’UTR和3’UTR內總數（5’UTR/3’UTR）為3 095個、5 912個、2個，InDel位于內含子區域內總數為106 672個，在基因間區域內總數為420 118個。

表3 ‘新露’和‘新萃豐’InDels數目與類型檢測結果

‘新萃豐’鑒定出606 188個InDel位點，其中插入297 690個，缺失311 822個，編碼區內移碼插入總數和移碼缺失數量為3 449個、2 015個。InDels大部分位于基因間區域，有389 525個（64.2%），次之為內含子 202 381個（16.89%），外顯子 8 048個（1.33%）。

2.5 DNA水平變異的基因分析

發生在外顯子區的變異可能會引起基因功能的異常改變，通過與參考基因組比對，檢測到‘新露’和‘新萃豐’基因組間發生非同義突變、同義突變和移碼突變等信息，分別獲得16 686個、19 829個變異SNPs以及1 687個、2 143個InDels，將這些變異位點比對到Gene ontology數據庫中。

2.5.1 ‘新露’核桃變異基因的GO富集分析

‘新露’核桃的SNP和InDel變異位點共注釋到4 865個基因，富集分析注釋結果分為生物過程、細胞組分和分子功能三類（如圖1）。注釋結果中，草酸代謝過程、囊泡運輸與錳離子結合分別在生物過程、細胞組分和分子功能類別中被顯著富集。在這些被顯著富集的三類反應中草酸代謝過程在植物體內廣泛存在，作用于植物抗逆過程中；高爾基體的囊泡運輸參與到植物細胞壁的形成，說明露仁性狀可能和次生細胞壁形成有關；而錳離子則是植物體細胞中許多酶的催化劑。

圖1 ‘新露’變異基因的GO注釋分類圖

2.5.2 ‘新萃豐’核桃變異基因的GO富集分析

‘新萃豐’核桃共注釋到6 131個基因，GO terms富集分析如圖2所示，在生物過程和細胞組分中微管過程和囊泡運輸被顯著富集，用于維持細胞骨架和信號傳導。在分子功能中富集到了天冬氨酸和苯丙氨酸的2種轉移酶，天冬氨酸轉氨酶，苯丙氨酸轉氨酶在植物體內分別作用于苯丙烷途徑前后，最終目的為合成木質素的單體。為進一步對變異基因的生物學功能進行了解，對富集到的天冬氨酸轉氨酶和苯丙氨酸氨基轉移后的裂解反應的基因進行挖掘，篩選出7個相關基因，結果見表4。

表4 木質素合成相關變異基因

圖2 ‘新萃豐’變異基因的GO注釋分類圖

3 討論

隨著測序技術的快速發展，蘋果［11］、柑橘［12］、梨［13］等多種果樹的全基因組序列陸續公布，不同果樹物種的基因組信息逐漸被完整破譯。通過全基因組重測序可以找到大量的SNP、拷貝數變異（CNV）、InDel、結構變異（SV）等變異信息，基于檢測到的變異位點能進一步研究該物種的特性［14］、群體進化［15］和定位目標性狀基因位點［16］。

SNP在基因組中廣泛存在，并含有豐富的遺傳信息，隨著測序技術的發展以及生物信息學分析能力的提升，SNP作為第三代分析標記應用于相關性狀的基因定位研究中［17］。SNP的最早研究用于人類基因組，最終繪制成142萬個SNP的人類基因組圖譜，這種高密度的SNP圖譜為整個基因組單倍型變異提供了依據［18］。將2個橄欖品種雜交鑒定獲得10 941個SNP，并使用其中多態性的SNP位點構建覆蓋基因組3 049 cM，25個連鎖群體的高密度遺傳連鎖圖譜，遺傳連鎖圖譜的構建可更有效地定位數量性狀和其他重要經濟性狀［19］。使用基因測序方法進行基因分型，對180種棗進行全基因組關聯分析（GWAS），鑒定出4 651個高質量SNP和45個與SNP性狀關聯基因，為標記輔助育種提供了數據基礎［20］。將SNP作為新型分子標記也應用于核桃分子育種中，將14 761個SNP進行全基因組關聯，確定了葉芽萌發和花芽分化等性狀的標記位點，為后續QTL鑒定奠定了基礎［21］。在自然適應和選擇性育種過程中，最常見且研究最多的SNP外還有InDel，InDel變異是不同個體間基因組同一位點的序列發生不等大小片段的插入與缺失，從而使同源序列比對產生空位（gap）的現象［22］。有研究者將SNP、InDel、SV以及簡單序列重復(SSR)組裝成核桃基因組變異圖譜，利用遺傳信息豐富的標記位點對中國五個核桃品種進行系統發育研究，將5個核桃種以100%的bootstrap(BS)值劃分為兩個已知區段核桃/胡桃組(Juglans/Dioscaryon 和核桃楸組 Cardiocaryon)［23］。在核桃上使用454焦磷酸測序技術開發鑒定獲得48 165個SNPs和1 037個InDels，轉換/顛換（ts/tv）比例為2.79:1，略高于本次新疆核桃測序的比例［24］。本研究選用具有代表性的新疆露仁品種與薄皮中的厚殼品種的‘新露’與‘新萃豐’，利用基因組重測序技術檢測核桃內果皮果殼發育不全，果仁裸露這一性狀的各種遺傳變異。‘新露’與‘新萃豐’的總SNP數目為4 700 220個，4 675 346個，總InDel數目為648 603個，606 188個，變異位點在基因編碼區，基因間區域和內含子區域都有分布，在全基因組水平上有著顯著的遺傳變異，2者的SNP數量顯著多于InDel的變異數量，說明外界環境對基因組的變化主要以單核苷酸多態性變異為主，這與梨［25］、馬鈴薯［26］等測序結果一致。本研究將檢測到的‘新露’與‘新萃豐’核桃基因組間所發生的非同義突變、同義突變等基因比對到GO數據庫中，發現在薄皮不露仁系列中的厚殼性狀品種‘新萃豐’富集到與天冬氨酸轉氨酶和苯丙氨酸轉氨酶的相關過程，通過與基因組序列比對分析，篩選到7個與木質素生物合成有關的基因。本試驗通過全基因組重測序技術對新疆2個核桃栽培種進行變異檢測分析，通過比對變異位點，挖掘變異的關鍵基因，為基因定位和克隆功能驗證以及為核桃分子育種的利用奠定了基礎。