吳志銘，馮源恒，楊章旗

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

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

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

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

1 材料與方法

1．1 材料

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

1．2 RNA提取與檢測(cè)

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

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

1．3 SNP位點(diǎn)挖掘

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

1．4 SNP位點(diǎn)Unigene的GO功能分析

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

1．5 SNP位點(diǎn)Unigene的KEGG代謝通路分析

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

2 結(jié)果與分析

2．1 馬尾松SNP位點(diǎn)的挖掘

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

2．2 馬尾松SNP位點(diǎn)所在Unigene的GO分類

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

為了解篩選出的含有SNP 的Unigene 的功能，對(duì)GO 注釋結(jié)果進(jìn)行進(jìn)一步分類。這些Unigene 被注釋到3大類41個(gè)功能區(qū)（圖3）。其中93條Unigene 參與生物過(guò)程（Biological process）的18 個(gè)功能區(qū)，49條Unigene 參與細(xì)胞成分（Cellular component）的14 個(gè)功能區(qū)，92 條參與分子功能（Molecular function）的9 個(gè)功能區(qū)。生物過(guò)程中參與代謝過(guò)程（metabolic process，73 條）、細(xì)胞過(guò)程（cellular process，68 條）和單一生物過(guò)程（single-organism process,64條）的基因最多，參與生物過(guò)程的正、負(fù)調(diào)節(jié)(positive regulation and negative regulation of biological process)和免疫系統(tǒng)過(guò)程(immune system process)的基因最少，均只有1 條；細(xì)胞成分中參與細(xì)胞（cell，25 條）、細(xì)胞組分（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ò)程相關(guān)。

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

2．3 馬尾松SNP 位點(diǎn)Unigene 的KEGG 代謝通路分析

對(duì)轉(zhuǎn)錄組數(shù)據(jù)中125 條含SNP 的Unigene 的KEGG注釋結(jié)果進(jìn)行處理，發(fā)現(xiàn)57條含K編號(hào)，已知KEGG 功能的基因有29 條被注釋到5 大類13 個(gè)通路中（圖4）。其中新陳代謝（Metabolism）有21 條Unigene，環(huán)境信息處理（Environmental information processing）有4 條，組織系統(tǒng)（Organismal systems）有2 條，遺傳信息處理（Genetic information processing）和細(xì)胞過(guò)程（Cellular processes）各1條。新陳代謝的Unigene 可分為9 個(gè)亞類，氨基酸代謝類（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位點(diǎn)所在Unigene的KEGG的代謝通路Fig.4 Metabolic pathway of KEGG in Unigene at P.massoniana SNP loci

3 討 論

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

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