滕 云 邱宗波 王穩(wěn)利 張雪如

（1. 信陽農林學院園藝學院，信陽 464000；2. 河南師范大學生命科學學院，新鄉(xiāng) 453007）

林木的次生生長是多年生木本植物維管形成層周期性分裂，向外分化出韌皮部，向內分化出木質部的過程［1］。維管形成層細胞的增殖和分化使樹木不斷增粗，但在不同類型以及不同地區(qū)的樹木中有著不同的活動樣式。對于溫帶的木本植物來說，冬季形成層由于低溫進入休眠，到次年春天氣溫回升時開始萌動，夏季形成層由于高溫進入旺盛的活躍期［2］。維管形成層隨季節(jié)由休眠期到萌動期再到活躍期，再由活躍期到休眠期的循環(huán)變化在樹木的生長和發(fā)育過程中起著關鍵的作用［2～3］。

隨著高通量測序技術的發(fā)展，轉錄組（RNASeq）測序已成為了解植物在某一發(fā)育階段或功能狀態(tài)下基因功能及基因結構的重要手段［3～4］。曹玉婷等［5］利用高通量RNA-Seq 測序技術對4 個不同發(fā)育階段的柳杉維管形成層進行轉錄組測序，獲得了大量與木材形成相關的重要功能基因。楊光等［6］對柚木邊材進行轉錄組測序，初步明確了部分編碼蛋白的功能、分類及代謝通路，為柚木分子育種工作的開展提供數(shù)據(jù)和參考。然而泡桐維管形成層及木質部區(qū)組織轉錄組測序的相關研究尚未見報道。

泡桐（Paulownia spp.）為玄參科（Scrophulariaceae）泡桐屬（Paulownia）樹種，具有生長快和用途廣等優(yōu)點，是我國重要的速生優(yōu)質用材樹種之一［7］。泡桐中已發(fā)現(xiàn)的基因大多是與抗病害相關的基因，與木材材質和產量形成相關基因的挖掘還相對比較薄弱［8］。因此，本研究利用HiSeqTM2500 測序平臺對泡桐休眠期、萌動期和活躍期形成層及木質部區(qū)組織的混合樣品進行了轉錄組測序和組裝，并利用生物信息學的方法，對相關基因進行功能注釋和分析，從而為泡桐功能基因的發(fā)掘利用，揭示木材形成的分子機制提供理論依據(jù)。

1 材料與方法

1.1 植物材料

實驗材料取自河南農業(yè)大學樹木園內生長良好樹齡大約為20 年的泡桐。取材時間分別為2018年1月27日（休眠期）、2018年3月18日（萌動期）、2018 年6 月4 日（活躍期），每個時間點選取3株生長一致的泡桐作為試材，分別在1.3 m 樹高處用鋒利的鑿子去除樹皮，然后用小刀刮取形成層，迅速裝入無RNA 酶的Eppendorf離心管中，用液氮速凍后儲存在-80℃冰箱中。

1.2 RNA的提取、文庫的制備及測序

采用Plant RNA Reagent 試劑（Invitrogen），按照說明書分別提取泡桐休眠期、萌動期和活躍期維管形成層及木質部區(qū)的總RNA?？俁NA 樣品的質量檢測是利用Agilent 2100 生物分析儀（Agilent Technologies，USA）進行的。休眠期、萌動期和活躍期泡桐形成層及木質部區(qū)組織的總RNA分別檢測合格后，將3 個時期的RNA 等量混合用于測序文庫構建。cDNA 文庫的制備參照賈新平等［9］方法。建好的文庫利用北京諾禾致源生物公司Illumina HiSeqTM2500測序平臺進行轉錄組測序。

1.3 數(shù)據(jù)組裝及功能注釋

