基于簡化基因組測序的不同干形云南松遺傳分析*

吳治洋，縱 丹,3，甘沛華，劉子騰，朱梅彩，何承忠,3

(1.西南林業(yè)大學，云南省高校林木遺傳改良與繁育重點實驗室，云南 昆明 650224；2.西南林業(yè)大學，西南地區(qū)生物多樣性保育國家林業(yè)和草原局重點實驗室，云南 昆明 650224；3.西南林業(yè)大學，西南山地森林資源保育與利用教育部重點實驗室，云南 昆明 650224)

云南松(Pinus yunnanensis)是中國荒山荒地造林的主要樹種，也是西南地區(qū)重要的用材樹種，具有生長快、材質(zhì)較好、耐干旱和木材用途廣泛等優(yōu)良特點[1-3]。云南松以云貴高原為起源中心和分布中心，分布于N23°～30°、E96°～108°，其輻射分布范圍包括云南大部分地區(qū)、西藏東南部、四川西南部、廣西西部和貴州西部[4-6]。云南省的云南松林地面積達344.84 hm2，占林分總面積的52%，是主要的森林植被覆蓋類型之一[4,7]，在經(jīng)濟、社會和環(huán)境的可持續(xù)發(fā)展中具有其他樹種不可替代的重要作用[8]。隨著氣候變遷、人類干擾活動加劇及自然災害頻發(fā)等因素，云南松的生存狀況不容樂觀。當前云南松林分中80%以上是純林，林分組成極其簡單[9]。絕大多數(shù)云南松林分由低矮、彎曲和扭曲等個體構(gòu)成，呈現(xiàn)為低質(zhì)低效林，林分表現(xiàn)出嚴重的衰退現(xiàn)象，優(yōu)良遺傳資源正在急劇減少[1,4,9-10]。云南松干形扭曲形成的機制還存在較大爭議，屬于遺傳因素決定，還是環(huán)境因素引起尚不明確。

簡化基因組測序(reduced-representation genome sequencing，RRGS)技術(shù)基于二代高通量測序技術(shù)，通過限制性內(nèi)切酶對基因組進行打斷，經(jīng)過添加接頭和過濾等篩選基因組特異性長度區(qū)域進行測序，進而開發(fā)大量遺傳標記[11-12]。該技術(shù)大大降低了基因組的復雜性，同時具有不依賴參考基因組而獲得大量遺傳標記、操作簡單和成本較低等優(yōu)點。目前，RRGS 技術(shù)主要被用于分子標記開發(fā)、群體遺傳及系統(tǒng)發(fā)生學、高密度遺傳圖譜構(gòu)建和QTL 定位等[13-14]。在裸子植物中，胚乳來源于雌配子體，存在于受精作用前，由大孢子葉直接發(fā)育而來，為單倍體，具有染色體來源于母本的特點。基于此，本研究利用RRGS 技術(shù)中基于測序的基因分型 (genotyping by sequencing，GBS)雙酶切(ddGBS)對不同干形云南松胚乳和針葉進行測序分析，根據(jù)單倍體胚乳和二倍體針葉的組內(nèi)及組間序列比對，在全基因組范圍中挖掘不同干形云南松SNP 標記并分析其特征，揭示不同干形云南松遺傳變異特點，為進一步探究云南松干形變異機制提供科學依據(jù)。

1 材料與方法

1.1 試驗材料

從云南省楚雄彝族自治州永仁縣方山云南松純林居群中，選取樹齡相近的直干形云南松5 株和扭紋角＞30o的扭干形云南松6 株，分別采集針葉和球果，同類干形的樣株相隔50 m 以上。針葉樣品裝入凍存管后立即放入液氮中凍存；球果帶回實驗室，人工脫粒后去除外種皮和合子胚，僅保留胚乳，按單株將人工剝?nèi)〉呐呷檠b入凍存管，置于液氮中凍存。樣本分為4 組，即直干形云南松胚乳組(ES01、ES02、ES03、ES04 和ES05)、直干形云南松針葉組(NS01、NS02、NS03、NS04和NS05)、扭干形云南松胚乳組(ET01、ET02、ET03、ET04、ET05 和ET06)以及扭干形云南松針葉組(NT01、NT02、NT03、NT04 和NT05、NT06)。

1.2 DNA 的提取及質(zhì)量檢測

