蔡年輝 鄧麗麗 許玉蘭 徐 楊 周 麗 王大瑋 田 斌 何承忠 段安安

(1.西南林業(yè)大學(xué)云南省高校林木遺傳改良與繁育重點(diǎn)實(shí)驗(yàn)室，昆明 650224； 2.西南林業(yè)大學(xué)國家林業(yè)局西南地區(qū)生物多樣性保育重點(diǎn)實(shí)驗(yàn)室，昆明 650224； 3.西南林業(yè)大學(xué)西南山地森林保育與利用省部共建教育部重點(diǎn)實(shí)驗(yàn)室，昆明 650224)

基于高通量測序的云南松轉(zhuǎn)錄組分析

蔡年輝1,2,3鄧麗麗1,2,3許玉蘭1,2,3徐 楊1,2,3周 麗1,2,3王大瑋1,2,3田 斌1,2,3何承忠1,2,3段安安1,2,3

(1.西南林業(yè)大學(xué)云南省高校林木遺傳改良與繁育重點(diǎn)實(shí)驗(yàn)室，昆明 650224；2.西南林業(yè)大學(xué)國家林業(yè)局西南地區(qū)生物多樣性保育重點(diǎn)實(shí)驗(yàn)室，昆明 650224；3.西南林業(yè)大學(xué)西南山地森林保育與利用省部共建教育部重點(diǎn)實(shí)驗(yàn)室，昆明 650224)

采用新一代高通量測序技術(shù)平臺(tái)Illumina Hiseq 2 000對(duì)云南松轉(zhuǎn)錄組測序，得到的數(shù)據(jù)進(jìn)行de novo組裝，獲得80 000條Unigenes，N50為1 881 nt、平均890 nt。與公共數(shù)據(jù)庫進(jìn)行比對(duì)，注釋到NR、NT、Swiss-Prot數(shù)據(jù)庫的Unigenes分別為43 434、46 415、29 418條。將Unigenes與COG數(shù)據(jù)庫比對(duì)，有14 792條Unigenes成功注釋，根據(jù)功能大致分成25類；與GO數(shù)據(jù)庫比對(duì)，有26 743條Unigenes獲得注釋，按功能分為細(xì)胞組分、分子功能和生物過程3大類55亞類，其中參與的生物過程較多；以KEGG數(shù)據(jù)庫參考，有25 873條Unigenes參與128條代謝途徑分支，以代謝相關(guān)的通路較為集中，并找到與木質(zhì)素合成關(guān)鍵酶的Unigenes。這些研究極大地?cái)U(kuò)充了云南松的基因資源，將有助于云南松基因的發(fā)掘與利用、分子標(biāo)記的開發(fā)及其種質(zhì)資源遺傳改良的研究等。

云南松；轉(zhuǎn)錄組；Illumina高通量測序

云南松(PinusyunnanensisFranch.)是分布于我國亞熱帶山地暖性氣候條件下的重要樹種，是云貴高原的主要鄉(xiāng)土樹種和西南地區(qū)特有的森林植被類型[1～3]，水平分布于東經(jīng)96°～108°、北緯23°～30°[1,4～5]，在云南省境內(nèi)垂直分布于海拔710～3 320 m，大多分布于1 500～2 500 m[6]。具有適應(yīng)性強(qiáng)、耐干旱瘠薄、喜光、木材用途廣泛等特點(diǎn)，是分布區(qū)域內(nèi)荒山造林的先鋒樹種和治理水土流失的重要樹種，在分布區(qū)域林業(yè)生產(chǎn)和生態(tài)經(jīng)濟(jì)建設(shè)中占有舉足輕重的地位，分別占云南省林分總面積和木材蓄積量的19.63%和14.28%[6～8]。目前云南松的研究以資源狀況、遺傳多樣性、培育等方面報(bào)道較多，鄧喜成等[6]通過對(duì)森林資源監(jiān)測，比較云南松林的資源動(dòng)態(tài)，表明云南松林可用資源數(shù)量仍在快速減少，材種結(jié)構(gòu)低質(zhì)化傾向加劇，迫切需要加強(qiáng)云南松林的撫育管理，科學(xué)開展森林經(jīng)營活動(dòng)；Xiao等[9]等提出云南松主要分布的區(qū)域，海拔2 500 m以下人類活動(dòng)對(duì)資源的影響較大；Wu等[10～11]對(duì)云南松林的土壤因子分析，為云南松資源的保護(hù)與改良提供借鑒；Wang等[12]以云南松分布區(qū)的群體為對(duì)象，探討云南松的遺傳分化及其可能的成因，為開展云南松的基礎(chǔ)研究提供指導(dǎo)，而轉(zhuǎn)錄組和基因組信息比較缺乏，對(duì)云南松分子標(biāo)記開發(fā)也相對(duì)滯后[13]。

