吳志銘，馮源恒，楊章旗

（1.廣西師范大學生命科學學院，廣西桂林 541006；2.廣西壯族自治區(qū)林業(yè)科學研究院廣西馬尾松工程技術(shù)研究中心 廣西優(yōu)良用材林資源培育重點實驗室，廣西南寧 530002）

馬尾松（Pinus massoniana）是松科（Pinaceae）松屬常綠針葉樹種，分布于秦嶺淮河以南和云貴高原以東等17 個省、自治區(qū)和直轄市[1]。馬尾松材脂兼用，木材可用于建筑、造紙和木纖維工業(yè)用材等領(lǐng)域；松脂加工產(chǎn)品是重要的工業(yè)原料，可應(yīng)用于醫(yī)藥、油漆和粘合劑等。中國70%以上的松脂產(chǎn)量來自馬尾松[2-4]。據(jù)第六次至第八次全國森林資源連續(xù)清查數(shù)據(jù)統(tǒng)計，馬尾松的森林面積不斷減少，加上不斷采割導(dǎo)致采脂樹木產(chǎn)脂能力降低及樹木死亡，致使馬尾松原脂產(chǎn)量出現(xiàn)萎縮。加快高產(chǎn)脂馬尾松選育已成為其遺傳改良的重要內(nèi)容。傳統(tǒng)選育通過連續(xù)幾年對處盛產(chǎn)期的馬尾松進行采脂測定，對其產(chǎn)脂能力做出評價，選出高產(chǎn)脂單株。其成本較高、耗時長、范圍窄且進展慢，開展分子輔助育種，能有效地縮短馬尾松育種年限[5]。

單核苷酸多態(tài)性（SNPs）是指單核苷酸在基因組水平上的突變引起的DNA 序列間的多態(tài)性，其突變包括單堿基的轉(zhuǎn)換、顛換和插入缺失[6]。SNP 因其在基因組中具有分布廣、多樣性高和易于分型等特點, 成為基因分型最理想的分子標記[7]。表達序列標簽（EST）是源于轉(zhuǎn)錄表達的特異功能基因的cDNA 片段。利用獲得的EST 序列進行SNP 標記開發(fā)，對未進行全基因組測序的動植物個體具有重要意義[8]。

本研究對馬尾松二代測序轉(zhuǎn)錄組數(shù)據(jù)中的SNP位點進行挖掘，并分析其功能和通路，為馬尾松產(chǎn)脂性狀關(guān)聯(lián)分析的開展提供可用的SNP標記。

1 材料與方法

1．1 材料

基于連續(xù)3年在南寧市林業(yè)科學研究所馬尾松高產(chǎn)脂種質(zhì)資源庫的采脂試驗結(jié)果，采集高產(chǎn)脂無性系桂GZ080B（15．96 g/10 cm）和普通產(chǎn)脂無性系桂GZ078B（8．18 g/10 cm）的3 個組織（頂芽、針葉和韌皮部）材料，迅速放入液氮中冷凍，送至北京諾禾致源科技股份有限公司進行高通量測序（Illumina HiSeqTM 2000/MiSeqTM），獲得轉(zhuǎn)錄組數(shù)據(jù)。

1．2 RNA提取與檢測

通過試劑盒法提取GZ080B 和GZ078B 的頂芽、針葉和韌皮部的RNA，分別使用1%瓊脂糖凝膠電泳監(jiān)測RNA 降解和污染，NanoPhotometer 分光光度計（IMPLEN，CA，USA）檢測RNA 的純度，Qubit 2．0 Flurometer（Life Technologies，CA，USA）中的Qubit RNA 分析試劑盒測量RNA 濃度，Agilent Bioanalyzer 2100 系統(tǒng)（Agilent Technologies，CA，USA）的RNA Nano 6000分析試劑盒評估RNA完整性（圖1）。

圖1 RNA凝膠電泳圖Fig.1 RNA gel electrophoresis

1．3 SNP位點挖掘