采用HiPureSFPlant DNA Kit 試劑盒分別提取針葉和胚乳的總DNA，DNA 的完整性、純度和質(zhì)量濃度采用1%的瓊脂糖凝膠電泳和核酸蛋白檢測儀共同進行檢測。合格的DNA 樣品稀釋至100 ng/μL，并置于-20 ℃冰箱保存，備用。

1.3 文庫構(gòu)建

ddGBS 文庫構(gòu)建及測序委托廣州基迪奧生物科技有限公司承擔完成。具體操作步驟：采用特定的2 種限制性內(nèi)切酶 (EcoR I+NlaIII)對質(zhì)檢合格的基因組DNA 進行消化；再進行末端修復，添加A-尾和Illumina 測序接頭；對300～400 bp的DNA 片段進行PCR 擴增富集；使用AMPure XP系統(tǒng)對PCR 產(chǎn)物進行純化，使用Agilent 2100 生物分析儀(Agilent，USA)對測序文庫進行檢測，并使用實時定量PCR 進行文庫定量，合格后采用Novaseq 6000 測序儀進行測序，測序策略為Illumina PE150。

1.4 高通量測序分析

Novaseq 6000 測序儀測序獲得的原始圖像數(shù)據(jù)經(jīng)Base Calling 轉(zhuǎn)化為序列數(shù)據(jù)，得到原始測序數(shù)據(jù)文件，結(jié)果以FASTQ 文件格式存儲，包含reads 的序列以及堿基的測序質(zhì)量。在測序數(shù)據(jù)下機之后，應用Fastp 軟件對下機數(shù)據(jù)進行質(zhì)量控制，過濾去除含adapter 的reads、含N 比例大于10%的reads 以及質(zhì)量值(Q)≤10且堿基數(shù)占整條read 的50%以上的低質(zhì)量reads，得到高質(zhì)量的數(shù)據(jù)(high quality clean reads)用于后續(xù)信息分析。

1.5 基因組SNP 檢測與開發(fā)

使用VAN DER AUWERA 等[15]的UnifiedGenotyper 模塊將處理好的比對文件進行多個樣本的變異檢測，以期獲得較高質(zhì)量的SNP 位點，對檢測所得的變異進行以下參數(shù)過濾：每15 bp 中存在超過3個SNP 位點、單位深度變異值(QD)＜4、比對質(zhì)量值離散程度(MQ)＜40、Fisher’s 檢測P值(FS)＞60 和基因型質(zhì)量值(GQ)＜20 的異常序列，過濾后再分別進行個體聚類以及群體聚類分析。

1.6 群體分析

根據(jù)所得的SNP 信息計算樣本之間的遺傳距離，對樣本進行聚類分析，從而推斷出樣本之間親緣關系遠近。利用treebest[16]鄰近法構(gòu)建進化樹。使用PLINK[17]和GCTA[18]進行主成分分析。采用admixture 軟件[19]進行群體結(jié)構(gòu)分析。使用perl 編程對云南松直、扭干形之間的期望雜合度(He)、觀測雜合度(Ho)、多態(tài)性信息含量(PIC)、觀測等位基因數(shù)(Na)、有效等位基因數(shù)(Ne)、基因多樣性指數(shù)(Nei)、Shannon 信息指數(shù)(I)、遺傳分化指數(shù)(Fst)和基因流(Nm)等指標進行計算。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)控和過濾結(jié)果

由于簡化基因組測序reads 為基因組DNA的酶切片段，其堿基組成分布檢查主要是檢測有無AT、GC 的分離現(xiàn)象?；诿盖衅蔚碾S機性和堿基互補配對原則，理論上每個測序循環(huán)的AT 和GC 含量相等，并在測序全過程中基本保持不變。本研究樣品ES01 的測序堿基組成分布如圖1 所示：AT 含量和GC 含量基本保持穩(wěn)定，模糊堿基N 所占比例較低，表明未知堿基數(shù)較少，且測序樣本受系統(tǒng)AT 偏好影響較小。結(jié)果表明：該樣品的文庫構(gòu)建質(zhì)量與測序質(zhì)量均可滿足后續(xù)相關分析。

圖1 樣品ES01 的堿基組成分布圖Fig.1 Distribution of base composition of sample ES01