對(duì)于缺乏基因組信息的物種，采用轉(zhuǎn)錄組測序可獲得大量的數(shù)據(jù)信息，有利于挖掘重要功能基因，是開展植物優(yōu)良性狀研究的重要手段[14～17]。隨著高通量測序技術(shù)的不斷發(fā)展，降低了測序的成本，同時(shí)獲得豐富的數(shù)據(jù)，利于轉(zhuǎn)錄組的分析，有助于基因結(jié)構(gòu)及其功能的研究，同時(shí)也為大量分子標(biāo)記的開發(fā)奠定基礎(chǔ)[18～22]。鑒于此，本研究采用Illumina HiSeq 2000高通量測序技術(shù)對(duì)云南松轉(zhuǎn)錄組測序分析，獲得大量的數(shù)據(jù)，對(duì)其進(jìn)行拼接與組裝，建立云南松轉(zhuǎn)錄組數(shù)據(jù)庫，結(jié)合生物信息學(xué)的方法對(duì)轉(zhuǎn)錄組數(shù)據(jù)庫序列進(jìn)行功能注釋、功能分類和代謝通路分析，以期為云南松基因組水平上的研究奠定基礎(chǔ)，同時(shí)也可為云南松EST-SSR分子標(biāo)記的開發(fā)和功能基因挖掘提供數(shù)據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)材料為永仁白馬河云南松母樹林采種育苗的5年生苗木，于2014年10月單株取樣，采集當(dāng)年生幼嫩的頂梢，然后迅速放入干冰中保存。

1.2 轉(zhuǎn)錄組測序

對(duì)采集的材料進(jìn)行RNA提取，構(gòu)建cDNA文庫，采用Illumina Hiseq 2000平臺(tái)測序。

1.3 序列組裝

轉(zhuǎn)錄組測序后統(tǒng)計(jì)原始數(shù)據(jù)的數(shù)量及其長度。對(duì)原始數(shù)據(jù)經(jīng)去除測序接頭、重復(fù)冗余序列及低質(zhì)量的序列數(shù)據(jù)等，獲得clean reads，統(tǒng)計(jì)clean reads的數(shù)量、總長度、Q20、N%、GC%等，Q20表示過濾后質(zhì)量不低于20的堿基的比例，N表示過濾后不確定的堿基的比例。采用Trinity軟件(http://trinityrnaseq.github.io/)進(jìn)行de novo組裝，首先通過序列之間的overlap將序列延伸為Contig，再根據(jù)序列的雙末端信息(paired-end reads)，將Contig連接，得到該是樣品的Unigene，分析Contig和Unigene的長度及其分布。

1.4 序列的注釋、功能分類和生物學(xué)通路分析

獲得的Unigene進(jìn)行生物信息學(xué)的分析，包括功能注釋及其功能分類。通過BLAST將拼接后的Unigene序列與蛋白數(shù)據(jù)庫進(jìn)行比對(duì)，以e-value<0.000 01為閥值獲取注釋，包括NCBI非冗余核酸數(shù)據(jù)庫(Non-redundant protein database，NR)和Swiss-Prot蛋白質(zhì)序列數(shù)據(jù)庫(Swiss-Prot protein database，Swiss-Prot)。比對(duì)到NR數(shù)據(jù)庫中，分析E-value、相似性及其物種的分布。根據(jù)NR注釋信息，使用Blast2GO得到GO注釋信息[23]，GO即基因本體論數(shù)據(jù)庫(Gene Ontology)，是一個(gè)國際標(biāo)準(zhǔn)化的基因功能分類數(shù)據(jù)庫，包括3個(gè)相對(duì)獨(dú)立的類別，即參與的生物過程(Biological process)、構(gòu)成的細(xì)胞組分(Cellular component)以及實(shí)現(xiàn)的分子功能(Molecular function)，使用WEGO對(duì)所有的Unigene進(jìn)行GO功能分類，從宏觀上認(rèn)識(shí)云南松基因功能分布特征[24]。對(duì)Unigene進(jìn)行蛋白質(zhì)直系同源數(shù)據(jù)庫(Cluster of orthologous groups，COG)比對(duì)，對(duì)其可能的功能分類統(tǒng)計(jì)[25]；同時(shí)也進(jìn)行基因組百科全書(Kyoto encyclopedia of genes and genomes，KEGG)代謝通路分析[26]，進(jìn)一步了解云南松的代謝途徑、生物學(xué)功能以及基因間的相互作用等。此外，按NR、Swiss-Prot、KEGG和COG的優(yōu)先順序?qū)nigene與這些蛋白庫進(jìn)行blastx比對(duì)(e-value<0.000 01)，將各Unigene比對(duì)到蛋白庫的CDS(Coding DNA Sequence)，未比對(duì)上的使用ESTScan預(yù)測出CDS。

2 結(jié)果與分析

2.1 云南松轉(zhuǎn)錄組測序及其de novo組裝

測序共獲得原始序列數(shù)據(jù)97 126 960條，對(duì)原始序列經(jīng)去除測序接頭、重復(fù)冗余、低質(zhì)量等過濾處理，獲得有效序列片段95 003 826條，其中，中間未知堿基序列的片段N為0，Q20高質(zhì)量序列占94.9%，GC含量占總堿基數(shù)的45.67%。說明云南松轉(zhuǎn)錄組Illumina Hiseq 2000平臺(tái)測序獲得序列的數(shù)量及其質(zhì)量均較高，為數(shù)據(jù)組裝提供較好的原始序列數(shù)據(jù)，以滿足后續(xù)的分析所用。

對(duì)測序的reads片段進(jìn)行denovo組裝，共獲得122 115條序列重疊群(Contig)，將這些Contig進(jìn)一步組裝，得到80 000條Unigenes，序列信息達(dá)71 215 401 nt，其中聚類(clusters)的Unigenes為25 226條，單獨(dú)(singletons)的Unigenes為54 774條。

