基于全基因組重測序的白化茶樹mSNP標記開發及驗證

劉浩然，張晨禹，龔洋，葉圓圓，陳杰丹，陳亮，劉丁丁*，馬春雷*

1.中國農業科學院茶葉研究所/農業農村部特種經濟動植物生物學與遺傳育種重點實驗室，浙江 杭州 310008；2.中國農業科學院研究生院，北京 100081

白化茶樹新梢相較于正常綠色茶樹品種具有更為獨特的外形及較高的氨基酸和較低的茶多酚含量，是用于制備名優綠茶的優質原料[1]。近年來，選育的白化茶樹品種日漸增多，如中黃1號、黃金芽等[2-4]。然而，受限于茶樹復雜的遺傳背景及較高的雜合度，有關茶樹白化產生的原因尚未解析清楚。同時，白化現象極易受環境影響，為白化茶樹品種的鑒定與選育工作增加了難度。因此，研究白化茶樹種質資源遺傳多樣性及親緣關系，并開發相應的分子標記，對白化茶樹種質資源鑒定與利用具有重要意義。

基于DNA水平的分子標記已在茶樹種質資源的遺傳多樣性分析以及親緣關系鑒定中得到廣泛應用[5-8]。其中單核苷酸多態性（Single nucleotide polymorphism，SNP）標記具有數量多、分布廣、易于分型等特點，目前已作為主流的分子標記用于茶樹種質鑒定、遺傳連鎖圖譜構建等研究中[9-11]。隨著高通量測序技術的快速發展，近年來一種包含多個SNP位點信息的多聚單核苷酸多態性（Multiple nucleotide polymorphism，mSNP）標記得到開發并應用，該標記可以在單個擴增子內檢測多個SNP，從而提高單對引物擴增所獲得的DNA片段信息[12]。此外，基于多重PCR擴增的Genoplexs與基于液相探針捕獲的GenoBaits的靶向測序基因型檢測（Genotyping by target sequencing，GBTS）技術也得到快速發展，該檢測技術通過對基因組DNA中選取的特定靶向位點進行測序和基因型檢測，能夠減少DNA測序量及數據分析量，從而提高檢測效率[13]。而mSNP標記可依托于任意一種GBTS技術實現對目標位點的靶向捕獲，有效地提高了位點的檢出效率及標記的利用率，相關研究已在玉米[14-15]、豬[16-17]等動植物種質資源鑒定、分子育種服務等工作中被廣泛報道。然而在茶學領域，目前尚未見有關mSNP標記與GBTS技術的應用。因此，本研究利用全基因組重測序對18份白化茶樹資源進行變異檢測，并基于變異檢測信息與GBTS技術開發了一套包含59個mSNP分子標記的液相芯片，旨在為GBTS技術與mSNP標記在茶樹資源鑒定、遺傳多樣性及親緣關系分析等工作中的應用提供支撐。

1 材料與方法

1.1 試驗材料

以保存于國家茶樹種質圃（杭州）內的18份白化茶樹資源作為重測序材料（表1），于2020年春季采摘初展的一芽二葉新梢。立即用液氮固定，之后保存于-80 ℃冰箱，用于后續的DNA提取。

表1 供試茶樹品種的名稱及產地Table 1 The name and origin of tested tea cultivars

1.2 試驗方法

1.2.1 DNA提取及重測序

采用天根植物基因組提取試劑盒提取18份白化茶樹資源的基因組DNA，DNA產物送至石家莊博瑞迪生物技術公司，用于全基因組重測序。利用Qubit儀器（Thermo Fisher Scientific，UA）檢測DNA的質量和濃度，樣品基因組DNA檢測合格后，利用Covaris System（Covaris M220）超聲波打斷儀將基因組DNA片段化，長度在250 bp左右，然后利用全基因組文庫構建試劑盒構建高通量測序文庫，并利用Ion Torrent S5高通量測序儀（Thermo Fisher Scientific，UA）進行測序。