對質(zhì)控后數(shù)據(jù)進行統(tǒng)計，共獲得高質(zhì)量序列數(shù)490 萬個，4 組分析樣本過濾后得到的高質(zhì)量測序數(shù)據(jù)在14.30～17.15 G 之間，GC 含量在44.61%～50.60%之間，堿基質(zhì)量Q30 比例在93.65%～95.19%之間(表1)。由此可知：所測序列的Q30 數(shù)據(jù)較高，表明堿基測序出錯率很低，GC 分布合理，數(shù)據(jù)量達到分析要求，建庫測序成功。

表1 云南松針葉和胚乳簡化基因組DNA 測序數(shù)據(jù)Tab.1 ddGBS sequencing of endosperm and needle of Pinus yunnanensis

2.2 SNP 和InDel 位點檢測

變異位點檢測結(jié)果顯示：共檢測出772 240個變異，其中SNP 變異有762 699個，InDel 變異9 541個。在SNP 變異中，顛換(38.28%)小于轉(zhuǎn)換(61.72%)；其中，顛換類型中，置換比例從小到大依次為T→A (2.40%)=A→T (2.40%)＜ T→G (3.80%)＜ A→C (3.88%)＜ G→T (4.16%)＜ C→A(4.19%)＜ C→G (8.61%)＜ G→C (8.83%)，而轉(zhuǎn)換類型中，置換比例從小到大依次為A→G (14.72%)＜T→C (14.83%)＜ G→A (15.96%)＜ C→T (16.21%)。在9 541個InDels 中，插入型變異共4 597個(48.18%)，缺失型變異為4 944個(51.82%)。In-Del 的長度從1～53 bp 均有分布，主要以小片段(1～10 bp)分布為主。其中，單堿基InDels 的數(shù)量最多，達到5 207個；2 堿基InDels 次之，有1 211個；3 堿基InDels 變異為480個(圖2)。

圖2 InDel 變異類型Fig.2 Variation types of InDel

2.3 遺傳變異分析

基于得到的762 699個SNP 標記，對云南松4個樣本組分別進行遺傳變異分析。結(jié)果(表2)顯示：云南松直干型胚乳組(ES)、直干型針葉組(NS)、扭干型胚乳組(ET)和扭干型針葉組(NT)的期望雜合度(He)為0.203 9～0.229 9、觀測雜合度(Ho)為0.100 2～0.145 7、多態(tài)性信息含量(PIC)為0.163 3～0.183 7、觀測等位基因數(shù)(Na)為0.340 5～0.442 3、有效等位基因數(shù)(Ne)為1.347 2～1.388 4，基因多樣性指數(shù)(Nei)為0.234 5～0.262 5、Shannon’s 信息指數(shù)(I)為0.304 3～0.344 8。由表2還可知：無論采用針葉材料(雙親型，二倍體)還是種子胚乳材料(單親型，單倍體)，云南松遺傳變異參數(shù)較為接近，表明云南松個體具有較高的雜合性。

表2 不同組云南松遺傳變異指標Tab.2 Genetic variation indexes of different groups of P.yunnanensis

4個樣本組之間的遺傳分化系數(shù)和基因流分析結(jié)果(表3)顯示：組間的Fst值在0.023 2～0.143 7 之間，而不同干形針葉比較組(NS 與NT)和不同干形胚乳比較組(ES 與ET)的Fst值分別為0.107 1 和0.143 7，表明直、扭干形為中等程度遺傳分化(0.05＜Fst＜0.15)。然而，不同組間基因交流豐富，Nm值均大于1，介于1.489 3～ 10.519 0之間。其中，直干形之間(Nm=6.363 4)及扭干形之間(Nm=10.519 0)的基因流值遠高于直干形和扭干形云南松之間的基因流值(Nm針葉=2.083 9，Nm胚乳=1.489 3)。

表3 云南松4個樣本組間遺傳分化系數(shù)和基因流Tab.3 Fst and Nm between four sample groups of P.yunnanensis

2.4 系統(tǒng)發(fā)育分析

采用treebest 對得到的762 699個SNP 位點構(gòu)建進化樹，把具有最低交叉驗證錯誤率的分群數(shù)(K值)定義為最優(yōu)分群數(shù)。當K值取值范圍為1～9 時，隨著K值增大，交叉驗證錯誤率呈現(xiàn)下降趨勢(圖3a)。當K=2 時，群體分為2個類群，但直干形和扭干形樣本區(qū)分不明顯；當K=9 時，具有最低交叉驗證錯誤率，樣本分為9個類群(圖3b)。由此可見，本研究中采集的22個云南松樣本間遺傳結(jié)構(gòu)比較復雜。