組裝的序列長度也可以反映組裝的質(zhì)量，對(duì)contig的序列長度分析可知(圖1:A)，平均長度462 nt，N50為1 240 nt，其中200～1 000 nt長度的序列占87.89%，1 000～2 000 nt長度的序列占7.15%，2 000～3 000 nt的占3.14%，≥3 000 nt的序列占1.83%。對(duì)Unigene的長度進(jìn)行統(tǒng)計(jì)(圖1:B)，平均長度為890 nt，N50為1 818 nt，有30%的序列長度大于1 000 nt，其中12%的大于2 000 nt。由此可以看出，通過轉(zhuǎn)錄組測序獲得大量的序列，經(jīng)組裝后Unigenes的長度明顯增加，平均長度增加近一倍。對(duì)于1 000 nt以上序列的分布也明顯增加，所占比例由12%增加至30%，表明組裝的效果較好，片段長度明顯增加，可進(jìn)一步開展后續(xù)分析。

2.2 云南松轉(zhuǎn)錄組Unigene的NR功能注釋

將組裝得到的80 000條Unigenes通過blast與NR庫進(jìn)行比對(duì)(evalue<0.000 01)，43 434條Unigenes在NR數(shù)據(jù)庫中找到相似序列，占總Unigenes的54.29%，結(jié)果如圖2所示。

圖1 云南松轉(zhuǎn)錄組組裝序列長度分布 A. Contigs的長度分布；B. Unigenes的長度分布Fig.1 Length distribution of assembly contigs and Unigenes of transcriptome for P.yunnanensis A. Length distribution of assembly contigs; B. Length distribution of assembly unigenes

圖2 云南松轉(zhuǎn)錄組Unigene的NR注釋分類 A. NR注釋的e-value分布；B. NR注釋的相似度分布；C. NR注釋的物種分布Fig.2 Category of NR annotation of transcriptome for unigenes of P.yunnanensis A. E-value distribution of NR annotation; B. Similarity distribution of NR annotation; C. Species distribution of NR annotation

E-value分布可以看出(圖2:A)，在注釋的43 434條Unigenes中，有近一半(49.2%)的E值分布于e-45～e-5，36.9%的分布于e-100～e-45，在e=0的情況下占14%。從匹配序列的相似度分布可以看出(圖2:B)，有29.1%的序列的相似度在60%～80%，有28.8%的相似度>80%，42.1%的相似度<60%。從注釋匹配的物種分布來看(圖2:C)，云南松43 434條Unigenes與其它物種已知基因具有不同程度的同源性，注釋序列分布較多的5個(gè)物種分別為北美紅杉(Piceasitchensis)、葡萄(Vitisvinifera)、桃(Amygdaluspersica)、火炬松(Pinustaeda)、蓖麻(Ricinuscommunis)，分別占46.2%、12.2%、3.6%、3.1%和3.1%，其余近1/3的分布于其它400多個(gè)物種中。從E-value和相似度分布來看，云南松在NCBI的NR庫中比對(duì)的匹配度較高，但由于缺乏云南松的轉(zhuǎn)錄組和基因組信息，部分Unigene在數(shù)據(jù)庫中無法匹配到已知的基因中。

2.3 云南松轉(zhuǎn)錄組Unigene的COG注釋及其分類

將獲得的Unigenes與COG數(shù)據(jù)庫進(jìn)行比對(duì)，預(yù)測Unigenes可能的功能，根據(jù)比對(duì)結(jié)果對(duì)Unigenes的功能分類并統(tǒng)計(jì)(圖3)。

結(jié)果表明，有14 792條Unigenes注釋到COG庫，占總Unigenes的18.49%，共獲得29 002個(gè)COG功能注釋信息，分布于25個(gè)功能分類，不同類別的基因表達(dá)豐度各不相同，以一般功能基因(General function prediction only)最多，占17.21%；其次是轉(zhuǎn)錄功能(Transcription)，占9.04%；復(fù)制、重組和修飾(Replication,recombination and repair，8.59%)，信號(hào)傳遞機(jī)制(Signal transduction mechanisms，7.51%)，翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(Posttranslational modification,protein turnover,chaperones，6.52%)以及碳水化合物轉(zhuǎn)運(yùn)和代謝(Carbohydrate transport and metabolism，5.94%)的功能基因也較豐富。而以核結(jié)構(gòu)(Nuclear structure)和細(xì)胞外結(jié)構(gòu)(Extracellular structures)的最少，分別僅有3個(gè)(0.01%)和26個(gè)(0.09%)。這些結(jié)果表明，云南松在轉(zhuǎn)錄、翻譯和信號(hào)傳遞等基因表達(dá)豐度較高。此外，有2 089個(gè)功能注釋信息未知，未確定其準(zhǔn)確的生物學(xué)功能，占所有功能注釋信息的7.20%。