1.2.2 基因組變異檢測和注釋

測序完成后對原始數據進行質控，去除接頭和低質量的reads，并利用BWA軟件將clean reads與舒茶早參考基因組序列（http://tpia.teaplant.org/download.html）進行比對?；赾lean reads在參考基因組上的定位結果，使用Picard軟件過濾冗余的 reads，保證檢測結果的準確性。最后利用GATK軟件進行SNP和Indel的檢測，并利用 SnpEff軟件對 SNP與Indel進行注釋。

1.2.3 mSNP位點篩選和引物設計

基于重測序變異數據文件，按照 SNP測序深度不低于5×，SNP雜合率＜10%，缺失比例＜20%，最小等位基因突變頻率MAF（Minor allele frequency）＞0.01的標準進行位點過濾，獲得候選SNP位點集；隨后以250 bp為1個滑動窗口，統計各個窗口內的SNP位點數目，獲得2＜SNP數＜10的目標候選區段集；進一步計算每個目標候選區段集的多態信息含量（Polymorphism information content，PIC）值，選擇PIC＞0.5，每個目標候選區段集中至少有1個SNP位點且MAF＞0.35的區段作為候選目標區段集。對開發的 mSNP位點利用 Life Technology公司多重引物在線設計頁面（https://ampliseq.com）進行引物設計，最終篩選出59對mSNP引物，用于液相芯片的開發（表2）。

表2 mSNP引物序列Table 2 mSNP primer sequences

續表2

1.2.4 液相芯片對測試樣品的檢測

向定量質檢合格的茶樹資源DNA中加入多重PCR Panel mix和多重PCR擴增酶體系進行PCR擴增。擴增產物利用羧基磁珠進行純化后，加入帶有Barcode的測序引物和高保真PCR反應體系再次進行 PCR擴增，不同的Barcode用于區分不同的樣品。經過羧基磁珠純化后擴增產物，完成多重PCR捕獲及文庫建庫。文庫構建完成后，使用Qubit 2.0熒光計（Thermo Fisher Scientific，UA）進行初步定量，使用qPCR的方法對文庫的有效濃度進行準確定量以保證文庫質量。文庫質檢合格后，進行上機測序。

1.3 數據處理及分析

利用MEGA軟件計算18個白化茶樹資源基于重測序數據和mSNP標記獲得的Nei’s遺傳距離，并構建Neighbor-Joining（NJ）系統進化樹，進行聚類分析。利用Excel統計重測序數據，并根據液相芯片檢測的基因分型數據，計算測試材料與對照材料的遺傳相似度。遺傳相似度按照測試材料與對照材料基因型相同的標記位點數目除以測試材料與對照材料均檢出的標記位點數目進行計算。利用TBTOOLS繪制18個白化茶樹資源的共有非同義突變基因的upset圖。利用GraphPad Prism 8.0.1軟件完成柱形圖繪制。

2 結果與分析

2.1 18份白化茶樹資源的基因組重測序分析

對18份白化茶樹資源進行全基因組重測序，測序數據質控分析結果顯示，除了 reads的前幾個堿基存在正常的不平衡現象之外，堿基質量分布基本無AT、GC分離現象，表明測序質量良好，可用于后續的分析。在測序深度為5×情況下，在各樣本中得到的原始測序數據在 165 975 048～249 259 370 bp，clean data在 163 009 946～245 996 848 bp，其中 BJG-F1-1樣本最少，BY-1樣本最多；Q30在76.01%～87.73%。18份白化茶樹資源與參考基因組的比對率為63.95%～98.64%，平均覆蓋深度為5.99×（表3）。結果表明，18份白化茶樹資源的測序數量充足，質量合格，可用于后續分析。

表3 18份白化茶樹品種重測序數據及比對參考基因組數據結果Table 3 Resequencing data of 18 albino tea cultivars and comparison results of reference genome

2.2 不同白化茶樹資源間的SNP變異檢測及注釋

