趙 陽，馮延芝，楊超偉，王保平，喬 杰，殷世雨，周海江，李芳東

（1.中國林業(yè)科學(xué)研究院 經(jīng)濟(jì)林研究所，河南 鄭州 450003；2.經(jīng)濟(jì)林種質(zhì)創(chuàng)新與利用國家林業(yè)和草原局重點實驗室，河南 鄭州 450003；3.信陽市南灣試驗林場，河南 信陽 464031）

泡桐Paulownia是玄參科Scrophulariaceae 泡桐屬Paulownia樹種的統(tǒng)稱，是原產(chǎn)我國的重要速生多用途樹種，栽培和利用歷史悠久[1]。泡桐具有材質(zhì)優(yōu)良、生長迅速和繁殖容易等優(yōu)良特性，在速生豐產(chǎn)林營造、城鄉(xiāng)居民區(qū)綠化和國民經(jīng)濟(jì)建設(shè)等領(lǐng)域都發(fā)揮了重要作用[2]。我國泡桐種質(zhì)資源豐富，自然分布范圍廣，為泡桐遺傳改良提供了有利條件[3]。為了高效開發(fā)和利用泡桐種質(zhì)資源，不斷提高新一輪次育種群體目標(biāo)性狀的相關(guān)基因頻率[4]，進(jìn)一步優(yōu)化交配設(shè)計[5]、實現(xiàn)增效基因的聚合[6]，從而促進(jìn)種質(zhì)創(chuàng)新與培育突破性良種，厘清泡桐屬植物的分類和系統(tǒng)發(fā)育關(guān)系是迫切需要解決的問題。

前期已有學(xué)者針對泡桐屬植物的分類問題開展了很多研究，除了基于形態(tài)學(xué)分類、數(shù)量分類和孢粉學(xué)分析等傳統(tǒng)方法外，還利用多種分子技術(shù)進(jìn)行分類[7-9]，試圖揭示泡桐屬植物的系統(tǒng)發(fā)育關(guān)系。然而，由于泡桐分布廣泛且多異花授粉，在漫長的自然遷移、頻繁的種間雜交和長期的人工栽培之后，屬內(nèi)種間親緣關(guān)系變得相當(dāng)復(fù)雜，不同研究方法取得的結(jié)果也不盡相同。比如，莫文娟等[3]利用ISSR 分子標(biāo)記分析了泡桐屬植物的親緣關(guān)系，并將其分為毛泡桐類（毛泡桐P.tomentosa、蘭考泡桐P.elongata、山明泡桐P.lampropylla、白花蘭考泡桐P.elongataf.alba、楸葉泡桐P.catalpifolia、圓冠泡桐Paulownia×henanensisC.Y.Zhang et Y.H.Zhao.Hybr.Nov.、宜昌泡桐P.ichangensis和亮葉毛泡桐P.tomentosavar.lucida）、白花泡桐類（白花泡桐P.fortunei、南方泡桐P.australis）和川泡桐類（臺灣泡桐P.taiwaniana、川泡桐P.fargesii）3 類；馬浩等[10]基于葉綠體DNA 的RFLP 分析，將泡桐屬15 個種分為南方泡桐組（南方泡桐和成都泡桐）、毛泡桐組（毛泡桐和蘭考泡桐）和白花泡桐組（白花泡桐、山明泡桐、海島泡桐、臺灣泡桐、建始泡桐、宜昌泡桐、興山泡桐、川泡桐、楸葉泡桐和白花蘭考泡桐）3 組；而盧龍斗等[11]基于RAPD 分析認(rèn)為白花泡桐和川泡桐的親緣關(guān)系在泡桐屬的7個桐種中最為接近。可以看出，前人的研究結(jié)果不盡相同，仍不能很好地解決泡桐屬內(nèi)部深層次的系統(tǒng)發(fā)育關(guān)系。