將測序得到的原始序列（raw reads）去除接頭、poly-N 序列、小于20 nt 的序列以及低質量測序片段后，得到高質量干凈數(shù)據(jù)（clean reads）。利用Trinity 軟件（參數(shù)設置K-mer=25）［10］對clean reads進行De novo 拼接組裝，得到長度不一的轉錄本（Transcripts），并從中選取最長的轉錄本作為泡桐的單基因簇（Unigene）。通過Blastx 軟件［11］將泡桐的Unigenes 序列分別與Nr（Non-redundant protein database）、GO（Gene ontology）、Nt（Non-redundant nucleotide sequences）、SwissProt（SwissProt protein database）、COG（Clusters of orthologous group）、Pfam（Protein families database）和KEGG（Kyoto encyclopedia of genes and genomes）數(shù)據(jù)庫進行比對，從而得到泡桐Unigenes 的功能注釋信息。

1.4 泡桐SSR位點分析

使用MISA 軟件（http：//pgrc.ipk-gatersleben.de/misa/misa.html）從泡桐Unigenes 中搜索簡單序列重復（SSR）位點，搜索標準參照許薔薇等［12］方法。

2 結果與分析

2.1 泡桐RNA-Seq測序及轉錄組組裝

為了獲得盡可能豐富的泡桐轉錄組信息，選取泡桐維管形成層的休眠期、萌動期和活躍期材料分別提取總RNA，然后將這3 個組織的RNA 等量混合后進行轉錄組測序。經(jīng)IllluminaTMHiSeq 2500 高通量平臺測序，總共獲得111 645 144 條原始序列（raw reads），將原始數(shù)據(jù)進行過濾后，得到109 021 918 條clean reads（約16.35 Gb）。原始數(shù)據(jù)已上傳到NCBI 中的GSE 數(shù)據(jù)庫（登錄號為GSE109323）。另外Q20為97.77%，Q30為94.30%，表明測序質量較好，可為后續(xù)的序列組裝提供高質量的原始數(shù)據(jù)。

測序得到的clean reads 利用Trinity 軟件進行De novo組裝，得到145 881個轉錄本（Transcripts），平均長度為888 nt（見表1）。其中，長度在1 000 nt以下的Transcripts 有104 177 條，占71.41%；長度在1 000～2 000 nt 的有24 155 條，占16.56%；而長度在2000 nt 以上的僅有17 549 條，占12.03%。進一步對轉錄本聚類去冗長得到104 432 個Unigenes，平均長度為662 nt（見表1）。其中，長度在1 000 nt 以下的Unigenes 有86 607 條，占82.93%；長度在1 000～2 000 nt 的有10 848 條，占10.39%；而長度在2 000 nt 以上的僅有6 977 條，占6.68%。Transcripts 和Unigenes 的N50 分 別 達 到1 655 和1 105 nt，表明組裝片段的完整性是比較高的。

表1 泡桐轉錄組組裝結果Table 1 Summary of transcriptome assembly for P.tomentosa

2.2 Unigene基因功能注釋及分類

為了獲得Unigene 的功能注釋信息，我們將組裝得到的104 432 個Unigenes 分別在Nr、Nt、SwissProt、GO、COG、KEGG和Pfam 數(shù)據(jù)庫中利用blastx軟件進行搜索比對。從表2 可以看出，Nr 數(shù)據(jù)庫注釋的Unigenes最多（40 789條），占總Unigenes的39.05%；僅有15 539 條Unigenes 在COG 數(shù)據(jù)庫得到注釋，占總Unigenes 的14.87%。Nt、SwissProt、GO、KEGG 和Pfam 分 別 注 釋 到31 675、30 499、29 168、16 316、28 828條Unigenes。另外，有48 214條（46.17%）泡桐Unigene 序列可以獲得至少1 個數(shù)據(jù)庫的功能注釋，而且在7個數(shù)據(jù)庫中都能得到注釋的Unigenes 共有7 386 個，占總Unigenes 數(shù)的7.07%。

2.3 Unigene的GO功能分析