通過對 18份白化茶樹資源進行 SNP檢測，SNP數量分布在 42 485 966～49 127 511個，BY-1中的SNP數目最少，BY-2中的SNP數目最多。通過對 SNP位點的突變類型進行分析，結果顯示，相較于舒茶早參考基因組，18份白化茶樹資源中發生轉換（Transtition，Ti）類型的SNP數量為113 490 437個，其中G突變為A發生最多；發生顛換（Transversion，Tv）類型的SNP數量有42 095 202個，其中T突變為A發生最多。整體 Ti/Tv的比率為3.05～3.10（表4）。進一步對SNP位點在基因組上的位置進行分析，結果表明，分布于基因間區的 SNP數量最多（圖 1A），在各樣本間的比例在81.20%～81.60%；其次是位于內含子區，比例分布在3.94%～4.05%，而位于CDS區的突變比例僅占總數的0.73%，其中非同義突變在CDS區域的比例較高，在各樣本中分布在60.77%～62.94%（圖1B）。

圖1 18份白化茶樹資源的SNP注釋結果Fig.1 SNP annotation results of 18 tea resources

表4 不同白化茶樹資源間的SNP變異類型統計Table 4 Statistic of SNP variation types in different albino tea resources

2.3 不同白化茶樹資源的聚類分析

基于18份白化茶樹資源的全基因組重測序數據，對18份白化茶樹資源進行聚類分析。如圖 2所示，18份白化茶樹資源可劃分為3類。類群Ⅰ包含建德白茶、黃金葉等6份種質資源，其中白葉1號與黃金葉都來自浙江安吉，而建德白茶、景白1號、中黃1號、中黃2號分別由浙江建德、景寧、天臺、縉云的群體種中選育而來，聚類結果說明它們遺傳基礎較為接近。值得注意的是，在Ⅰ類聚類結果中，來自不同產地的建德白茶與黃金葉，中黃1號與中黃2號具有最為接近的遺傳關系，說明他們可能有著極為相似的遺傳背景。類群Ⅱ中包含BY-3、黃金菊等4份種質資源，其中江西的黃金菊與本課題組選育的白化株系花月遺傳距離最為接近，根據王松琳等[18]的觀察，兩者具有相似的表型性狀。類群Ⅲ包含天臺白茶、BY-1等8份種質資源，其中BJG-F1-1、BJG-F1-2、中茶129、白雞冠有著明確的親緣關系，聚為一類，中黃3號與黃金芽聚為一類。而來自于浙江天臺的天臺白茶及本課題組發現的芽變單株BY-1雖同聚在類群Ⅲ中，但與白雞冠、黃金芽親緣關系較遠。

圖2 基于重測序的18份白化茶樹資源的Neighbor-joining聚類Fig.2 Neighbor-joining clustering of 18 albino tea resources based on different markers

2.4 不同白化茶樹資源間DNA水平變異基因分析

編碼區的非同義突變可能會導致基因功能的改變。本研究篩選了不同白化茶樹資源間可能存在功能差異的基因。在18份白化茶樹資源中共篩選到17 056個共有非同義突變基因（圖3）。

圖3 18份白化茶樹資源共有非同義突變基因統計Fig.3 Statistics of non-synonymous mutant genes in 18 albino tea resources

進一步的GO（Gene ontology）富集分析結果顯示，在生物學過程、細胞組分、分子功能中共注釋到12 800個共有基因，其中在生物學過程中，非同義突變基因主要富集在代謝過程和細胞過程中；在細胞組分中，主要富集在膜和膜成分中；在分子功能中，主要富集在結合、催化活性部分（圖 4）。KEGG（Kyoto encyclopedia of genes and genomes）分析結果顯示，17 056個基因共富集到 135個代謝途徑，其中在代謝途徑和次級代謝產物的生物合成途徑中富集到的基因數量較多，分別為1 178個和654個。部分非同義突變基因顯著富集到植物病原相互作用、ABC轉運蛋白、甘油酯代謝等途徑中（圖4）。