圖3 基于SNP 分析群體結(jié)構(gòu)Fig.3 Analysis of population structure by SNPs

系統(tǒng)進化樹(圖4)顯示：22個樣本聚為2 大類，即針葉型類群和胚乳型類群。進一步分析發(fā)現(xiàn)：在針葉型類群中，除NT01 被聚類到胚乳組中，其余不同干形樣本沒有明顯的分群現(xiàn)象；而在胚乳類群中，除ET06 樣本外，不同干形的樣本具有較好的分群現(xiàn)象?？梢?，采用胚乳材料進行不同干形云南松的分型比較可行。

圖4 基于SNP 變異位點構(gòu)建的系統(tǒng)進化樹Fig.4 Phylogenetic tree based on SNPs

進一步以直干形樣本SNP 位點的相同堿基為參照，篩選扭干形云南松樣本的SNP 變異位點，共挑選出27 條序列中的39個SNP 位點。采用鄰近法對39個SNP 位點繪制系統(tǒng)進化樹，結(jié)果(圖5)顯示：不同干形的樣本能被明確地區(qū)分，說明這些位點與云南松干形變異的發(fā)育調(diào)控相關。對這27 條序列進行數(shù)據(jù)庫比對注釋，僅有5 條序列被NR 和GO 數(shù)據(jù)庫注釋(圖6)，注釋結(jié)果依次為多向耐藥性蛋白(pleiotropic drug resistance protein 3，PDR3)、乙烯響應轉(zhuǎn)錄因子(ethylene-responsive transcription factor，ERF017 和RAP2-9)及MYB 相關蛋白340 (MYB-related protein 340)，表明這些因子可能參與了云南松不同干形發(fā)育的調(diào)控作用。

圖5 基于篩選的39個SNP 位點構(gòu)建的系統(tǒng)進化樹Fig.5 Phylogenetic tree based on selected 39 SNPs

圖6 含有SNP 位點的5個序列片段注釋信息Fig.6 Annotation of five SNP location fragments

2.5 主成分分析

基于762 699個SNP 位點對云南松22個樣本進行主成分分析，結(jié)果(圖7)顯示：不同干形云南松個體交錯分布，沒有產(chǎn)生明顯的分群現(xiàn)象，表明云南松個體具有較高的雜合性，從而導致云南松干形變異的復雜性。

圖7 22個云南松樣本主成分分析Fig.7 Principal component analysis of 22 Pinus yunnanensis samples

3 討論

簡化基因組是基于二代高通量測序技術(shù)通過特定的限制性內(nèi)切酶對基因組DNA 進行隨機切割，降低基因組復雜性，進而展現(xiàn)基因組序列特征[20-22]，其具有高通量、高準確性、高特異性和標記數(shù)量多等優(yōu)點[23]，尤其是適用于缺少相關參考基因組的物種，能有效簡化基因組的復雜性。當前簡化基因組技術(shù)廣泛地應用于植物分子育種和遺傳多樣性研究[24]。段義忠等[25]基于SLAF-seq技術(shù)獲得沙冬青56 295個多態(tài)性SLAF 標簽，并對其遺傳變異進行了分析；張序等[21]采用簡化基因組技術(shù)對寬杯杜鵑進行單核苷酸多態(tài)性位點挖掘和遺傳多樣性分析，共開發(fā)出103 133個高質(zhì)量SNP 位點，并對寬杯杜鵑的遺傳多樣性進行了分析。本研究分別以單倍型胚乳和二倍型針葉為材料，采用簡化基因組技術(shù)開發(fā)出762 699個云南松高質(zhì)量SNP 位點，表明云南松個體差異極豐富、遺傳背景極復雜。

豐富的遺傳多樣性是一個物種應對環(huán)境變化的基礎[26]，觀測雜合度、期望雜合度和多態(tài)性信息含量等參數(shù)被用來表征群體中遺傳多樣性的大小[27]。本研究共檢測出SNP 位點762 699個和InDel 位點9 541個。云南松直干型胚乳組、云南松直干型針葉組、云南松扭干型胚乳組和云南松扭干型針葉組的期望雜合度介于0.203 9～0.229 9、觀測雜合度介于0.100 2～0.145 7、多態(tài)性信息含量介于0.163 3～0.183 7，表明該研究中云南松樣本的遺傳多樣性較低，這是由于本研究旨在初步揭示云南松基于簡化基因組的SNP 位點分布特點，采用樣本數(shù)較少。在云南松不同干形針葉組和胚乳組中，F(xiàn)st值分別為0.107 1 和0.143 7 (0.05＜Fst＜0.15)，表明直干形云南松和扭干形云南松之間存在中等程度遺傳分化，該結(jié)果與遺傳變異分析、系統(tǒng)發(fā)育分析及主成分分析結(jié)果基本一致。這一結(jié)果也與采用SRAP[28]、AFLP[29]和SSR[30-31]標記的結(jié)果基本吻合。同時，不同干形云南松之間的基因交流值均大于1，表明不同干形云南松之間不存在生殖隔離或基因交流障礙，扭干形云南松僅為云南松的一種類型。