全基因組重測序（Whole genome re-sequencing）是基于DNA 分子水平對已知基因組序列的物種的不同個體進(jìn)行高通量測序，可以通過將測試品種的序列與已知參考序列進(jìn)行比對，確定測試品種遺傳構(gòu)成的差異[12]。單核苷酸多態(tài)性（Singlenucleotide polymorphisms，SNP）是指基因組中某個特定位置上發(fā)生單個核苷酸變異所引起的DNA 序列多態(tài)性，與ISSR（Inter-simple sequence repeat）、RFLP（Restriction fragment length polymorphism）、RAPD（Random amplified polymorphic DNA）和SCAR（Sequence characterized amplified region）等分子標(biāo)記相比，具有數(shù)量多、分布密度高、快速且高通量分型等優(yōu)點[13-14]。隨著測序技術(shù)的快速發(fā)展、擁有參考基因組物種數(shù)量的不斷增多以及基因分型成本的下降，SNP 標(biāo)記日漸廣泛地被應(yīng)用到植物群體遺傳學(xué)研究中[15-18]。白花泡桐全基因組測序的完成推動了SNP 標(biāo)記在泡桐研究中的應(yīng)用。李文楊等[19]基于篩選到的SNP 位點對由毛泡桐和白花泡桐為親本，連同正反交獲得的181 個F1代個體構(gòu)成的作圖群體進(jìn)行基因型分析，構(gòu)建了高密度的連鎖遺傳圖譜，為泡桐分子育種研究提供了數(shù)據(jù)參考。但是，尚未見到有關(guān)SNP 標(biāo)記在泡桐系統(tǒng)發(fā)育分析和分類地位研究的報道。因此，本研究在對泡桐屬11 個桐種的典型株進(jìn)行全基因組重測序的基礎(chǔ)上開發(fā)SNP 標(biāo)記，解析泡桐屬內(nèi)的親緣關(guān)系和遺傳結(jié)構(gòu)，以期為泡桐核心育種群體構(gòu)建、骨干親本篩選和種質(zhì)創(chuàng)新提供依據(jù)。

1 材料與方法

1.1 試驗材料

本試驗選用的11 份材料均采自中國林科院經(jīng)濟(jì)林研究所所屬、位于湖北省鐘祥市盤石嶺林場的泡桐種質(zhì)資源庫，涵蓋了泡桐屬所有種的典型株。采集葉片后用去離子水清洗，放入液氮中進(jìn)行速凍，置于冰箱中-80 ℃下保存?zhèn)溆?。試驗材料詳?xì)信息見表1。

表1 供試試驗材料Table 1 Materials used in the experiment

1.2 試驗方法

1.2.1 DNA 提取、建庫和測序

采用改良的CTAB 法提取葉片樣品的基因組DNA[20]，經(jīng)1%瓊脂糖凝膠電泳和紫外分光光度計檢測質(zhì)量和濃度后選取合格的DNA 樣品進(jìn)行GBS 文庫構(gòu)建。將DNA 用MseI 限制性內(nèi)切酶隨機(jī)打斷成短的DNA 片段后進(jìn)行末端修復(fù)，然后在DNA 片段兩端連接dA 尾并連接測序接頭。對加上接頭的DNA 片段經(jīng)過AMPure XP 磁珠純化后選擇300～400 bp 范圍的片段進(jìn)行PCR 擴(kuò)增和電泳回收純化，最后在Illumina HiSeq 4000 平臺上機(jī)進(jìn)行雙末端PE150 測序。

1.2.2 泡桐基因組SNPs 挖掘

測序獲得原始數(shù)據(jù)（Raw data）后，去除含接頭序列的reads、單端測序序列中未檢出的堿基超過該序列長度10%的序列和低質(zhì)量（≤5）堿基數(shù)超過該序列長度50%的序列后，得到高質(zhì)量的測序數(shù)據(jù)（Clean reads）。

以白花泡桐為參考基因組（Paulowniafortunei,https://www.ncbi.nlm.nih.gov/genome/?term=paulownia），采用BWA（Burrows-wheeler-alignment tool）軟件[21]將Clean reads 比對到參考基因組，然后使用picard 對比對結(jié)果標(biāo)記并去除重復(fù)序列。采用GATK 軟件進(jìn)行群體SNP 檢測，選擇貝葉斯模型檢測群體中的多態(tài)性位點，通過有效過濾得到高質(zhì)量SNPs 結(jié)果用于后續(xù)分析[22]。