圖3 云南松轉(zhuǎn)錄組Unigene的COG功能注釋分布Fig.3 COG functional annotation distribution of unigenes of transcriptome for P.yunnanensis A.RNA processing and modification; B.Chromatin structure and dynamics; C.Energy production and conversion; D.Cell cycle control,cell division,chromosome partitioning; E.Amino acid transport and metabolism; F.Nucleotide transport and metabolism; G.Carbohydrate transoport and metabolism; H.Coenzyme transport and metabolism; I.Lipid transport and metabolism; J.Translation,ribosomal structure and biogenesis; K.Transcription; L.Replication,recombination and repair; M.Cell wall/membrane/envelope biogenesis; N.Cell motility; O.Posttranslational modification,protein turnover,chaperones; P.Inorganic ion transport and metabolism; Q.Secondary metabolites biosynthesis,transport and catabolism; R.General function prediction only; S.Function unknown; T.Signal transduction mechanisms; U.Intracellular trafficking,secretion,and vesicular transport; V.Defense mechanisms; W.Extracellular structures; Y.Nuclear structure; Z.Cytoskeleton

2.4 云南松轉(zhuǎn)錄組Unigene的GO注釋及其分類

根據(jù)NR注釋進(jìn)一步進(jìn)行GO功能注釋，按照這些基因參與的生物過程、構(gòu)成的細(xì)胞組分以及實(shí)現(xiàn)的分子功能進(jìn)行分類統(tǒng)計(jì)，從宏觀上了解云南松基因功能的分布特征，便于理解基因所代表的生物學(xué)意義(圖4)。

分析可知(圖4)，共有26 743條Unigenes(占總Unigenes的33.43%)注釋194 586個(gè)GO功能，對(duì)這些功能注釋進(jìn)行分類，生物過程、細(xì)胞組分和分子功能3大類分別注釋94 703、70 538和29 345個(gè)基因，其中以生物過程注釋的居多，占48.67%，其次是細(xì)胞組分，占36.25%，以分子功能較少(15.08%)。3個(gè)功能分類又可細(xì)分為55個(gè)功能亞類，分別包括22、17和16個(gè)亞類。在生物過程包含的22功能亞類中，以細(xì)胞過程(cellular process)、代謝過程(metabolic process)和單一有機(jī)體過程(single-organism process)較多，分別占該類型的15.93%、15.47%和11.25%，以運(yùn)動(dòng)(locomotion)最低，僅占生物過程全部的0.03%。細(xì)胞組分中以細(xì)胞(cell)和細(xì)胞組分(cell part)較多，均占該類別的25.52%，其次是細(xì)胞器(organelle)，占該類別的20.00%，而以病毒體(virion)和病毒體組分(virion part)較低，兩者均占該類別的0.003%。分子功能中以結(jié)合(binding)和催化活性(catalytic activity)居多，分別占該類型的47.74%和37.89%，而以通道調(diào)節(jié)活性(channel regulator activity)、翻譯調(diào)節(jié)活性(translation regulator activity)、金屬伴侶蛋白活性(metallochaperone activity)、蛋白標(biāo)簽(protein tag)較低，分別占該類別的0.003%、0.007%、0.014%和0.017%。這些結(jié)果顯示了云南松新梢基因表達(dá)的總體情況，不難看出，以代謝活動(dòng)相關(guān)的基因量較多，表明云南松代謝能力較強(qiáng)。

2.5云南松轉(zhuǎn)錄組Unigene的KEGG代謝通路分析

將Unigene比對(duì)到KEGG數(shù)據(jù)庫，根據(jù)注釋信息分析代謝通路，了解基因產(chǎn)物在細(xì)胞中的代謝途徑及其基因產(chǎn)物的功能，結(jié)果表明共有25 873個(gè)Unigenes(32.34%)獲得注釋，對(duì)其可能參與或涉及的代謝通路統(tǒng)計(jì)分析，可將云南松Unigenes歸為5個(gè)類別(Level 1)、20個(gè)亞類(Level 2)、128個(gè)代謝通路(Pathway)(圖5)。

圖4 云南松轉(zhuǎn)錄組Unigene的GO功能分類Fig.4 GO classification of unigenes of transcriptome for P.yunnanensis

圖5 云南松轉(zhuǎn)錄組Unigene的KEGG分類Fig.5 KEGG classification of unigenes of transcriptome for P.yunnanensis

由圖5可以看出，5大類別中(Level 1)，以代謝(Metabolism)相關(guān)的通路所占比例最多，為64.38%，其次為遺傳信息處理(Genetic Information Processing)相關(guān)的通路，占18.43%，而以細(xì)胞過程(Cellular Processes)相關(guān)的通路最少，僅占3.43%，環(huán)境信息處理(Environmental Information Processing)和生物系統(tǒng)(Organismal Systems)相關(guān)的通路分別占5.85%和7.92%。進(jìn)一步細(xì)分為亞類(Level 2)，其中代謝相關(guān)的通路可細(xì)分為11亞類，包括氨基酸代謝(Amino acid metabolism)、其他次生物質(zhì)代謝(Biosynthesis of other secondary metabolites)、碳水化合物代謝(Carbohydrate metabolism)、能量代謝(Energy metabolism)、全局整體映射(Global map)、糖生物合成和代謝(Glycan biosynthesis and metabolism)、脂類物質(zhì)代謝(Lipid metabolism)、輔助因子和維生素代謝(Metabolism of cofactors and vitamins)、其它氨基酸代謝(Metabolism of other amino acids)、萜類化合物和聚酮化合物代謝(Metabolism of terpenoids and polyketides)和核苷酸代謝(Nucleotide metabolism)。而細(xì)胞過程僅包括運(yùn)輸和代謝(Transport and catabolism)1個(gè)亞類，其余遺傳信息處理、環(huán)境信息處理和生物系統(tǒng)相關(guān)的通路分別包括4、2、2個(gè)亞類，共20個(gè)亞類。這些亞類中，所注釋的Unigenes從高到低的依次為：全局整體映射(25.12%)、碳水化合物代謝(8.63%)、環(huán)境適應(yīng)(Environmental adaptation，7.71%)、翻譯(Translation，6.79%)、脂類物質(zhì)代謝(6.33%)、其他次生物質(zhì)代謝(5.77%)、信號(hào)傳導(dǎo)(Signal transduction，4.96%)、折疊、分類和降解(Folding,sorting and degradation，4.78%)、氨基酸代謝(4.68%)、轉(zhuǎn)錄(Transcription，4.43%)、核苷酸代謝(Nucleotide metabolism，4.18%)、運(yùn)輸和代謝(3.43%)等，其中，以代謝相關(guān)的通路居多，表明云南松在這一時(shí)期具有較強(qiáng)的代謝活動(dòng)。