當前云南松林中存在廣泛且豐富的干形扭曲現(xiàn)象，其成因存在一定的爭議。樹木干形的生長發(fā)育與其分布的生境極其相關，干形的多樣性極大地反映出林木對周邊環(huán)境的適應性，同時也是進化的基礎[1-2,7,32]。陳守常等[32]對四川省集中分布的云南松居群進行了長期觀測調(diào)查，認為干形扭曲是云南松所具有的一種特性，林分條件(水分、土壤、年齡和疏密度)尤其疏密度是樹木干形扭曲的決定性因素，而風在一定程度上也是促進和強化扭曲的外在因子。姜漢僑[33]結(jié)合云南松干形調(diào)查資料提出：樹木干形扭曲是云南松的一種原有性狀，具有遺傳性，并認為在適宜的生境條件下，干形的扭曲度可以在一定程度上減小，而且干形的扭曲度變化是一種數(shù)量性狀。SHELBORUNE 等[34]對火炬松進行了深入研究，認為樹干通直度可遺傳，直干型樹木繁殖所得后代大部分也是樹干通直，其遺傳力為0.39。周蛟等[35]對云南松的遺傳改良試驗表明：優(yōu)良林分速生林生長量遺傳力為0.76，遺傳增益為0.36，其子代種子可以改變云南松林分中樹干彎曲及木材紋理扭曲的不良特性。由于云南松為雌雄同株植物，存在自花授粉和異花授粉，從而導致云南松林的個體之間存在較大差異[36]。同時，云南松花粉在風力作用下的水平傳播距離可達數(shù)百公里，同時在垂直面上還有大量擴散[37]。當直干形和扭干形花粉或種子通過風落入同一林分時，就可能出現(xiàn)同一林分中直、扭干形共存現(xiàn)象。本研究材料采集于直干形和扭干形云南松共存的同一林分內(nèi)，且部分不同干形云南松樣本相鄰分布。基于針葉材料開發(fā)獲得的SNP 位點不能有效地區(qū)分直、扭干形云南松，但基于胚乳材料開發(fā)得到的SNP 位點能夠較好地將2 種干形云南松分型。此外，對篩選出包含39個SNP 位點的27 條序列注釋表明：雖然僅有5 條序列被注釋，但注釋預測的乙烯響應轉(zhuǎn)錄因子ERF017 和RAP2 均屬于AP2/ERF植物轉(zhuǎn)錄因子家族，可通過響應乙烯、細胞分裂素、生長素及脫落酸和油菜素內(nèi)酯等內(nèi)源激素直接或間接調(diào)控植物生長發(fā)育及抗逆性[38-39]。基于此，認為云南松干形扭曲變異主要受遺傳因素調(diào)控，但環(huán)境因素的協(xié)同作用增強了該特征的顯著性。

4 結(jié)論

直、扭2 種干形云南松之間存在較為頻繁的基因交流，但不同干形云南松間存在中度遺傳分化，且個體差異極豐富、遺傳背景極復雜?；卺樔~材料開發(fā)得到的SNP 位點不能有效區(qū)分直、扭干形云南松，但基于胚乳材料開發(fā)得到的SNP 位點能夠較好地將2 種干形云南松分型；對篩選出的包含39個SNP 位點的27 條序列注釋分析表明：僅有5 條序列被注釋到多向耐藥性蛋白(PDR3)、乙烯響應轉(zhuǎn)錄因子 (ERF017，RAP2)及MYB 相關蛋白340 (MYB-related protein 340)，可能在云南松直干形和扭干形的發(fā)育中發(fā)揮著調(diào)控作用，故認為云南松干形扭曲變異主要受遺傳因素調(diào)控，但環(huán)境因素的協(xié)同作用增強了該特征的顯著性。