圖4 18份白化茶樹資源中共有非同義突變基因的GO和KEGG富集分析Fig.4 GO and KEGG enrichment analysis of non-synonymous mutant common genes in different albino tea resources

研究表明，植物的白化現象與葉綠素合成代謝途徑密切相關[19-21]。為探究葉綠素代謝相關基因在不同白化茶樹資源間的序列差異，本研究基于18份白化茶樹資源的重測序數據提取了73個葉綠素合成相關基因的序列信息，突變位點比對分析顯示，在L-谷氨酸-tRNA、糞卟啉原-Ⅲ氧化酶、原卟啉還原酶等 14個與葉綠素合成相關的基因中共有98個錯義突變位點，這些基因中錯義位點的改變可能是導致部分白化茶樹資源葉色變異的重要原因（圖5）。

圖5 14個葉綠素合成相關基因錯義突變位點匯總Fig.5 Summary of missense mutation sites of 14 chlorophyll synthesis related genes

2.5 基于重測序的茶樹mSNP標記位點的篩選

對前期鑒定到的mSNP區段，按照在染色體上均勻分布的原則，優先選擇CDS區、PIC值高的區段作為目標mSNP分子標記位點，初步篩選到64個mSNP目標區段，區段內共包括259個SNP位點；然后提取這些目標區段，設計多重靶向PCR引物對18份白化茶樹樣本進行位點測試，過濾掉穩定性與多態性較差的SNP位點，最終獲得了 1個包含 222個 SNP位點、59個mSNP標記的液相芯片。進一步基于獲得的mSNP位點變異信息對18份白化茶樹資源進行遺傳多樣性分析。18份材料檢出的mSNP標記的PIC值主要分布在0.4～0.8區間（圖6A），表明這些標記具有較高的多態性。將18份資源基于59個mSNP位點的進化樹與基于全基因重測序數據的進化樹進行比對，發現二者的聚類結果相似，都可將18份種質資源分為3類，其中Ⅱ類的聚類結果完全一致，將 BY-2、BY-3、花月、黃金菊聚為一類，表明本研究開發的包含59個mSNP位點的液相芯片可高效檢測這18份白化茶樹資源的遺傳距離。

2.6 59個mSNP標記液相芯片的品種鑒別能力驗證

為檢驗該液相芯片的穩定性及品種鑒別能力，以開發液相芯片時所用的18份白化茶樹資源為對照，重新采集了9份開發芯片時所用資源（白葉1號、黃金菊、黃金葉、花月、景白1號、建德白茶、中黃1號、中黃2號、中黃3號）及4份新資源（翠峰、嘉茗1號、茂綠、福鼎大白茶）作為測試品種。利用液相芯片對其進行基因型檢測，將測試品種的基因型與對照品種的基因型進行比較，判斷該芯片對品種的鑒別能力。檢測結果表明，相同品種間一致性位點數在206～217個（圖6B）。遺傳分析結果表明，13份測試材料與18份對照材料兩兩比較結果差異較為明顯，其中當測試品種與對照品種為同一材料時，兩者間的遺傳相似度均高于90%，且主要位于96%～98%（圖6C）；當檢測品種與對照品種為不同品種時其遺傳相似度大多低于 84%，且主要分布在70%～80%（圖6C）。以上結果表明：（1）本研究開發的液相芯片能夠準確區分測試品種是否為對照品種，可實際應用于相關品種的資源鑒定；（2）利用該液相芯片進行品種鑒別時，可以以測試品種與對照品種之間的遺傳相似度是否高于90%為標準，用于判斷測試品種與對照品種是否為同一品種。

圖6 mSNP標記開發與檢驗分析Fig.6 mSNP marker development and testing analysis

3 討論

