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

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

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

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

1 材料與方法

1.1 試驗材料

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

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

1.2 試驗方法

1.2.1 DNA提取及重測序

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

1.2.2 基因組變異檢測和注釋

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

1.2.3 mSNP位點篩選和引物設(shè)計

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

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

續(xù)表2

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

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

1.3 數(shù)據(jù)處理及分析

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

2 結(jié)果與分析

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

對18份白化茶樹資源進行全基因組重測序，測序數(shù)據(jù)質(zhì)控分析結(jié)果顯示，除了 reads的前幾個堿基存在正常的不平衡現(xiàn)象之外，堿基質(zhì)量分布基本無AT、GC分離現(xiàn)象，表明測序質(zhì)量良好，可用于后續(xù)的分析。在測序深度為5×情況下，在各樣本中得到的原始測序數(shù)據(jù)在 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）。結(jié)果表明，18份白化茶樹資源的測序數(shù)量充足，質(zhì)量合格，可用于后續(xù)分析。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討論

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

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

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