基于GBS測序的古藺野生茶樹遺傳多樣性與群體遺傳結構分析

王小萍,劉 飛,李明紅,張 廳,賴 謙,劉 曉,熊元元,唐曉波,李春華,王 云

(1.四川省農業(yè)科學院茶葉研究所,成都 610066;2.瀘州市經(jīng)濟作物站,四川 瀘州 646000;3.古藺縣農業(yè)農村局,四川 古藺 646500)

【研究意義】茶樹屬多年生異花授粉經(jīng)濟作物,起源于我國云貴高原,栽培歷史悠久,主要種植于我國華南、西南、江南以及江北茶區(qū)。四川盆地位于我國西南腹地,位于26°03′～34°19′ N,97°21′～108°31′ E,屬西南茶區(qū),地勢西高東低,氣候、生態(tài)類型復雜多樣,造就了區(qū)域內豐富的茶樹種質資源,包括地方品種、有性群體種、野生茶樹等。四川省野生茶樹主要分布在金沙江沿岸以及龍門、邛崍山脈一帶[1]。野生茶樹具有品質特異、抗性強的特點,是茶樹育種重要的物質基礎,因此研究野生茶樹的演化進程、群體遺傳結構,明確野生茶樹的遺傳多樣性和親緣關系,對優(yōu)良茶樹品種的選育具有重要意義。【前人研究進展】單核苷酸多態(tài)性(SNP)是基于簡化基因組測序(GBS)技術發(fā)展起來的第3代最具優(yōu)勢的分子標記技術,具有簡單、快速、特異性強、穩(wěn)定遺傳、便于檢測等優(yōu)點[2]。張鵬等[3]利用SNP芯片技術研究107份云南玉米自交系,發(fā)現(xiàn)云南地區(qū)玉米種質資源遺傳基礎豐富,可創(chuàng)建新的具有良好應用潛力的雜交優(yōu)勢群。孫倩等[4]采用SNP和Indel標記揭示了巴西木薯的遺傳多樣性。周新桐等[5]利用GBS測序技術構建了甘藍型油菜連鎖圖譜,并對油菜粉色花性狀進行了QTL定位。Hou等[6]基于GBS測序技術對信陽藍殼雞面部色素性狀進行了全基因組關聯(lián)分析,定位到10個與面部色素性狀緊密關聯(lián)的SNP標記。周艷春等[7]利用SNP標記對93份無芒雀麥種質資源產量性狀進行全基因組關聯(lián)分析,篩選出33個核心SNP(P<0.001),有利于縮短無芒雀麥的育種進程。由此可見,GBS測序技術研究結果可靠,已廣泛應用于動植物遺傳多樣性分析、群體遺傳結構分析、QTL定位、遺傳圖譜構建以及全基因組關聯(lián)分析等研究領域,且研究結果可靠?！颈狙芯壳腥朦c】古藺縣位于四川省南部邊緣,位于27°41′～28°20′ N,105′34′～106°20′ E,與貴州畢節(jié)縣、赤水縣相鄰。前期調查發(fā)現(xiàn),古藺縣蘊含豐富的野生茶樹,多以喬木、小喬木為主,樹徑多大于0.3 m,最大可達2.0 m,且這些野生茶樹具有較強的抗性或優(yōu)異的品質性狀,是茶樹新品種選育的寶貴資源?！緮M解決的關鍵問題】為分析古藺當?shù)匾吧铇涞倪z傳多樣性與群體遺傳結構,揭示其遺傳分化特征,本研究以65份野生茶樹為研究對象,利用GBS測序技術,基于獲得的SNP分子標記研究參試茶樹的遺傳多樣性和群體遺傳結構,以期為四川省茶樹種質創(chuàng)新、優(yōu)良品種選育等提供理論依據(jù),同時也為今后優(yōu)良性狀基因深度挖掘提供參考。

1 材料與方法

1.1 試驗材料

供試野生茶樹共65份(表1),其中古藺野生茶樹49份(以G1、G2、G3......表示)、敘永野生茶樹5份(以X1、X2、X3......表示)、合江野生茶樹1份(以H1表示)、崇州枇杷茶4份(以CP1、CP2、CP3......表示)、雷波野生茶4份(以LB1、LB2、LB3......表示)、云南野生茶樹2份(以YD1、YD2表示)。于2019年9月取參試材料新鮮成熟葉片,-80 ℃液氮冷凍保存用于GBS測序。

表1 65份野生茶樹資源

1.2 茶樹基因組DNA提取及建庫測序