泡桐Unigene 通過與GO 數(shù)據(jù)庫比對分析，可將注釋的29 168個Unigene分為分子功能（molecu-lar function）、生物過程（biological process）和細胞組分（cellular component）3 大類55 個功能組（見圖1）。其中，歸于生物過程的Unigene 有69 943 個，38 276個Unigene歸于細胞組分，33 316個Unigene歸于分子功能，這一Unigene 功能注釋分類結果顯示了泡桐維管形成層發(fā)育過程中基因表達的總體情況。在生物過程這一功能分類中，細胞過程（cellular process）（73 324 個）、代謝過程（metabolic process）（162 580 個）和單一有機體過程（single-organism process）（28 373個）功能組注釋到的Unigene 數(shù)量最多，分別占69 943 條生物學過程Unigene的22.74%、22.02%和16.68%，而生物相（biological phase）注釋到的Unigene 最少，僅占0.03%。在細胞組分這一功能分類中，細胞（cell）（19 186 個）和細胞組分（cell part）（30 231 個）功能組注釋到的Unigene 數(shù)量最多，分別占38 276 個細胞組分Unigene 的20.05%和20.04%，突觸部分（synapse part）、突觸（synapse）、胞外基質成分（extracellular matrix component）、共質體（symplast）4 個功能組注釋到的Unigene 均占38 276 個細胞組分Unigene的0.01%。在分子功能這一功能分類中，結合（binding）（49 522 個）和催化活性（catalytic activity）（12 555個）功能組注釋到的Unigene 數(shù)量最多，分別占33 316 個分子功能Unigene 的45.80%和37.85%。

表2 泡桐Unigenes序列的功能注釋Table 2 Functional annotation of P.tomentosa Unigenes

2.4 Unigene的COG功能分析

將104 432條泡桐Unigene與COG數(shù)據(jù)庫進行比對，泡桐維管形成層中有15 539 個Unigenes 得到注釋。從圖2 可以看出，注釋到COG 數(shù)據(jù)庫中的Unigenes 共分為25 類。其中，注釋數(shù)目最多的一類是2 351 個一般功能預測的Unigene，占總數(shù)的15.13%。接下來的一類是2 247個翻譯后修飾、蛋白質折疊、分子伴侶注釋的Unigene 和1 804 個翻譯、核糖體結構和生物發(fā)生注釋的Unigene，分別占總數(shù)的14.16%和11.61%。而核結構注釋的Unigene（52 個）、細胞外結構注釋的Unigene（43個）和細胞活性注釋的Unigene（8 個）數(shù)目較少，占總數(shù)的比例分別為0.33%、0.28%和0.05%。另外590個未知功能注釋的Unigene占總數(shù)的3.80%。

2.5 Unigene的KEGG通路注釋

將104 432 條泡桐Unigene 采用Blastx 比對到KEGG 數(shù)據(jù)庫，共有16 316 條Unigenes 得到注釋，涉及到130條代謝通路中。Unigene注釋序列最多的20個代謝途徑列于表3。由表3可知，最有代表性的通路是“核糖體”（1 216 條）、“植物病原體相互作用”（646條），“碳代謝”（638條），其次是“氨基酸的生物合成”（559條）、“蛋白質在內質網(wǎng)中的加工”（533 條）和“剪接”（396 條）。在苯丙素生物合成（Phenylpropanoid biosynthesis）代謝途徑中，共發(fā)現(xiàn)266 條Unigenes。其中與木質素合成相關的酶有苯丙氨酸氨裂解酶（Phenylalanine ammonia lyase，PAL）、肉桂酸-4-羥化酶（Cinnamic acid-4-hydroxylase，C4H）、肉桂醇脫氫酶（Cinnamyl alcohol dehydrogenase，CAD）、肉桂酰輔酶A 還原酶（Cinnamoyl-CoA reductase，CCR）和咖啡酸輔酶A-O-甲基轉移酶（Caffeic acid coenzyme A-O-methyltransferase，CCoAOMT）。這些Unigenes 功能的注釋為后續(xù)研究木材形成的分子機制起到極其重要的作用。