5個(gè)Level 1、20個(gè)Level 2的Unigene定位到128個(gè)具體的KEGG代謝通路。按基因獲得的注釋量高低排序，將前20個(gè)代謝通路列于表1。從表1可以看出，以代謝途徑(Metabolic pathway)的最多，占全部的16.42%，其次為次生代謝產(chǎn)物合成和植物病原互作(Plant-pathogen interaction)，分別占8.70%和6.99%。

此外，為探討與木質(zhì)素合成相關(guān)Unigenes，本研究涉及到的苯丙類生物合成(Phenylpropanoid biosynthesis)的Unigenes有619條。進(jìn)一步對(duì)該代謝通路進(jìn)行分析，找到與木質(zhì)素生物合成相關(guān)的酶(表2)，包括苯丙氨酸氨裂解酶(phenylalanine ammonia-lyase，PAL)、4-香豆酰輔酶A連接酶(4-coumarate-CoA ligase，4CL)、莽草酸/奎寧酸羥基肉桂酰轉(zhuǎn)移酶(shikimate O-hydroxycinnamoyltransferase，HCT)、咖啡酸輔酶A-O-甲基轉(zhuǎn)移酶(caffeoyl-CoA O-methyltransferase，CCoAOMT)、肉桂酰輔酶A還原酶(cinnamoyl-CoA reductase，CCR)、阿魏酸-5-羥基化酶(ferulate-5-hydroxylase，F(xiàn)5H)、咖啡酸3-O-甲基轉(zhuǎn)移酶(caffeic acid 3-O-methyltransferase，COMT)等。

表1云南松轉(zhuǎn)錄組Unigene的KEEGpathway注釋

Table1TheKEEGpathwayannotationofunigenesoftranscriptomeforP.yunnanensis

代謝通路Pathway代謝通路編號(hào)PathwayID注釋Unigene數(shù)Numnerofunigene比例Percentage(%)代謝途徑Metabolicpathwayko01100515716．42次生代謝產(chǎn)物合成Biosynthesisofsecondarymetabolitesko0111027348．70植物病原互作Plant?pathogeninteractionko0462621966．99植物激素信號(hào)傳導(dǎo)Planthormonesignaltransductionko0407513684．35剪接體Spliceosomeko030408512．71RNA轉(zhuǎn)運(yùn)RNAtransportko030138322．65嘧啶代謝Pyrimidinemetabolismko002406572．09嘌呤代謝Purinemetabolismko002306572．09淀粉和蔗糖代謝Starchandsucrosemetabolismko005006412．04苯丙生物合成Phenylpropanoidbiosynthesisko009406191．97內(nèi)吞作用Endocytosisko041445271．68甘油磷脂代謝Glycerophospholipidmetabolismko005644901．56內(nèi)質(zhì)網(wǎng)蛋白加工Proteinprocessinginendoplasmicreticulumko041414741．51黃酮類化合物生物合成Flavonoidbiosynthesisko009414721．50RNA聚合酶RNApolymeraseko030204521．44mRNA監(jiān)控途徑mRNAsurveillancepathwayko030154351．38真核生物核糖體生物合成Ribosomebiogenesisineukaryotesko030084001．27RNA降解RNAdegradationko030183871．23泛素介導(dǎo)的蛋白降解途徑Ubiquitinmediatedproteolysisko041203801．21核糖體Ribosomeko030103531．12

表2云南松轉(zhuǎn)錄組中參與木質(zhì)素合成的酶

Table2TheenzymerelatedtoligninbiosynthesisoftranscriptomeforP.yunnanensis

EC號(hào)酶名稱Enzymedefinition酶的縮寫Enzymeabbr．編號(hào)EntryEC：4．3．1．24苯丙氨酸氨裂解酶Phenylalanineammonia?lyasePALK10775EC：6．2．1．124?香豆酰輔酶A連接酶4?coumarate?CoAligase4CLK01904EC：1．2．1．44肉桂酰輔酶A還原酶Cinnamoyl?CoAreductaseCCRK09753EC：2．1．1．68咖啡酸3?O?甲基轉(zhuǎn)移酶Caffeicacid3?O?methyltransferaseCOMTK13066EC：2．1．1．104咖啡酸輔酶A?O?甲基轉(zhuǎn)移酶Caffeoyl?CoAO?methyltransferaseCCoAOMTK00588EC：1．14．?．?阿魏酸?5?羥基化酶Ferulate?5?hydroxylaseF5HK09755EC：2．3．1．133莽草酸/奎寧酸羥基肉桂酰轉(zhuǎn)移酶ShikimateO?hydroxycinnamoyltransferaseHCTK13065