從高產(chǎn)脂無性系和普通產(chǎn)脂無性系的頂芽、針葉和韌皮部3 個比較組成的維恩圖中得到2 329 個差異表達基因，該圖可直觀展現(xiàn)各種組合間的差異表達基因數(shù)量。為篩選出含SNP位點的Unigene，從轉(zhuǎn)錄組數(shù)據(jù)中找到各比較組中的差異表達基因統(tǒng)計表（表格有7大數(shù)據(jù)庫注釋），根據(jù)NR注釋結(jié)果挑選已知功能的Gene ID，通過Novofinder（北京諾禾致源科技股份有限公司測序結(jié)果自帶）輸入Gene ID搜索含SNP 位點的Gene ID 的Unigene，并對其突變類型數(shù)量及所在密碼子位置進行統(tǒng)計，最后根據(jù)含SNP 位點的Gene ID 搜索其Unigene 的GO 功能注釋和KEGG 通路注釋結(jié)果。從中選出3 個比較組中共有的Unigene、只存在于針葉和韌皮部的Unigene 和韌皮部獨有的Unigene進行分析。

1．4 SNP位點Unigene的GO功能分析

根據(jù)含SNP 位點的Unigene 的GO 注釋結(jié)果，經(jīng)過數(shù)據(jù)處理，通過OmicShare 在線軟件的動態(tài)GO 富集分析（https://www．omicshare．com/tools/home/report/goenrich．html），把含SNP 位點的Unigene 的GO 注釋結(jié)果進行功能分類，最后對GO功能進行統(tǒng)計分析。

1．5 SNP位點Unigene的KEGG代謝通路分析

根據(jù)含SNP 位點的Unigene 在KEGG 數(shù)據(jù)庫的注釋結(jié)果，剔除沒有K 編號的Unigene，經(jīng)過數(shù)據(jù)處理，通過OmicShare 在線軟件的動態(tài)KEGG 富集分析（https://www．omicshare．com/tools/home/report/koenrich．html），對含編號的Unigene進行統(tǒng)計分析。

2 結(jié)果與分析

2．1 馬尾松SNP位點的挖掘

將高產(chǎn)脂無性系與普通產(chǎn)脂無性系的韌皮部、針葉與頂芽的轉(zhuǎn)錄組數(shù)據(jù)進行兩兩對比（圖2）。在發(fā)現(xiàn)的2 329 個差異表達基因中進行兩次篩選，第1次根據(jù)NR 數(shù)據(jù)庫的注釋結(jié)果篩選出366 條Unigene，第2 次篩選出含SNP 位點的Unigene 125 條，共656 個SNP 位點。對656 個SNP 位點的突變類型進行統(tǒng)計，發(fā)現(xiàn)轉(zhuǎn)換類型有4 種，顛換類型有8 種。發(fā)生轉(zhuǎn)換突變頻率較高，其中T/C 和C/T 轉(zhuǎn)換占總SNP位點的33．54%，A/G 和G/A 占28．05%；發(fā)生顛換類型的各種突變頻率較低，分別為11．59%（C、G）、9．60%（G、T）、8．69%（A、T）和8．53%（A、C）。對SNP位點所在密碼子位置進行統(tǒng)計時，發(fā)現(xiàn)只有47．99%的SNP 位點在密碼子上，在第一位置、第二位置和第三位置發(fā)生的突變比例分別為2∶1∶2。為進一步了解SNP 位點的信息，分析每個個體該位點的基因型和突變后的基因型，根據(jù)支持該位點的reads個數(shù)和GATK3 軟件得到的該位點的基因型，若/兩邊堿基相同，則為純合位點，若不同，則為雜合位點；基因型在不同部位中相同但在產(chǎn)脂能力不同的馬尾松中不同。根據(jù)統(tǒng)計結(jié)果發(fā)現(xiàn)，純合突變32 個，雜合突變121個。

2．2 馬尾松SNP位點所在Unigene的GO分類

圖2 差異表達基因Venn圖Fig.2 Venn map of differentially expressed genes