2.6 泡桐SSR特征分析

通過對泡桐Unigene的搜索，104 432個Unigenes 序列共檢測出16 118 個SSR 位點（見圖3），包含了9 001個單核苷酸、5 024個二核苷酸、2 048個三核苷酸、17 個四核苷酸、14 個五核苷酸和14 個六核苷酸重復序列。其中，單核苷酸中又以A/T重復基元為主（4 427 個）占全部單核苷酸重復的49.18%。AG/CT 類型在二核苷酸重復中最多，三核苷酸重復以AAG/CTT居多。

3 討論

近年來，RNA-Seq 高通量測序技術已成為非模式生物遺傳背景解析的重要手段，已經(jīng)用于木本 植 物 杉 木（Cunninghamia lanceolata）［1］、柳 杉（Cryptomeria fortunei）［5］、柚木（Tectona grandis）［6］、

泡桐（Paulownia tomentosa）［8］和泡泡刺（Nitraria sphaerocarpa）［13］等轉錄組的研究。本研究利用Illumina 高通量測序技術對泡桐維管形成層及木質部區(qū)進行轉錄組測序，獲得109 021 918 條clean reads，堿基Q30為94.30%，當Q30在80%以上就認為測序質量非常可靠［9］。對clean reads 片段進行拼接組裝后共獲得104 432 個Unigenes，平均長度為662 bp。其中N50 值為1 105 bp，說明本研究所得序列中長片段較多，組裝效果較好，有利于進一步的研究分析。

表3 泡桐轉錄組Unigene的KEGG代謝途徑分類Table 3 Classification of unigenes KEGG metabolic pathways in P.tomentosa transcriptome

將泡桐維管形成層及木質部區(qū)組裝得到的Unigenes 與NR、GO、SwissProt、KEGG 和COG 等公共數(shù)據(jù)庫進行序列比對，共有48 214（46.16%）個Unigenes成功獲得注釋，但仍有56 218（53.84%）個序列定位不清楚，未得到功能注釋，這與許多木本植物杉木（Cunninghamia lanceolata）［1］、泡桐（Paulownia tomentosa）［8］、云南松（Pinus yunnanensis）［14］和穗花杉（Amentotaxus argotaenia）［15］均出現(xiàn)大量未獲得注釋Unigene 的結果類似。出現(xiàn)這一現(xiàn)象的原因之一可能是由于這些Unigene 片段較短或是本身為非編碼RNA 序列［16］。研究表明，轉錄組中Unigene 的序列越長，獲得注釋信息的可能性就越高［9］。另外也有可能是這些大量未獲得注釋Unigene 是泡桐特有的新基因，這對后續(xù)開展泡桐維管形成層特異性分析提供重要信息。

RNA-Seq 高通量測序技術不僅使以較低的成本研究非模式生物基因組成為可能，同時也為SSR 位點的挖掘提供了大量序列信息［16］。曹玉婷等［5］在柳杉維管形成層轉錄組序列中鑒定獲得3 772 個SSR 位點，占比最高的為三核苷酸重復SSR。王興春等［17］利用MISA軟件分析了連翹轉錄組，從Unigenes 中獲得3 199 個SSR 位點。其中，占比最高的為單核苷酸重復SSR，其次是二核苷酸重復和三核苷酸重復。本研究利用MISA 軟件分析了泡桐形成層及木質部區(qū)轉錄組，從長度為1 000 nt 以上的Unigenes 中獲得16 118 個SSR 位點，其中9 001 個是單核苷酸重復的SSR，接下來是5 024 個二核苷酸重復和2 048 個是三核苷酸重復SSR，這與連翹［17］、杜仲［18］和蓮霧［19］中發(fā)現(xiàn)的SSR 結果基本一致。這些結果將為泡桐SSR 標記開發(fā)提供重要序列信息，也為泡桐種質資源遺傳多樣性分析以及分子標記輔助育種等工作提供了重要基礎。