圖6 云南松轉(zhuǎn)錄組Unigene的CDS序列長度分布Fig.6 CDS length distribution of transcriptome for P.yunnanensis

2.6 云南松轉(zhuǎn)錄組Unigene的CDS預(yù)測

按NR、Swiss-Prot、KEGG和COG的優(yōu)先級(jí)順序?qū)nigene序列與以上的蛋白庫進(jìn)行比對(duì)(evalue<0.000 01)，不能比對(duì)的用ESTScan預(yù)測其編碼區(qū)，結(jié)果有42 668條Unigenes能比對(duì)到蛋白庫，另預(yù)測3 697條的CDS，其長度分布如圖6所示。從圖6A可以看出，比對(duì)的CDS序列長度1 000 nt以上的占26%，其中1 000～2 000 nt的占總數(shù)的19.33%，2 000～3 000 nt的占4.56%，3 000 nt以上的占2.15%。預(yù)測的CDS集中分布于200～1 000 nt，占96.29%，有極少數(shù)(0.27%)的序列在3 000 nt以上(圖6B)。

3 討論

研究采用Illumina高通量測序的數(shù)據(jù)量大，信息較多，是獲得大量EST的有效途徑，對(duì)于至今未開展基因組測序的物種而言，高通量測序是挖掘該物種生長發(fā)育過程中表達(dá)的重要基因可取途徑，已應(yīng)用于多種物種的基因研究，如輻射松(P.radiata)[14]、油松(P.tabuliformis)[15]、地中海白松(P.halepensisMill.)[16]、高山松(P.densata)[17]等。本研究通過云南松轉(zhuǎn)錄組測序，經(jīng)數(shù)據(jù)組裝后獲得80 000條Unigenes，與龐大的松樹全基因組相比[27]，本研究獲得的信息量大，可為后續(xù)的基因功能分析、基因克隆、分子標(biāo)記的開發(fā)等方面奠定基礎(chǔ)。通過序列分析，組裝后獲得的Unigene序列長度平均為890 nt，與其它松樹如馬尾松[28～29]、油松[15]、地中海白松[16]等相比，本研究云南松轉(zhuǎn)錄組序列拼接平均長度較好，與此同時(shí)，組裝后的N50達(dá)到1 818 nt，序列較長。總體來看，從序列的數(shù)量、序列的Q20、序列的長度分布、N值等方面比較，獲得的序列數(shù)量和質(zhì)量都比較高，有利于后續(xù)的研究分析。測序拼接后的序列clusters有25 226條、singletons有54 774條，singleton的序列比較多，可能是拼接時(shí)嚴(yán)格的參數(shù)選擇或測序時(shí)低表達(dá)基因的出現(xiàn)等，其中大量的singleton也獲得注釋，說明這些singleton也代表轉(zhuǎn)錄組信息，為高質(zhì)量有用的reads，這在柳樹(Salixspp.)的研究中也提到[30]。

采用生物信息學(xué)分析方法，對(duì)云南松轉(zhuǎn)錄組的序列進(jìn)行功能注釋及其功能分類。通過NR數(shù)據(jù)庫的比對(duì)，有43 434條序列比對(duì)上，占總Unigene的54.29%。通過COG數(shù)據(jù)庫比對(duì)，獲得了大量多方面的表達(dá)基因，這些基因反映了云南松在一定時(shí)期的表達(dá)情況。按其功能可把14 792條序列包含的基因分為25類，其中以一般功能、轉(zhuǎn)錄較為豐富，而有7.2%(2 089條)的Unigenes功能未知，難以確定其生物功能，可能是注釋信息不完善造成的，在很多林木轉(zhuǎn)錄組分析中也出現(xiàn)同樣的情況，如泡桐(Paulownia)[31]、遼東櫟(Ouercusliaotungensis)[32]、錐栗(Castaneahenryi)[33]、柳樹[30]等。轉(zhuǎn)錄組中的Unigenes根據(jù)GO功能分為生物過程、細(xì)胞組分和分子功能3大類55個(gè)亞類，以分布于細(xì)胞過程、代謝過程、細(xì)胞、細(xì)胞組分、細(xì)胞器、結(jié)合和催化活性等方面較為集中；根據(jù)KEGG代謝通路分析，共有25 873個(gè)Unigenes(32.34%)注釋到對(duì)應(yīng)的128個(gè)代謝通路中，其中以代謝相關(guān)通路的占多數(shù)，表明云南松具有比較豐富的代謝活動(dòng)。通過注釋的信息量來看，云南松基因含量豐富，為其適應(yīng)性強(qiáng)提供了分子生物學(xué)方面的解釋與支持。盡管這些Unigenes并不覆蓋云南松整個(gè)蛋白編碼區(qū)，但這些數(shù)據(jù)也提供較為豐富的表達(dá)信息，為今后云南松功能基因的挖掘、利用、分子標(biāo)記的開發(fā)奠定基礎(chǔ)?？偟奈醋⑨孶nigenes有27 825條，占總Unigene的34.78%，這些序列可能是因?yàn)楸旧頌榉蔷幋a的RNA序列，或者是因長度不足未包含蛋白質(zhì)功能域的信息，也有可能是數(shù)據(jù)庫中的基因信息不足或云南松特有的新基因未能匹配上[29]。