1.2.3 數(shù)據(jù)分析

基于檢測獲得的SNP，采用VCF2Dis 軟件計算遺傳距離矩陣，采用MEGA11 軟件中鄰近法（Neighbor-joining method）構(gòu)建進(jìn)化樹，采用FigTree 軟件對進(jìn)化樹進(jìn)行可視化。采用plink 軟件計算特征向量及特征值[23]，然后使用Seaborn 可視化庫繪制PCA 分布圖，最后通過fastSTRUCTURE的默認(rèn)參數(shù)篩選K值并推斷種群結(jié)構(gòu)[24]。以每kb的變異數(shù)代表不同基因編碼序列（CDS）的突變率，采用四分位法將基因分組，采用clusterProfiler 包對各組基因進(jìn)行富集分析。

2 結(jié)果與分析

2.1 測序質(zhì)量

采用Illumina 測序平臺對泡桐屬11 個種的典型株進(jìn)行全基因組重測序，過濾掉不符合要求（含有接頭和低質(zhì)量的reads）的數(shù)據(jù)后，共獲得高質(zhì)量測序數(shù)據(jù)63.08 GB，平均每個樣本獲得5.87 GB 的數(shù)據(jù)。經(jīng)過與參考基因組比對，堿基錯誤率在1%以下（Q20）的占比95.75%～96.94%，平均為96.54%；錯誤率低于0.1%（Q30）的占比89.13%～92.00%，平均為90.97%。GC 含量32.98%～38.99%，平均為36.69%。與參考基因組的比對率為72.54%～95.26%，平均為81.87%，測序深度最高為9.69×，最低為6.97×，平均測序深度為8.42×，＞1×堿基覆蓋度平均為95.26%。上述結(jié)果表明，各樣本與參考基因組的相似度滿足重測序分析要求，可進(jìn)行后續(xù)分析。各樣本數(shù)據(jù)見表2。

表2 重測序數(shù)據(jù)概況Table 2 Overview of the whole genome resequencing data

2.2 SNP 類型與分布

采用GATK 進(jìn)行SNP 檢測，過濾后最終獲得7 492 966 個高質(zhì)量SNPs 用于后續(xù)分析。進(jìn)一步統(tǒng)計上述SNP的變異類型，發(fā)現(xiàn)在6種變異類型中，A/G 和C/T 變異類型最多，分別占比32.00%和31.96%（圖1A），即4 791 902 個SNPs 發(fā)生了堿基轉(zhuǎn)換（Ts），2 701 064 個SNPs 發(fā)生了堿基顛換（Tv），Ts/Tv=1.78。進(jìn)一步統(tǒng)計分析他們在基因組的分布情況，發(fā)現(xiàn)有57.8%的SNP 位于基因間區(qū)，20.9%的SNP 位于內(nèi)含子區(qū)，6.4%的SNP 位于外顯子區(qū)，基因上游區(qū)域和下游區(qū)域內(nèi)的SNP分別占比5.7%和5.0%（圖1B）。

圖1 泡桐屬11 個桐種變異位點及注釋信息的統(tǒng)計Fig.1 Statistics of location and annotation for the SNP and InDel variants of 11 Paulownia samples

2.3 SNP 的種間變異

不同品種泡桐個體中檢測到的總變異數(shù)各不相同，毛泡桐具有最多的變異數(shù)（4 900 517）（圖1C），和較低的堿基轉(zhuǎn)換/堿基顛換比值（TS/TV=1.77，在11 種泡桐中僅高于臺灣泡桐和川泡桐），表明毛泡桐的基因交換事件發(fā)生的較少。不同種泡桐非同義突變與同義突變的比值（Nonsyn/Syn）相差不大，白花泡桐最高，達(dá)到1.03，其他桐種為0.94～0.95（表3，圖1D）。