采用改良CTAB法提取茶樹基因組DNA,使用Qubit和Nanodrop檢測DNA質量。質檢合格的基因組DNA由基迪奧生物科技公司進行30×的GBS測序,利用Novaseq6000測序儀平臺,按照GBS實驗流程,對基因組DNA進行限制性內切酶酶切,加上帶barcode的測序接頭,對樣本進行混合,構建小片段文庫(250～550 bp),采用PE125雙末端測序。

1.3 SNP與InDel標記檢測

用 FASTP (版本 0.18.0)[8]對 Illumina 平臺的原始數(shù)據(jù)進行過濾,過濾標準如下:①去除含有未知核苷酸(N)≥10%的 reads;②去除低質量reads(質量值Q≤10 的堿基數(shù)占整條read的50%以上);③去除含接頭的reads。過濾后的clean reads 用于組裝分析。使用比對軟件 BWA(v0.7.12)[9]采用mem算法將過濾后的reads 比對到參考基因組上,比對參數(shù)為-k32-M;比對結束后使用軟件 picard(版本 1.129)對結果進行標記,并使用變異檢測軟件 GATK(版本 3.4-46)[101]進行群體變異檢測,利用ANNOVAR(版本2)[11]對檢測出的變異進行功能注釋。

1.4 野生大樹茶遺傳多樣性與遺傳結構分析

采用POPGENE軟件計算各地區(qū)茶樹遺傳多樣性參數(shù):多態(tài)性位點數(shù)(PLN)、多態(tài)位點百分比(PPL)、有效等位基因數(shù)(Ae)、觀測等位基因數(shù)(Na)、最小等位基因頻率(Maf)、觀測雜合度(Ho)、期望雜合度(He)以及多態(tài)信息含量(PIC)和Shannon’s指數(shù);利用VCFtools軟件[12]計算各居群的群體分化指數(shù)(Fst)。采用Plink和Admixture軟件[13]進行群體結構分析,估算出最佳的群體亞群數(shù)。