4 結(jié)論

通過轉(zhuǎn)錄組測序，獲得的Unigenes與公共數(shù)據(jù)庫進(jìn)行比對(duì)，有65.22%的Unigenes獲得注釋信息，其中COG分析將其分為25個(gè)不同功能分類，以一般功能基因和轉(zhuǎn)錄功能的較多；GO分析得到55個(gè)不同功能分類，以生物學(xué)過程占多數(shù)，各亞類中以細(xì)胞過程、代謝過程等類別的基因含量豐富，表明云南松具有強(qiáng)的代謝活力及其代謝機(jī)制；KEGG的128個(gè)代謝通路，以涉及代謝途徑和次生代謝產(chǎn)物合成較為豐富，并找到涉及到木質(zhì)素合成的關(guān)鍵酶多個(gè)，有助于將來開展云南松材質(zhì)的研究。研究結(jié)果揭示了云南松基因豐富，并分析基因的分布特征，確定基因代謝通路的分布，為云南松較強(qiáng)的適應(yīng)性提供了轉(zhuǎn)錄組水平上的支持，也為后續(xù)云南松基因組學(xué)、功能基因、分子標(biāo)記的開發(fā)等方面的研究奠定基礎(chǔ)。

1.金振洲,彭鑒.云南松[M].昆明:云南科技出版社,2004:1-66.

2.中國科學(xué)院昆明植物研究所編著.云南植物志(第4卷):種子植物[M].北京:科學(xué)出版社,1986:54-57.

3.中國科學(xué)院中國植物志編輯委員會(huì).中國植物志:第7卷[M].北京:科學(xué)出版社,1978:122-282.

4.陳飛,王健敏,孫寶剛,等.云南松的地理分布與氣候關(guān)系[J].林業(yè)科學(xué)研究,2012,25(2):163-168.

5.陳飛,王健敏,陳曉鳴,等.基于Kira指標(biāo)的云南松氣候適宜性分析[J].林業(yè)科學(xué)研究,2012,25(5):576-581.

6.鄧喜慶,皇寶林,溫慶忠,等.云南松林資源動(dòng)態(tài)研究[J].自然資源學(xué)報(bào),2014,29(8):1411-1419.

7.鄧喜慶,皇寶林,溫慶忠,等.云南松林在云南的分布研究[J].云南大學(xué)學(xué)報(bào):自然科學(xué)版,2013,35(6):843-848.

8.Zhang L,Xu W H,Ouyang Z Y,et al.Determination of priority nature conservation areas and human disturbances in the Yangtze River Basin,China[J].Journal for Nature Conservation,2014,22(4):326-336.

9.Xiao X Y,Haberle S G,Shen J,et al.Latest Pleistocene and Holocene vegetation and climate history inferred from an alpine lacustrine record,northwestern Yunnan Province,southwestern China[J].Quaternary Science Reviews,2014,86(2):35-48.

10.Wu J X,Wang Y X,Chen Q B,et al.Soil improvement of Pinus yunnanensis forest at different age in Central Yunnan Plateau[J].Advanced Materials Research,2014,864-867:2565-2568.

11.Lin Y M,Cui P,Ge Y G,et al.The succession characteristics of soil erosion during different vegetation succession stages in dry-hot river valley of Jinsha River,upper reaches of Yangtze River[J].Ecological Engineering,2014,62(1):13-26.

12.Wang B S,Mao J F,Zhao W,et al.Impact of Geography and Climate on the Genetic Differentiation of the Subtropical PinePinusyunnanensis[J].PLoS ONE,2013,8(6):e67345.

13.Xu Y L,Zhang R L,Tian B,et al.Development of novel microsatellite markers forPinusyunnanensisand their cross amplification in congeneric species[J].Conservation Genetics Resources,2013,5(4):1113-1114.

14.Li X,Yang X,Wu H X.Transcriptome profiling of radiata pine branches reveals new insights into reaction wood formation with implications in plant gravitropism[J].BMC Genomics,2013,14(22):547-554.

15.Niu S H,Li Z X,Yuan H W,et al.Transcriptome characterisation ofPinustabuliformisand evolution of genes in thePinusphylogeny[J].BMC Genomics,2013,14(1):202-207.

16.Pinosio S,González-Martínez S C,Bagnoli F,et al.First insights into the transcriptome and development of new genomic tools of a widespread circum-mediterranean tree species,PinushalepensisMill.[J].Molecular Ecology Resources,2014,14(4):846-856.

17.Wan L C,Zhang H,Lu S,et al.Transcriptome-wide identification and characterization of mirnas fromPinusdensata[J].BMC Genomics,2012,13(1):132-142.

18.Fang P,Niu S,Yuan H,et al.Development and characterization of 25 EST-SSR markers inPinussylvestrisvar.mongolica(Pinaceae)[J].Applications in Plant Sciences,2014,2(1):1300057.