表3 不同桐種的變異情況Table 3 Summary of variants from all Paulownia samples

2.4 群體進(jìn)化樹與主成分分析

基于SNP 標(biāo)記，采用VCF2Dis 計算各桐種間的遺傳距離，結(jié)果見圖2。泡桐屬植物種間的遺傳距離為0.15～0.59，其中，毛泡桐與白花泡桐、臺灣泡桐與白花泡桐的遺傳距離分別為0.59 和0.58，屬于親緣關(guān)系較遠(yuǎn)的2 組桐種。宜昌泡桐與山明泡桐，鄂川泡桐與山明泡桐的遺傳距離均為0.15，屬于親緣關(guān)系較近的2 組桐種。根據(jù)種間遺傳距離，構(gòu)建泡桐屬植物的系統(tǒng)發(fā)育樹，結(jié)果表明11 個桐種可以被分為3 大支（圖3）。其中，臺灣泡桐和川泡桐為一支（Ⅰ），毛泡桐單獨為一支（Ⅱ），白花泡桐、楸葉泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、華東泡桐、建始泡桐和蘭考泡桐為一支（Ⅲ）。通過主成分分析方法再次分析11 個桐種間的親緣關(guān)系，根據(jù)主成分1 和主成分2 同樣將11 份材料分為3 個類群：毛泡桐為一個類群，臺灣泡桐和川泡桐為一個類群，其他桐種為一個類群（圖4），結(jié)果與系統(tǒng)發(fā)育分析結(jié)果高度一致。這表明，本研究鑒定到的SNPs 在泡桐屬植物親緣關(guān)系評價中可靠性較高。

圖2 基于SNP 計算的11 個桐種的遺傳距離Fig.2 Genetic distances of 11 Paulownia samples based on SNPs

圖3 基于SNP 構(gòu)建的11 個桐種的Neighbor Joining 進(jìn)化樹Fig.3 Neighbor Joining phylogenetic tree of 11 Paulownia samples based on SNPs

圖4 基于SNP 構(gòu)建的11 個桐種的主成分分析Fig.4 Principal component analysis of 11 Paulownia samples based on SNPs

2.5 群體遺傳結(jié)構(gòu)分析

群體遺傳結(jié)構(gòu)分析是遺傳關(guān)系分析的一種重要手段[18]。為了進(jìn)一步解析11 個桐種間的親緣關(guān)系，采用FastSTRUCTURE 軟件對他們進(jìn)行群體遺傳結(jié)構(gòu)分析。通過假定不同的祖先來源數(shù)目K，可以將不同群體之間是否含有相同祖先成分表示出來。在Structure 群體遺傳結(jié)構(gòu)圖中，每一列豎格代表一個樣本的遺傳背景，每種顏色的色塊代表估計的一個祖先，每種色塊占該豎格的比例代表該祖先對該樣品遺傳背景的貢獻(xiàn)比例[25]。為了確定合適的分群數(shù)量，假定祖群值K為1～9，將不同群體之間是否含有相同的祖先成分表示出來。fastSTRUCTURE 中chooseK 程序可根據(jù)種群結(jié)構(gòu)邊際似然近似值（Marginal likelihood）自動篩選最佳K值，結(jié)果顯示，3 是邊際似然值達(dá)到最大時的最小K值，即最合理的祖先種群數(shù)為3，表明11 個桐種可能來自于3 個祖先（圖5）。該結(jié)果與系統(tǒng)發(fā)育樹和主成分分析結(jié)果相互吻合，即來源相同的3 個類群可以獨立成群，類群Ⅰ包括毛泡桐；類群Ⅱ包括川泡桐和臺灣泡桐；類群Ⅲ包括白花泡桐、楸葉泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、華東泡桐、建始泡桐和蘭考泡桐。值得一提的是，類群Ⅲ中的華東泡桐顯示了來自類群Ⅱ的混合遺傳成分，說明類群Ⅱ和類群Ⅲ可能存在基因交流（圖6）。