為了解篩選出的含有SNP 的Unigene 的功能，對GO 注釋結(jié)果進行進一步分類。這些Unigene 被注釋到3大類41個功能區(qū)（圖3）。其中93條Unigene 參與生物過程（Biological process）的18 個功能區(qū)，49條Unigene 參與細胞成分（Cellular component）的14 個功能區(qū)，92 條參與分子功能（Molecular function）的9 個功能區(qū)。生物過程中參與代謝過程（metabolic process，73 條）、細胞過程（cellular process，68 條）和單一生物過程（single-organism process,64條）的基因最多，參與生物過程的正、負調(diào)節(jié)(positive regulation and negative regulation of biological process)和免疫系統(tǒng)過程(immune system process)的基因最少，均只有1 條；細胞成分中參與細胞（cell，25 條）、細胞組分（cell part，25 條）和膜（membrane，24 條）的基因最多，其次是大分子復(fù)合物（macromolecular complex，18條）；分子功能中參與催化活性（catalytic activity，71 條）和結(jié)合活性（binding，62 條）的基因最多，參與抗氧化活性（antioxidant activity，1 條）和分子功能調(diào)節(jié)器（molecular function regulator，1 條）的基因最少。含有SNP 的Unigene 主要與馬尾松的代謝過程相關(guān)。

圖3 馬尾松SNP位點所在Unigene的GO分類Fig.3 GO classification of Unigene SNP loci of P.massoniana

2．3 馬尾松SNP 位點Unigene 的KEGG 代謝通路分析

對轉(zhuǎn)錄組數(shù)據(jù)中125 條含SNP 的Unigene 的KEGG注釋結(jié)果進行處理，發(fā)現(xiàn)57條含K編號，已知KEGG 功能的基因有29 條被注釋到5 大類13 個通路中（圖4）。其中新陳代謝（Metabolism）有21 條Unigene，環(huán)境信息處理（Environmental information processing）有4 條，組織系統(tǒng)（Organismal systems）有2 條，遺傳信息處理（Genetic information processing）和細胞過程（Cellular processes）各1條。新陳代謝的Unigene 可分為9 個亞類，氨基酸代謝類（Amino acid metabolism）最多，有8 條；其次是碳水化合物代謝類（Carbohydrate metabolism）和其他次生代謝產(chǎn)物合成（Biosynthesis of other secondary metabolites），分別有7 條和6 條；涉及萜類和聚酮化合物代謝（Metabolism of terpenoids and polyketides）的有4 條。結(jié)果表明，KEGG 代謝通路分析與GO 分類得出的結(jié)果均與代謝相關(guān)。

圖4 馬尾松SNP位點所在Unigene的KEGG的代謝通路Fig.4 Metabolic pathway of KEGG in Unigene at P.massoniana SNP loci

3 討 論

本研究根據(jù)馬尾松高產(chǎn)脂無性系桂GZ080B和普通產(chǎn)脂無性系桂GZ078B的轉(zhuǎn)錄組測序結(jié)果，在差異表達分析結(jié)果中根據(jù)NR數(shù)據(jù)庫注釋和含有SNP位點兩大特點，從2 329個Unigene中篩選出374條Unigene，共2 192個SNP位點，根據(jù)試驗要求從中選出3個比較組中共有的基因、只存在于針葉和韌皮部的基因和韌皮部獨有的基因。在對突變的類型進行統(tǒng)計時，發(fā)現(xiàn)C轉(zhuǎn)換為T的頻率最高（33．54%），原因是CG中的C常為甲基化狀態(tài)，自發(fā)脫氨后成為胸腺嘧啶。

綜合GO分類和KEGG代謝通路分析，從馬尾松轉(zhuǎn)錄組中篩選出的含SNP 的基因主要是與生物體代謝和分子功能相關(guān)，與萜類化合物和聚酮類化合物合成相關(guān)的基因較少。SNP標記位點源自編碼序列，通過EST 數(shù)據(jù)庫可以直接開發(fā)出與功能基因相關(guān)的SNP 標記，為進一步的功能基因研究提供依據(jù)[9]。本研究可為馬尾松產(chǎn)脂相關(guān)基因與標記開發(fā)等研究提供參考。