根據(jù)所得到的SNP信息,計算樣品間距離,對樣品進行聚類分析,利用treebest軟件[14]通過鄰接法(Neighbor-joining methods)構建SNP進化樹。并基于個體間的SNP差異程度,利用R語言包(http://www.r-project.org/)進行主成分分析,根據(jù)各個樣本在第一主成分(PC1)和第二主成分(PC2)兩個綜合指標中的數(shù)值大小,繪制二維坐標圖。

2 結果與分析

2.1 野生茶樹基因組測序結果與SNP注釋

對65份野生茶樹基因組DNA進行分析,構建GBS文庫及酶切測序,共獲得77.30 G數(shù)據(jù),過濾后獲得74.48 G數(shù)據(jù),540 633 846條高質量reads,平均每份茶樹資源獲得1.15 G數(shù)據(jù),8317 444條高質量reads。65份野生茶樹Q20≥93.30%,Q30≥91.84%,表明測序質量較好,GC含量在40.02%～41.71%,分布正常。將高質量reads比對到茶樹參考基因組上,65份野生茶樹共獲得799 217個變異,其中769 893個SNPs,29 324個Indels。

利用 ANNOVAR軟件對769 893個高質量的SNPs進行檢測(圖1-a),692 570個SNPs位于基因間區(qū)(Intergenic region),占比89.96%;31 173個SNPs位于內含子區(qū)(Intronic region),占比4.05%;19 754個SNPs位于外顯子區(qū)(Exonic region),占比2.57%;8994個SNPs位于轉錄終止位點下游區(qū)域(Downstream region),占比1.17%;7703個SNPs位于轉錄終止位點上游區(qū)域(Upstream region),占比1.0%;1801個SNPs位于3′端UTR變異(UTR3),占比0.23%;949個SNPs位于5′端UTR變異(UTR5),占比0.12%;812個SNPs位于轉錄終止位點(Downstream; Upstream),占比0.11%;137個SNPs位于可變剪切位點(Splicing region),占比0.02%。

a:SNPs位置;b:SNPs功能;c:SNPs變異類型。a: Location of SNPs;b:Function of SNPs;c:Types of SNPs.圖1 SNPS注釋結果Fig.1 The results of SNP annotation

對編碼區(qū)SNPs進行功能注釋(圖1-b),其中11 791個SNPs屬非同義突變(Nonsynonymous),7640個SNPs屬同義突變(Synonymous),311個SNPs屬無義突變(Stop gain),12個SNPs屬起始缺失(Stop loss)。

對SNP變異類型進行注釋(圖1-c),76.08%的SNPs屬于堿基轉換(Transition,ts),其中G->A占比26.47%,C->T占比26.53%,T->C占比11.65%,A->G占比11.61%;23.92%的SNPs屬于堿基顛換(Transversion,tv),其中A->C與A->T分別占比2.57%與3.48%,G->C與G->T分別占比2.11%與3.76%,C->A與C->G分別占比3.81%與2.11%,T->A與T->G則分別占比3.52%與2.57%。表明SNP變異類型以堿基轉換為主,且ts/tv值為3.18,表明SNP變異類型多以同類型核苷酸之間的變異為主。

2.2 野生茶樹群體遺傳多樣性與遺傳結構分析

從表2可知,6個群體平均有效等位基因數(shù)(Ne/Ae)為1.1284個,平均多態(tài)性位點數(shù)(PLN)為9475 456,最高為Gulin群體(2828 344),多態(tài)性位點比例(PPL)為57.33%,最低為合江野生茶樹(172 911),多態(tài)性位點比例為8.53%;各群體觀測等位基因數(shù)(Na)分別在0.4267～0.9147,平均為0.7677;最小等位基因頻率(Maf)介于0.1514～0.5000,均值為0.3213;基因多樣性指數(shù)(Nei)最高是古藺野生茶樹(0.1328),最低是云南野生茶樹(0.0632),均值為0.0979;各群體觀測雜合度(Ho)與期望雜合度(He)僅古藺地區(qū)變化范圍稍大,為0.0636～0.1233,其余變化范圍較小;多態(tài)性信息含量(PIC)則介于0.0312～0.1025,均值為0.0596;香農指數(shù)(Shi)最高的為古藺地區(qū)(0.199),最低的為云南地區(qū)(0.0576),均值為0.1122。

表2 參試野生茶樹群體遺傳多樣性

從表3可知,不同地區(qū)野生茶樹群體間遺傳分化指數(shù)Fst變化范圍為-0.7042～0.4372,云南野生茶樹與崇州枇杷茶、古藺野生茶樹、雷波野生茶樹、敘永野生茶樹和合江野生茶樹群體間遺傳分化指數(shù)分別是0.3898、0.3375、0.3673、0.4372和0.3094,表明群體間存在非常強的遺傳分化,其中與敘永群體間遺傳分化最高,為0.4372;合江野生茶樹與崇州枇杷茶、古藺野生茶樹、雷波野生茶樹、敘永野生茶樹種群間遺傳分化均為負數(shù),表明這些群體間存在普遍的基因交流,遺傳親緣關系較近,其中與古藺群體遺傳分化值最小,為-0.7042;古藺野生茶樹與崇州枇杷茶、雷波野生茶樹、敘永野生茶樹種群間遺傳分化指數(shù)分別為0.04373、0.11740和-0.03481,表明古藺群體與崇州枇杷茶群體間存在遺傳分化,但分化較弱,與雷波群體間存在中等程度的遺傳分化,而與敘永群體間則存在基因交流。

表3 野生茶樹群體遺傳分化指數(shù)(Fst)

利用769 893個高質量的SNPs變異,采用Plink和Admixture軟件對65份野生茶樹進行群體結構分析。將亞群數(shù)K值范圍設定為2～9,計算不同K值下的交叉驗證錯誤率(Cross-validation error,CV error),由圖2-a可知,當K值為2時,CV值最小,因此可將參試的65份茶樹種質分成2個亞群,再根據(jù)Q值的大小,將野生茶樹劃分到最大Q值所在的亞群(圖2-b)。2個亞群中CP4、G19、LB5、YD1、YD2共5份野生茶樹歸為亞群Ⅰ(紅色),其余60份野生茶樹歸為亞群Ⅱ(藍色)。

a:計算K為2～9時的CV error。b:K=2時的群體結構,在該群體結構中,紅、藍分別代表亞群1～2,根據(jù)顏色占比推斷該野生茶樹屬于哪個亞群。a:The cross-validation error was calculated when K was 2 to 9.b:The Population structure at K=2,in this Population structure,red,blue represented subgroups 1 to 2,and the subgroup to which the individual belonged was inferred based on the proportion of the color.圖2 65份野生茶樹群體遺傳結構Fig.2 Population genetic structure of 65 wild tea resources

2.3 野生茶樹的SNP系統(tǒng)進化分析

利用769 893個高質量SNPs計算樣品間距離,并對樣品進行聚類分析,從而推斷出種群間的親緣關系遠近。利用treebest軟件鄰接法(Neighbor-joining methods)構建進化樹(圖3)。65份野生茶樹可以明顯分為兩大群,其中YD2、YD1、CP4、LB5、G19、G10為Ⅰ類群,其余59份為Ⅱ類群,這與遺傳結構分析結果存在差異,系統(tǒng)進化分析將G10聚在Ⅰ類群,而在遺傳結構分析中G10與大部分野生茶樹聚為一類。同時,Ⅱ類群可以分為3個亞群(Ⅱ-a、Ⅱ-b和Ⅱ-c),其中亞群Ⅱ-a僅包含材料G65,Ⅱ-b包含G67、G69、G70、G71和G725 5份野生茶樹,均源自相近地理位置,Ⅱ-c則包含其余53份野生茶樹。

紫色為古藺野生茶樹;紅色為敘永野生茶樹;暗紅色為雷波野生茶樹;淺藍色為合江野生茶樹;藍色為崇州枇杷茶;淺綠色為云南野生茶樹。Purple was wild tea resources in Gulin;Red was wild tea resources in Xuyong;Dull red was wild tea resources in Leibo; Light blue was wild tea resources in Hejiang;Blue was Chongzhou Pipa tea; Light green was wild tea resources in Yunnan.圖3 65份野生茶樹系統(tǒng)發(fā)育進化樹Fig.3 Phylogenetic tree of the 65 wild tea resources

2.4 PCA主成分分析

基于個體間SNPs差異程度,利用R軟件進行主成分分析,分別利用PCA1與PCA2、PCA1與PCA3以及PCA2與PCA3繪制散點圖(圖4),發(fā)現(xiàn)PCA1與PCA2最能解釋65份野生茶樹的遺傳變異情況,PCA1vPCA3和PC2vPC3次之。由圖4-a可知,65份野生茶樹能明顯區(qū)分為4個亞群,亞群1包含42份古藺野生茶樹、5份敘永野生茶樹、1份合江野生茶樹、3份雷波野生茶樹以及3份崇州枇杷茶,亞群2包含3份茶樹茶樹,分別來自古藺、雷波、崇州,間于云南與大部分古藺茶樹集中點間,亞群3包含6份古藺野生茶樹,分布在Y軸上方,與亞群1大部分古藺野生茶樹集中點距離較遠,亞群4則是2份云南野生茶樹,和其余茶樹區(qū)分明顯。

a:PC1vPC2; b:PC1vPC3;c:PC2vPC3.圖4 65份野生茶樹SNP基因型主成分分析Fig.4 Principal components analysis of SNP genotype of the 65 wild tea resources

3 討 論

SNP是最常見的一種遺傳變異,具有數(shù)量多、分布廣的特點,通過比較分析種群內個體的SNP位點,可研究種群的遺傳多樣性、親緣關系以及群體遺傳結構。王麗鴛等[15]研究發(fā)現(xiàn)茶樹SNP的分布規(guī)律,即茶樹SNP分布頻率高,變異類型主要是G/A與C/T的堿基轉換為主。郭燕等[16]利用GBS技術對貴州久安100份古茶樹資源開展全基因組簡化測序分析,共獲得548 597個高質量SNP,其中G/A和C/T的SNP變異類型最多,且ts/tv為2.95。本研究利用GBS技術對65份野生茶樹進行簡化基因組測序分析,獲得769 893個高質量SNPs,變異類型以G/A與C/T為主,這與上述結果一致,但ts/tv為3.18,較上述研究高,表明參試茶樹SNP變異類型較豐富。通過對外顯子區(qū)域(基因編碼區(qū))19 754個SNP進行功能注釋,SNP導致的非同義突變與同義突變的比例為1.54,表明參試茶樹在進化過程中受到自然環(huán)境的正選擇效應,環(huán)境適應能力強,同時這些變異為后續(xù)SNP標記的開發(fā)、功能基因的研究提供基礎。

遺傳多樣性是作物群體內個體基因水平的變異與環(huán)境等多因素共同作用的結果,是物種適應環(huán)境、保持遺傳潛力的基礎,對種群的進化發(fā)展具有重要意義[17]。研究表明,雷波野生茶樹[18]、云南野生茶樹[19-20]、江華苦茶群體種[21]、貴州野生茶樹[22]以及廣西野生茶樹[23]均具有較高的遺傳多樣性,而本研究中基于獲取的SNP對參試茶樹的遺傳多樣性進行分析,結果表明茶樹群體的遺傳多樣性參數(shù)與之相比處于較低水平(基因多樣性指數(shù)Nai=0.0979、期望雜合度 He=0.0745、多態(tài)信息含量PIC=0.0596、Shannon信息指數(shù)=0.1122)。推測原因可能與供試茶樹的數(shù)量和種群來源有關,本研究的65份野生茶樹源于四川省古藺縣(49份)、合江縣(1份)、敘永縣(5份)、雷波縣(4份)、崇州(枇杷茶,4份)和云南省德宏芒市(2份),其中古藺、合江和敘永地理位置相鄰,氣候環(huán)境較一致,在文中后續(xù)遺傳分化研究中也發(fā)現(xiàn)這3個地方野生茶樹在遺傳演化的過程中,存在頻繁的基因交流,遺傳親緣關系近,而源于這3個地區(qū)野生茶樹達55份,占比85%,這可能是導致本研究總體遺傳多樣性低的原因之一。另外前人研究多選用多態(tài)性豐富的SSR標記對茶樹遺傳多樣性進行研究,而多態(tài)性豐富的標記獲得的等位基因變異數(shù)較多,檢測到的有效等位基因、雜合度、香農指數(shù)等參數(shù)也較高[24-25],這可能也是研究獲得較高遺傳多樣性參數(shù)的原因。何旭東等[26]利用SLAF-seq測序技術研究舒瑪櫟遺傳多樣性時也獲得了較低的遺傳多樣性參數(shù),且遺傳多樣性指數(shù)遠低于利用SSR標記的研究結果。

群體結構分層、等位基因分布不均容易導致基因型與表型分析結果假陽性或假陰性,因此對參試群體進行群體結構研究是基因與表型關聯(lián)分析的基礎。本研究采用Plink和Admixture軟件對65份野生茶樹進行群體結構分析,當K=2時,CV值最小,由此將65份野生茶樹分成2個亞群,該結果與系統(tǒng)進化樹分析結果基本相符,而與主成分分析存在差異。主成分分析不僅將遺傳距離較遠的YD1、YD2劃分在一起,也將其余川內野生茶樹劃分為3個亞群,其中亞群2分布介于亞群1與亞群4之間,且劃分在亞群2的3份野生茶樹在群體結構和系統(tǒng)進化分析中則和YD1、YD2聚在一起,亞群3則是從古藺野生茶樹中劃出,因此主成分分析結果是對群體結構和系統(tǒng)進化分析的補充,更好地解釋了野生茶樹的遺傳結構。

我國西南云貴高原一直被認為是茶樹的原產地,四川省在地理位置上處于云貴高原的過渡地帶。本研究結果顯示65份野生茶樹資源來自2個血統(tǒng),四川崇慶、雷波、古藺3個地區(qū)各1份材料與云南材料血統(tǒng)一致,表明四川與云南茶樹資源存在演化進程;而17份古藺地區(qū)和1份雷波地區(qū)的茶樹組成另一血統(tǒng)且和前者區(qū)別明顯,無任何交集,表明四川省存在原產的野生大茶樹,與鐘渭基[1]的四川野生大茶樹系本地原產一致,也與周斌等[18]的“四川與云南的野生茶樹群落相互獨立”一致。本研究對四川古藺、雷波、敘永、合江、崇州5個茶樹群體結構研究中,四川純正血統(tǒng)的野生茶樹集中在古藺地區(qū),其余茶樹存在2個血統(tǒng)雜合現(xiàn)象,但主要以四川血統(tǒng)為主,表明古藺地區(qū)存在更原始的茶樹血統(tǒng),其余區(qū)域的野生茶樹均以此血統(tǒng)為主與另一血統(tǒng)雜交。筆者經(jīng)過近幾年實地考察發(fā)現(xiàn),古藺地區(qū)野生大茶樹數(shù)量較多,以原桂花鄉(xiāng)和黃荊鄉(xiāng)為中心廣泛分布,據(jù)此推測古藺地區(qū)及周邊可能是四川野生大茶樹的發(fā)源中心,這與周斌等[18]的雷波地區(qū)可能是川內野生茶樹的發(fā)源中心觀點不一致。因此,針對此觀點,需在今后研究中進一步擴大各地區(qū)試驗樣本數(shù)進行驗證。

4 結 論

65份古藺野生茶樹具有較低的遺傳多樣性,親緣關系較近,存在頻繁的基因交流;遺傳結構分析可將65份野生茶樹分為兩個亞群。