圖5 基于fastSTRUCTURE 計算不同K 值時種群結(jié)構(gòu)的邊際似然值Fig.5 Marginal likelihood value of Paulownia samples with different K values calculated by fastSTRUCTURE

圖6 基于fastSTRUCTURE 推斷的11 個桐種的群體分層Fig.6 Population stratification based on fastSTRUCTURE analysis of 11 Paulownia samples

2.6 SNP 相關(guān)基因的功能分析

蛋白質(zhì)編碼序列的突變可能會影響基因功能和不同的生物學(xué)過程。采用變異密度法（每kb變異數(shù)）計算泡桐屬植物蛋白質(zhì)編碼基因的突變率（Mutation rates，MR）[26]，根據(jù)突變率分布的四分位距，這些基因被分為高、較高、較低和低突變率4 組。高突變率組包括4 597 個基因（MR ≥6.0），較高突變率組包括4 033 個基因（6.0 ＞MR ≥4.0），較低突變率組包括8 395 個基因（4.0 ＞MR ≥2.0），低突變率組包括9 878個基因（MR ＜2.0）。進(jìn)一步GO 富集分析表明，進(jìn)化速度不同的基因，其功能類別也存在差異。高突變率的基因被富集到細(xì)胞壁代謝途徑，包括木葡聚糖代謝調(diào)控、次生細(xì)胞壁和木栓素生物合成過程相關(guān)基因，如β-己糖胺酶、葡聚糖半乳糖基轉(zhuǎn)移酶、GMA12/MNN10 家族半乳糖基轉(zhuǎn)移酶和UDP-糖基轉(zhuǎn)移酶等33 個基因；花粉代謝途徑，包括花粉壁組件、花粉管和花藥發(fā)育相關(guān)基因，如ABORTED 轉(zhuǎn)錄因子、ABC 轉(zhuǎn)運蛋白和PHDfinger 等31 個基因；次生代謝物代謝通路，包括苯丙烷生物合成過程、次生代謝物生物合成過程和次生代謝過程相關(guān)基因，如苯丙氨酸解氨酶、GST 轉(zhuǎn)運蛋白和AP2 等35 個基因；形態(tài)建成相關(guān)通路，包括植物形態(tài)器官發(fā)生調(diào)控、幼苗發(fā)育調(diào)控和側(cè)根形成相關(guān)基因，如銨轉(zhuǎn)運蛋白、ABC-2型轉(zhuǎn)運蛋白和MYB 等53 個基因（圖7）。

圖7 高突變率基因集中所包含細(xì)胞壁合成、花粉、次生代謝和形態(tài)建成相關(guān)基因的GO 富集分析Fig.7 GO enrichment analysis for the genes involved in the pathways of cell way,pollen,secondary metabolite and morphogenesis with high mutation rates (MR)

3 討 論

我國泡桐種質(zhì)資源豐富，但很多材料缺乏可追溯的譜系，遺傳背景不明，給泡桐屬植物核心育種群體的建立、親本選配和種質(zhì)創(chuàng)新等帶來諸多不便。李芳東等[2]在對全國泡桐屬種質(zhì)資源調(diào)查收集過程中，從形態(tài)特征、生態(tài)習(xí)性、分布狀況和性狀穩(wěn)定性等方面對泡桐屬植物進(jìn)行了全面比較分析，提出泡桐屬主要包括11 個桐種的認(rèn)識（表1），推動了泡桐的遺傳演化及分類學(xué)研究。

本研究共鑒定到7 492 966 個高質(zhì)量SNPs，其中堿基轉(zhuǎn)換與顛換比為1.78，說明堿基轉(zhuǎn)換是泡桐屬植物主要的SNP 變異類型，該結(jié)果與茶樹Camelliasinensis(L.) Kuntze[27]、葡萄Vitis[28]和毛竹Phyllostachysedulis[29]中SNP 的鑒定結(jié)果基本一致。研究表明，堿基轉(zhuǎn)換的SNP 對蛋白質(zhì)氨基酸序列的改變小于堿基顛換的SNP[30]。本研究鑒定的大多數(shù)SNP 分布在基因組的基因間區(qū)和非編碼區(qū)，該結(jié)果與咖啡和茶樹中SNP 的注釋結(jié)果相似[31-32]。雖然這些SNP 可能不影響基因功能，但卻能夠體現(xiàn)植物群體間的遺傳和進(jìn)化關(guān)系[18]，作為遺傳標(biāo)記應(yīng)用于泡桐屬植物遺傳關(guān)系分析。