19.Feng Y H,Yang Z Q,Wang J,et al.Development and characterization of SSR markers fromPinusmassonianaand their transferability toP.elliottii,P.caribaeaandP.yunnanensis[J].Genetics & Molecular Research Gmr,2014,13(1):1508-1513.

20.Iwaizumi M G,Tsuda Y,Ohtani M,et al.Recent distribution changes affect geographic clines in genetic diversity and structure ofPinusdensifloranatural populations in Japan[J].Forest Ecology & Management,2013,304(4):407-416.

21.Lesser M R,Parchman T L,Buerkle C.Cross-species transferability of SSR loci developed from transciptome sequencing in lodgepole pine[J].Molecular Ecology Resources,2012,12(3):448-455.

22.Liu J J,Hammett C.Development of novel polymorphic microsatellite markers by technology of next generation sequencing in western white pine[J].Conservation Genetics Resources,2014,6(3):647-648.

23.Conesa ,G?tz S,García-Gómez J M,et al.Blast2GO:A universal tool for annotation,visualization and analysis in functional genomics research[J].Bioinformatics,2005,21(18):3674-3676.

24.Ye J,Fang L,Zheng H,et al.WEGO:a web tool for plotting GO annotations[J].Nucleic Acids Research,2006,34(2):293-297.

25.Tatusov R L,Galperln M Y,Natale D A,et al.The COG database:A tool for genome-scale analysis of protein functions and evolution[J].Nucleic Acids Research,2000,28(1):33-36.

26.Kanehisa M,Goto S,Kawashima S,et al.The KEGG resource for deciphering the genome[J].Nucleic Acids Research,2004,32(22):277-280.

27.Zimin A,Stevens K A,Crepeau M W,et al.Sequencing and assembly of the 22-gb loblolly pine genome[J].Genetics,2014,196(3):875-890.

28.王曉峰,何懷衛(wèi)龍,蔡衛(wèi)佳,等.馬尾松轉(zhuǎn)錄組測序和分析[J].分子植物育種,2013,11(3):385-392.

29.Bai T D,Xu L,Xu M,et al.Characterization of masson pine(PinusmassonianaLamb.) microsatellite DNA by 454 genome shotgun sequencing[J].Tree Genetics & Genomes,2014,10(2):429-437.

30.鄭紀(jì)偉.柳樹轉(zhuǎn)錄組高通量測序及SSR標(biāo)記開發(fā)研究[D].南京:南京林業(yè)大學(xué),2013.

31.鄧敏捷,董焱鵬,趙振利,等.基于Illumina高通量測序的泡桐轉(zhuǎn)錄組研究[J].林業(yè)科學(xué),2013,49(6):30-36.

32.劉玉林,李偉,張志翔.基于高通量測序的遼東櫟轉(zhuǎn)錄組學(xué)研究[J].生物技術(shù)通報(bào),2014(7):119-124.

33.張琳,范曉明,林青,等.錐栗種仁轉(zhuǎn)錄組及淀粉和蔗糖代謝相關(guān)酶基因的表達(dá)分析[J].植物遺傳資源學(xué)報(bào),2015,16(3):603-612.

TranscriptomeAnalysisforPinusyunnanensisBasedonHighThroughputSequencing

CAI Nian-Hui1,2,3DENG Li-Li1,2,3XU Yu-Lan1,2,3XU Yang1,2,3ZHOU Li1,2,3WANG Da-Wei1,2,3TIAN Bin1,2,3HE Cheng-Zhong1,2,3DUAN An-An1,2,3

(1.Southwest Forestry University,Key Laboratory for Forest Genetic and Tree Improvement & Propagation in Universities of Yunnan Province,Kunming 650224；2.Key Laboratory of Biodiversity Conservation in Southwest China,State Forestry Administration,Southwest Forestry University,Kunming 650204；3.Key laboratory for Forest Resources Conservation and Use in the Southwest Mountains of China,Ministry of Education,Southwest Forestry University,Kunming 650224)

The transcriptome ofPinusyunnanensiswas sequenced by using Illumina Hiseq 2 000. In total 80 000 Unigene with an average length of 890 nt and N50 of 1 881 nt were obtained by de novo assembly. Of the Unigene, 43 434, 46 415 and 29 418 Unigenes had significant similarity with known data bank in NR, NT and Swiss-Prot, respectively. 14 792 Unigenes were annotated in clusters of orthologous groups of proteins(COG) and assigned to 25 clusters. 26 743 Unigenes were annotated in gene ontology(GO) and grouped into biological processes, cellular components and molecular function three functional categories, 55 sub-categories. The biological processes were most commonly existed. A total of 25 873 Unigenes were divided into 128 Kyoto Encyclopedia of Genes and Genomes(KEGG) pathways whose functions focused on metabolism. We found some Unigenes related to lignin biosynthesis. The sequence data forP.yunnanensiswiill be helpful for the gene discovery and utilization, molecular marker development and genetic improvement in the further research.

Pinusyunnanensis;transcriptome;Illumina high throughput sequencing technology

國家自然科學(xué)基金項(xiàng)目(31360189)和(31260191)

蔡年輝(1975—)，男，碩士，講師，主要從事森林培育研究。

2015-06-12

S791.257

A

10.7525/j.issn.1673-5102.2016.01.011