基于全基因組重測序的白化茶樹mSNP標(biāo)記開發(fā)及驗(yàn)證

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

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

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

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

1 材料與方法

1.1 試驗(yàn)材料

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

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

1.2 試驗(yàn)方法

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）進(jìn)行測序。

1.2.2 基因組變異檢測和注釋

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

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

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

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

續(xù)表2

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

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

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

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

2 結(jié)果與分析

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

對18份白化茶樹資源進(jìn)行全基因組重測序，測序數(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份白化茶樹資源進(jìn)行 SNP檢測，SNP數(shù)量分布在 42 485 966～49 127 511個，BY-1中的SNP數(shù)目最少，BY-2中的SNP數(shù)目最多。通過對 SNP位點(diǎn)的突變類型進(jìn)行分析，結(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）。進(jìn)一步對SNP位點(diǎn)在基因組上的位置進(jìn)行分析，結(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份白化茶樹資源進(jìn)行聚類分析。如圖 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)系較遠(yuǎn)。

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

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

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

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

進(jìn)一步的GO（Gene ontology）富集分析結(jié)果顯示，在生物學(xué)過程、細(xì)胞組分、分子功能中共注釋到12 800個共有基因，其中在生物學(xué)過程中，非同義突變基因主要富集在代謝過程和細(xì)胞過程中；在細(xì)胞組分中，主要富集在膜和膜成分中；在分子功能中，主要富集在結(jié)合、催化活性部分（圖 4）。KEGG（Kyoto encyclopedia of genes and genomes）分析結(jié)果顯示，17 056個基因共富集到 135個代謝途徑，其中在代謝途徑和次級代謝產(chǎn)物的生物合成途徑中富集到的基因數(shù)量較多，分別為1 178個和654個。部分非同義突變基因顯著富集到植物病原相互作用、ABC轉(zhuǎn)運(yù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)基因的序列信息，突變位點(diǎn)比對分析顯示，在L-谷氨酸-tRNA、糞卟啉原-Ⅲ氧化酶、原卟啉還原酶等 14個與葉綠素合成相關(guān)的基因中共有98個錯義突變位點(diǎn)，這些基因中錯義位點(diǎn)的改變可能是導(dǎo)致部分白化茶樹資源葉色變異的重要原因（圖5）。

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

2.5 基于重測序的茶樹mSNP標(biāo)記位點(diǎn)的篩選

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

2.6 59個mSNP標(biāo)記液相芯片的品種鑒別能力驗(yàn)證

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

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

3 討論

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

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

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