系統(tǒng)進(jìn)化樹、主成分分析以及遺傳結(jié)構(gòu)分析的結(jié)果高度一致，均支持3 個亞群的分類模式。其中，毛泡桐為類群Ⅰ；川泡桐和臺灣泡桐為類群Ⅱ；白花泡桐、楸葉泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、華東泡桐、建始泡桐和蘭考泡桐為類群Ⅲ，該結(jié)果與范永明[33]基于花序形狀、花蕾形狀和花色等22 個性狀采用數(shù)量分類學(xué)方法取得的結(jié)果基本一致。本研究還發(fā)現(xiàn)，地理位置相近的桐種并沒有因此聚到同一類群，說明3 個類群可能的基因來源與地理位置沒有明確的相關(guān)性[34]。

本研究中，高突變率基因（MR ≥6.0）被富集到細(xì)胞壁代謝、花粉代謝、次生代謝物代謝和形態(tài)建成等生物學(xué)過程，共同作用導(dǎo)致了不同桐種在遺傳、基因功能和表型上的多樣性。比如，木葡聚糖代謝調(diào)控、次生細(xì)胞壁和木栓素生物合成過程等相關(guān)基因影響細(xì)胞壁的代謝過程。樹干形成層向內(nèi)分化形成具有生活力的木質(zhì)部細(xì)胞，失活后殘留的細(xì)胞壁是影響木材品質(zhì)的重要成分[35]。因此，這些SNP 可能在木材形成過程中相關(guān)基因的功能調(diào)控中發(fā)揮作用，需要后續(xù)進(jìn)一步研究。再如，花粉壁組件、花粉管和花藥發(fā)育相關(guān)基因影響花粉的代謝過程，這些SNP 可能是泡桐屬在長期進(jìn)化過程中出現(xiàn)楸葉泡桐和山明泡桐2 個雄性不育種[1]的誘因之一。楸葉泡桐是泡桐屬11 個種中唯一兼具適宜密植和材質(zhì)優(yōu)良特性的桐種，是培育裝飾材的首選桐種[36]，雄性不育致使其只能用作母本進(jìn)行雜交育種，限制了雜種優(yōu)勢的利用。對此類關(guān)鍵SNP 位點的篩選和鑒定可為通過分子育種手段恢復(fù)楸葉楸葉泡桐的育性提供參考。

本研究以泡桐屬11 個桐種的典型單株為對象，利用WGRS 技術(shù)挖掘覆蓋全基因組的SNP 標(biāo)記，對解析泡桐屬的種間親緣關(guān)系及遺傳結(jié)構(gòu)具有重要意義，但也存在一些局限，下一步還需要增加種內(nèi)不同基因型或不同種源的試驗材料參與分析，開展更深入的理論研究，更有效地支撐泡桐種質(zhì)創(chuàng)新工作。

4 結(jié) 論

本研究首次報道了泡桐屬11 個桐種間的全基因組遺傳變異。通過全基因組重測序鑒定出7 492 966 個高質(zhì)量SNPs。這些SNP 不僅有助于進(jìn)行泡桐種質(zhì)資源遺傳變異分析和起源分化探索，同時還對鑒定具有表型效應(yīng)的等位基因變異具有重要意義。本研究基于這些SNP 對泡桐屬植物的親緣關(guān)系和遺傳結(jié)構(gòu)進(jìn)行了系統(tǒng)分析，將為泡桐種質(zhì)創(chuàng)新、挖掘優(yōu)良性狀相關(guān)基因進(jìn)而培育突破性良種奠定基礎(chǔ)。