應用全基因組重測序技術可獲取全基因組范圍內的位點變異信息，目前已廣泛應用于基因定位、遺傳圖譜構建、突變位點鑒定、遺傳進化研究[22-23]。本研究通過對18份白化茶樹資源進行重測序，在測序深度 5×情況下，共獲得155 585 639個SNP位點。利用這些覆蓋全基因組的SNP位點對18份白化材料進行遺傳多樣性分析，結果表明，18份白化茶樹資源可聚為3類，其中具有親緣關系的白雞冠及其子代聚為一簇，建德白茶、黃金葉、白葉1號等表現出較為接近的進化關系，聚類結果與本課題組前期采用SSR標記的研究結果相似，總體呈現出有明確親緣關系或原產地地理位置較近的品種資源最先聚類的趨勢[18]。進一步對檢測到的所有SNP位點進行功能注釋及非同義突變分析，最終在18份白化資源中共獲得17 056個含有非同義突變位點基因。對含有非同義突變位點基因的功能富集進行分析，結果表明，這些基因主要富集在茶樹的代謝途徑和次級代謝產物的生物合成途徑中，特別是葉綠素合成途徑中存在多個非同義突變基因。葉綠素合成途徑中包含27種基因編碼的15種酶[24-25]，已研究表明原葉綠素酸脂氧化還原酶、葉綠素酶、谷氨酸-tRNA還原酶等基因突變可能是導致部分白化茶樹品種產生葉色變異的重要因素[24,26-27]。因此，為探究這些突變位點是否與18份茶樹資源的白化現象有關，本研究初步分析了14個葉綠素合成關鍵基因在18份白化茶樹CDS區上存在的98個錯義突變位點，如L-谷氨酸-tRNA中存在7個錯義突變位點、糞卟啉原-Ⅲ氧化酶存在3個錯義突變位點等（圖5），這些基因上的突變位點信息可為相關基因的功能研究提供參考。

目前，已有研究利用GBTS技術開發過多個高密度液相芯片，且已實際應用于作物品種選育及分子育種設計等工作中，如大豆40K液相芯片涵蓋了49個已知功能基因的85個探針位點，能夠有效地檢測功能基因的關鍵變異類型[28]；小麥中分別開發了16K、40K、60K的液相芯片，并挖掘出多個具有與小麥性狀相關的QTL位點信息[29-31]。除此之外，玉米、豬等動植物中也均開發了不同密度的液相芯片用于基因型的檢測[14,32-33]。這些芯片的開發極大地幫助了相關動植物在分子標記保護、品種判別等方面的工作[34]。適時將相關研究技術引入茶學領域，將有助于茶樹種質資源鑒別及品種保護工作?；诖?，本研究首次在茶樹中開發mSNP標記，通過對18份白化茶樹資源重測序檢測到的近一億六千個SNP標記進行篩選，前期先后篩選3套含260、390、650個mSNP標記的組合以供后續試驗選擇，經檢驗18份白化茶樹資源利用3套標記反映出的親緣關系與重測序得到的親緣關系基本一致，驗證了標記開發的可靠性。為進一步篩選變異位點實現對18份白化茶樹資源的鑒別，以260個組合的mSNP標記為基礎，通過進一步的區段篩選和mSNP檢測，最終獲得1個含59個mSNP標記的液相芯片。利用該液相芯片檢測13份茶樹資源的基因型信息，結果表明同一品種之間的遺傳相似度在92%～98%，不同品種間的遺傳相似度則在84%以下。由于檢測樣本偏少，因此初步將該芯片判斷樣品是否為同一品種的閾值線設為90%。

然而，由于茶樹基因組內存在大量的DNA重復序列[35-36]，以及測序技術自身存在的檢測誤差導致目前相同品種間仍然存在部分差異無法消除的情況[37]，后續有待通過測序技術的提升以及芯片密度的增加進一步降低檢測誤差?；诖?，本研究認為此次開發的含 59個 mSNP液相芯片能夠準確完成對 18份白化茶樹資源及其相關資源的品種鑒別工作，證實了mSNP的芯片可實際應用于茶學領域，后續將通過收集更多茶樹品種資源開發更高密度的液相芯片，用于完成對不同茶樹品種的鑒別與分析工作。