日本落葉松木質部發育相關基因篩選及其共表達網絡構建*

李 慧 代新仁 周再知 李全梓

(1.中國林業科學研究院林木遺傳育種國家重點實驗室 北京 100091；2.中國林業科學研究院熱帶林業研究所 廣州510520)

轉錄組學(transcriptome)是在RNA水平上檢測生物體中基因轉錄情況及其調控規律，并結合生物表型闡明基因功能的一門學科。隨著測序成本的降低及測序效率的提高，大規模并行測序技術(Illumina)RNA-seq成為當今發掘新基因和探明非模式生物基因表達模式的最為快捷、有力的手段(Luoetal., 2016; Songetal., 2015)。目前轉錄組測序及分析也廣泛應用于林木研究，在林木次生生長復雜分子機制研究中，利用該技術確定大量在木材形成過程中的相關候選基因及通路。最早的木材形成轉錄組研究是通過切向冷凍切片技術獲得不同木材形成階段的植物材料進行轉錄組測定，結合cDNA微矩陣芯片技術，最終得到2 995個探針，涵蓋了大約10%的楊樹(Populus)基因(Hertzbergetal., 2001)，隨后大量的相似研究獲得更高的基因覆蓋度(Moreauetal., 2009; 2005)。轉錄組技術的應用使得轉錄圖譜更加完善。通過該技術，在楊樹中確定了大量參與次生細胞壁成分(包括纖維素、木聚糖、葡甘露聚糖及木質素)形成過程的酶類編碼基因，且通過相關試驗證實這些基因確實在木材形成過程中發揮作用(Yeetal., 2015；Wangetal., 2019a；Yanetal., 2019)。針葉樹種的轉錄組研究起步相對較晚，最早是通過EST(Expressed Sequence Tag)測序(Kirstetal., 2003)，近幾年對針葉樹種也開展了RNA-seq方面的研究，例如Cronn 等(2017)對19株花旗松(Pseudotsugamenziesii)進行泛轉錄組測序，以此為參考，對179個針葉樣品二代測序結果進行比對、組裝，探討光周期和季節性變化對其生長和休眠的影響; Jokipii-Lukkari等(2018)按季節對歐洲云杉(Piceaabies)發育木質部、韌皮部和形成層進行轉錄組測序，用以揭示形成層的季節性變化規律及晚材的形成過程。

二代測序技術不僅能提供參考基因組，而且也能為序列重新組裝提供參考(Torreetal., 2019; Jokipii-Lukkarietal., 2018)，但是由于二代測序結果不完整及低質量的轉錄本限制了可變剪切版本的分析及剪切體的基因功能注釋(Metzker, 2010)。基于此類問題，Pacfic Biosciences 公司研發的 single-molecule real-time(SMRT)測序技術便應運而生(Eidetal., 2009; Carruthersetal., 2018; Duetal., 2017)，該測序技術能獲得全長或較長讀長的轉錄組，其中包括大量長讀長的轉錄本，平均讀長大于10 kb, 有些甚至可達60 kb，由此獲得的轉錄本包含完整的編碼序列及基因家族的共有特征(Seccoetal., 2013; Sharonetal., 2013; Tilgneretal., 2015)。但是SMRT-seq技術也存在缺陷，例如，由其提供的基因信息不太準確，低的基因覆蓋率導致較高錯誤率。因此結合二代測序技術(Illumina RNA-seq)及環形一致性序列(circular-consensus reads)對三代測序結果進行校正能獲得更為可靠的結果(Renetal., 2018; Wangetal., 2019b)。

針葉樹種基因組較大，一般可達10～40 Gb，重復序列較高(Chanoetal., 2017)，造成測序和組裝難度大；同時由于轉化體系不夠成熟及育苗周期長等因素的限制，導致基因功能研究無法順利開展。日本落葉松(Larixkaempferi)的生長速度在落葉松屬(Larix)中屬中等偏快，種植面積廣泛，且具有耐寒、樹干通直、適應性和抗病性強等特點(周賢武等，2014)，是我國主要的針葉紙漿用材樹種(孫曉梅等，2005)。目前關于日本落葉松轉錄組的研究僅見Zhang等(2012)對體細胞胚胎轉錄組的研究，以及Li等(2017b)對整個莖部轉錄組的研究。本研究以日本落葉松為研究材料，利用二代和三代轉錄組學結合生物信息學手段，對其木質部特異表達的基因進行篩選，通過GO、KEGG分析及與擬南芥(Arabidopsisthaliana)同源比對，獲得與木材形成相關的基因；通過WGCNA分析，篩選出與木材形成相關基因共表達的核心基因，為開展日本落葉松及其他針葉樹種基因功能研究提供參考。

1 材料與方法

1.1 材料收集及RNA提取

日本落葉松材料取自遼寧省清原縣大孤家鎮種質植物園(42°22′N, 124°51′E)。試材植株7年生，取材時先將樹皮剝掉，輕刮樹皮內側和樹干外層除去形成層，從去形成層樹干外層刮取木質部，從去形成層樹皮內層刮取韌皮部，選取健康、成熟針葉。材料收集后立即包入無RNAase的錫紙后置于液氮中速凍保存。每組織設置3個生物學重復。

總RNA的提取采用Lorenz等(2010)的改良CTAB法。經DNaseI去除DNA，無水乙醇沉淀，RNase-free水溶解，得到的RNA經Agilent 2100 Bioanalyzer和NanoDrop測定RNA濃度、RIN值、28S/18S比值、評估片段大小及OD260/OD280值等指標后，確定RNA的完整度和純度，保存于-80 ℃超低溫冰箱備用。

1.2 PacBio文庫構建、測序及組裝

將木質部、韌皮部和針葉組織樣品的RNA等量混樣后，取1 μg 混樣RNA，先通過Clontech SMARTer PCR cDNA synthesis Kit(cat.no.634926)試劑盒的鉚釘引物olig(dT)30合成第1鏈cDNA。通過Advantage 2 PCT kit(Clonetech, cat.No.639206)試劑盒的long PCR(LD-PCT)合成第2鏈。采用BluePinppin size selection系統產生1～2 kb、2～3 kb、3～6 kb 和 5～10 kb的 cDNA片段。對于長度大于3 kb轉錄本，使用BluPinppin再篩選1次。產生的片段采用Pacific Biosciences SMRTbell template Prep Kit 1.0(part 100-259-100)試劑盒構建文庫。用SMRT Cell 8 Pac v3(part100-171-800)在PcBio RSII real-time(RT)測序平臺上測序，每個樣品7個SMRT cells：其中1～2 kb和2～3 kb文庫各使用3個SMRT cells，大于3 kb的片段使用1個SMRT cell。使用SMRT Analysis Server(v2.3)獲得三代全長轉錄本序列：首先將獲得的插入片段進行分類(Reads Classify)，該步驟可去掉來自SMRT cell cDNA分子的Reads of Insert cDNA引物和polyA尾，并將Reads of Insert分成全長(full-length)或非全長(non full-length)和嵌合(chimeric)或非嵌合(non-chimeric)類別；接著在LSC軟件中，用二代測序數據對全長非嵌合(FLNC)轉錄本進行校正，具體方法參照Au等(2012)的研究；最后使用CD-HIT-EST(v4.5.4-2011-03-07)對不同片段文庫的FLNC進行合并去冗余，最終獲得三代全長轉錄本序列(Lietal., 2006)。完成以上步驟后使用Cogent軟件對轉錄本進行重構，獲得UniTransModels(Lietal., 2017a)，并將其比對到NR(ftp:∥ftp.ncbi.nlm.nih.gov/blast/db)、NT(ftp:∥ftp.ncbi.nlm.nih.gov/blast/db)、COG(http:∥www.ncbi.nlm.nih.gov/COG)、KEGG(http:∥www.genome.jp/kegg)和Swiss-Prot(http:∥ftp.ebi.ac.uk/pub/databases/swissprot)數據庫得到注釋信息。通過NR數據庫中的注釋，用Blast2GO(https:∥www.blast2go.com)獲得GO數據庫(http:∥geneontology.org/)注釋(Conesaetal., 2005)。

1.3 mRNA文庫構建、組裝及差異基因篩選

取檢測合格的總RNA 1 μg，用含有Oligo(dT)的磁珠富集mRNA，并用片段化緩沖液將其打成小片段。以此mRNA為模板，采用六聚體隨機引物合成cDNA第1鏈，再加入緩沖液、dNTPs、RNaseH和DNApolymeraseI合成第2鏈cDNA。得到的cDNA文庫經純化后加入接頭A。經PCR擴增后，產物用Illumina HiSeq 2500測序平臺進行測序。最終得到100 bp的成對讀長。

測序完成后，利用NGS QC Toolkit 過濾原始數據得到待分析數據。再用Hierarchical Indexing for Spliced Alignment of Transcripts(HISAT)(Kimetal., 2010)將二代測序數據匹配到三代測序所得轉錄本上。使用Cufflinks 軟件計算每個樣本中所有基因的FPKM值。對不同重復及樣品進行Pearson相關性檢驗。通過DEseq2(用來分析差異表達基因的R軟件包)獲得不同組織間的差異表達基因。

1.4 木質部特異表達基因GO功能富集及KEGG通路富集分析

為了解篩選出來的木質部特異的高表達和低表達基因在生物體內的功能，本研究利用Goseq和topGO(Youngetal., 2010)進行基因本體論(Gene Ontology，GO)富集分析。通過將所有mRNA比對到GO分析數據庫(http:∥geneontology.org/)中的GO terms, 統計富集到每個GO條目中的mRNA數目并進行超幾何檢驗，將其比對到參考基因背景上。同時使用京都基因與基因組百科全書數據庫(Kyoto Encyclopedia of Genes and Genomes，KEGG)(https:∥www.kegg.jp/)對木質部特異高表達和低表達的mRNA進行代謝及信號轉導途徑富集分析并通過超幾何檢驗將其比對到整個參考基因背景。利用R軟件中的phyper計算GO和KEGG的富集顯著性P，計算公式如下：

式中：N為具有GO和KEGG注釋所有背景基因的數目；n為N中差異表達基因的數目；M為所有基因中注釋到某特定GO 條目和KEGG代謝通路中的基因的數目；m為注釋到某特定GO 條目和KEGG代謝通路中的差異表達基因的數目。

1.5 木質部特異上、下調基因權重共表達網絡分析

通過R軟件中的WGNCA包對以上篩選的木質部特異上、下調的基因進行共表達網絡分析。對得到的不同模塊也進行GO和KEGG分析，具體方法參照1.3。結合GO和KEGG 分析結果，篩選出木材形成相關的基因富集的模塊。參照以上篩選出的模塊中基因關系對文件中的權重系數，選排名前1 000的基因對，用Cytoscape中cytoHubba的12種算法：最大團中心性(Maximal Clique Centrality，MCC)、最大鄰域分量密度(Density of Maximum Neighborhood Component，DMNC)、最大鄰域組成(Maximum Neighborhood Component，MNC)、度值法(Degree method，Deg)、邊緣滲流分量(Edge Percolated Component，EPC)、瓶頸值(BottleNeck，BN)、偏心度(EcCentricity，EC)、緊密度(Closeness，Clo)、發散性(Radiality，Rad)、中介性(Betweenness Centrality，BC)、應力(Stress，Str)和集聚系數(Clustering coefficient，Cc)(Aokietal., 2007; Chinetal., 2014)，選出每種算法中排名前30的基因，12種算法中排名前30的基因的交集為本研究篩選出來的核心基因。

2 結果與分析

2.1 測序數據質控檢驗

對不同組織、不同重復間基因的表達量進行pearson相關性檢驗(圖1A)。同一組織不同重復間的相關性較高，相關系數均大于0.9，表明試驗中3個重復間的重復性較好；而不同組織間的相關系數較低，表明不同組織間的基因表達存在較大差異。可見，本試驗選取3個組織的3個生物學重復滿足數據分析基本要求。對木質部、韌皮部和針葉中基因的表達水平分布統計結果(圖1B)顯示，3個組織標準化表達量的峰值在0.1～0.2處較高，其中表達量中位數位于0.4～0.5之間。表達量中位數值在韌皮部中最高，在針葉中最低。

圖1 不同組織及相同組織不同重復間相關性分析及表達量統計

對日本落葉松的的二代轉錄組測序分析(表1)顯示，木質部、韌皮部和針葉中分別得到66 463 409、77 166 823 和 63 449 442條測序總讀長。將原始讀長比對到已校正的三代測序文庫，其中在木質部、韌皮部和針葉中能比對上的測序讀長分別為24 451 888、31 954 781和34 865 468，經計算3個組織的比對率分別為36.79%、41.41%和54.95%。最終經過拼接和合并轉錄本，得到在木質部中表達的基因21 535條，在韌皮部表達的基因22 986條，在針葉中表達的基因23 233條。

表1 3個組織中測序結果統計

2.2 木質部特異上調、下調基因篩選

為了解不同組織基因的表達差異，對3個組織在設定閾值log2(FC)≥1或log2(FC)≤-1，Q≤ 0.05下進行兩兩比較得到差異表達基因。木質部(X)相對于韌皮部(P)得到4 836個差異表達基因，木質部相對于針葉(L)得到6 213個差異表達基因(圖2)。通過log2(X vs P)≥1 且Q≤ 0.05篩選閾值篩選木質部相對韌皮部高表達的基因，通過log2(X vs L)≥1 且Q≤ 0.05篩選閾值篩選木質部相對針葉高表達基因，二者交集篩選木質部特異上調基因。通過log2(X vs P)≤-1 且Q≤ 0.05篩選木質部相對韌皮部低表達基因，通過log2(X vs L)≤-1且Q≤ 0.05篩選木質部相對針葉低表達基因，二者交集篩選木質部特異表達下調基因。最終共得到1 648個木質部特異表達上調基因和948個木質部特異表達下調基因(圖2)。

圖2 木質部特異上調基因和下調基因篩選

2.3 木質部特異上、下調基因GO和KEGG功能富集分析

為了解篩選出的2 596個木質部特異表達基因的功能，用GO和KEGG進行功能和代謝通路預測。GO分析主要集中在生物學過程(biological process)、細胞成分(cellular component)及分子功能(molecular function)3方面。對得到的初步結果，在每個GO分類中選擇排名前5的GO term進行可視化。本研究篩選出的2 596個木質部特異表達上調和表達下調的差異基因富集于39個GO 類別中，其中在生物學過程中，排名前3的GO 類別為代謝過程、細胞過程和定位，富集的基因數目分別為900、893和204(圖3A)；在細胞成分中，排名前3的GO 類別為膜、細胞和細胞組件(cell part)，富集的基因數目分別為615、561和557；在分子功能中，排名前3的GO類別為催化活性、結合位點(binding)和轉運活性，富集的基因數目分別為1 081、812和113。在KEGG分析結果中，以Q從小到大選擇排序前15的代謝通路進行可視化。2 596個基因富集于112個代謝通路中，排名前3的代謝通路分別為淀粉和蔗糖代謝(starch and sucrose metabolism)、類黃酮生物合成(flavonoid biosynthesis)及代謝途徑(metabolic pathways)(圖3B)，富集因子分別為0.196、0.431 和0.110。其中本研究關注的苯丙烷代謝途徑、淀粉和蔗糖代謝途徑富集的基因數目分別為38和196。

圖3 木質部特異基因GO和KEGG分析

2.4 核心基因與木材形成相關基因共表達網絡的構建及驗證

為了解木質部特異上調和下調基因內部的共表達關系，對篩選出的基因進行WGCNA分析。經過30次迭代計算(圖4A)，2 596個木質部特異上調和下調的基因被分為22個動態模塊。經合并后，最終得到3個合并模塊：深藍色模塊、深綠色模塊和綠色模塊(圖4B)。根據WGCNA相關理論(Lukensetal., 2012)，相關性越低，表明基因間獨立性越高。深藍色模塊及其中的基因與其他模塊及其中的基因相關性均偏低(圖4C、D)。結合模塊內基因的KEGG及BLASTN結果，發現在深藍色模塊中木材形成相關的基因較多，為此選取該模塊內的部分基因構建共表達網絡。

圖4 木質部特異上、下調基因WGCNA分析

根據WGCNA分析的基因對權重系數，對基因模塊由高到低排序，選出排名前1 000的基因對。通過Cytoscape中cytoHubba的12種算法分別篩選排名前30的核心基因，取交集后即為本模塊中的核心基因。經初步篩選，共獲得18個核心基因(圖5)。將深藍色模塊的邊權重文件中排名前50 000的基因關系對導入Cytoscape中，構建核心基因與木材形成相關基因共表達網絡(圖6A)。通過建立網絡，最終在深藍色模塊中篩選出17個核心基因，這17個基因與纖維素、半纖維素合成相關的11個基因和木質素合成相關的14個基因都具有較高的關聯度。網絡中的42個基因中除Lkgene1125和Lkgene4328外，其余基因在木質部中表達量均高于在韌皮部和針葉中的表達量(圖6B)。

為驗證本研究構建的共表達網絡中的共表達關系，將共表達網絡中的所有基因比對到毛果楊(Populustrichocarpa)基因組上。除Lkgene15627在毛果楊中無對應同源基因，其余41個均可找到對應同源基因。其中Lkgene2505和Lkgene1394與毛果楊中的Potri.003G183900同源；Lkgene540和Lkgene816與毛果楊中Potri.006G033300同源；Lkgene979、Lkgene560和Lkgene961與毛果楊中的Potri.004G089300同源；Lkgene8640、Lkgene8812、Lkgene6881、Lkgene19434和 Lkgene28309與毛果楊中的Potri.005G145300同源；Lkgene166和Lkgene231與毛果楊中的Potri.013G157900同源；Lkgene18758和Lkgene27841與毛果楊中的Potri.002G202300同源。其余25個基因與毛果楊中的基因一一對應(表2)。將此41個基因對應的31個毛果楊同源基因輸入AspWood網站(http:∥aspwood.popgenie.org)驗證共表達關系，結果表明28個同源基因在木質部(≥300 μm)中均共表達(圖6C)，說明它們參與木質部的形成，表明構建的共表達網絡有效。

圖6 共表達網絡內基因表達豐度及共表達關系驗證

表2 日本落葉松與毛果楊的基因對應

3 討論

目前已在擬南芥和楊樹中鑒定出大量木質素和纖維素形成相關的基因(楊立, 2016; Xieetal., 2011; Zhongetal., 2008)。而日本落葉松由于還無參考基因組，相關方面研究仍相對滯后。在Li等(2017b)日本落葉松的莖部轉錄組的研究中鑒定到EXLA1、EXPB17、LRX3、BRU1和CSLD3等參與木材形成的基因。本次研究鑒定到Lkgene10000與擬南芥的CSLA9同源，參與半乳甘露聚糖生物合成代謝(表3)，同Li等(2017b)鑒定到的CSLD3基因均屬于CSL(cellulose synthase-like)超基因家族成員。研究表明在CSL家族中存在6個亞型家族，包括CslA、CslB、CslC、CslD、CslE和CslG(Somervilleetal., 2000)。在植物細胞壁中纖維素主要由CSLA基因家族成員編碼相關纖維素合成酶(徐宗昌等，2017)。CSL家族亞型內部成員功能也有差異，例如，在楊樹中發現的PtCslA1、PtCslA3和PtCslA5，其中PtCslA1和PtCslA3可以編碼甘露聚糖合成酶，PtCslA1能進一步編碼葡甘聚糖合成酶，而PtCslA5對木聚糖或葡甘聚糖底物無任何催化活性(Suzukietal., 2006)。Lkgene25783、Lkgene11296和Lkgene7496與擬南芥的F22D1.120同源，預測參與木葡聚糖生物合成過程，半乳甘露聚糖合成代謝和木葡聚糖合成代謝均與植物體內半纖維素合成相關。本研究還鑒定出纖維素合成相關基因，其中Lkgene946、Lkgene1213和Lkgene1637與擬南芥的DEC同源，Lkgene5019與CEL1同源，二者均為內切葡聚糖酶(endoglucanase)家族成員，該家族成員可參與纖維素微纖絲、細胞壁纖維素合成和次生細胞壁形成過程(Nicoletal., 1998)。6個基因(表3)為蔗糖磷酸合成酶(sucrose phosphate synthase)家族成員，與擬南芥的SPS3基因同源，可調控細胞中纖維伸長過程(Lutfiyyaetal., 2007)。

6個基因與擬南芥中的PAL4同源(表3)，主要參與苯丙烷代謝的第1步：由苯丙氨酸到肉桂酸的轉化過程。苯丙氨酸在PAL的脫氨基作用下形成反式肉桂酸(trans-cinnamic acid)及氨基離子。同時PAL是連接苯丙氨酸代謝途徑和芥草酸的重要催化酶，抑制PAL可減少木質素單體合成物的前體底物(楊立，2016)。Lkgene3952與擬南芥的CCR1同源，在木質素合成的后期起到關鍵作用。另外，在毛果楊中有報道表明，PtrCCR2與PtrCAD1可形成異源二聚體來提高毛果楊木質素單體生物合成過程中的酶活性進而促進木質素單體的形成(Wangetal., 2019a; Yanetal., 2019)。而在日本落葉松中，是否也存在該類蛋白復合體調控木質素單體合成還需進一步試驗證實。Lkgene5584、Lkgene540、Lkgene826和Lkgene816與擬南芥的CYP98A3同源，為細胞色素催化單氧酶基因家族成員，能夠催化對香豆酯(p-coumaric esters)3′端羥基化，進而形成木質素單體(Schochetal., 2001)。Lkgene19432、Lkgene4199和 Lkgene8671為PER(peroxidase)家族基因，參與木質素的生物合成和降解過程。近來在對小麥(Triticumaestivum)的過氧化物酶活性試驗中發現，該類基因同時可能參與細胞壁中半纖維素的合成(Lorenzoetal., 2019)。Lkgene1782與擬南芥中的HXK2基因同源，參與細胞程序化死亡(Kimetal., 2006)。

本研究鑒定到C4H、HCT和COMT3類可影響S型木質素單體的含量的基因。其中Lkgene231和 Lkgene166與擬南芥的C4H同源，可調控合成木質素的基礎代謝——碳代謝流(表3)。C4H為血紅素P450細胞色素催化單氧酶，在生物反應中可充當還原劑，肉桂酸經C4H的作用形成對香豆酸(p-coumaric acid)。相關研究表明，C4H可參與由內硅酸到羥基苯丙烯酸的生物合成過程。在煙草(Nicotiana)中，抑制C4H能降低植株中木質素含量且影響S型與G型木質素的比例(Kajitaetal.,1997)。Lkgene1394和Lkgene2505與擬南芥中的HCT基因同源，參與木質素合成過程(表3)。HCT可與C3H協同催化由香豆酰輔酶A到咖啡酰輔酶A的轉化，能控制3種木質素單體在生物體內的比例分配，沉默HCT可增加H型木質素，降低S型木質素含量(楊立，2016)。Lkgene3070、Lkgene834和Lkgene8363與COMT1同源，可調控木質素單體甲基化過程。COMT是位于F5H催化反應下一步的咖啡酸-O-甲基轉移酶，其可以5-羥基松柏醇、5-羥基松柏醛和咖啡酸為底物生產芥子酸、芥子醛和阿魏酸。抑制COMT表達時可降低S型木質素含量(Osakabeetal., 1999)。在本研究中，鑒定到的C4H、HCT和COMT1在日本落葉松的木質部中表達量都較高(圖6B)，暗示在日本落葉松木質部中S型木質素含量較高。

核心基因在整個后轉錄過程中起到關鍵作用，它們可能是整個調控網絡中的關鍵樞紐。在本研究中鑒定到36個基因與擬南芥中的LAC12同源(表3)，其中Lkgene26047被確認為核心基因(圖6)，預測可能參與木質素降解過程。目前關于LAC12在植物中功能的報道較少，僅在酵母中證實參與乳糖和半乳糖的轉運過程(Rigamonteetal., 2011)。23個基因與擬南芥的LAC17同源(表3)，其中Lkgene10602被確認為核心基因(圖6)，可參與木質素單體的聚合過程和維管束間的G型木質素沉降過程。LAC17可參與木質素合成，當抑制其表達時可顯著降低莖部木質素含量(Berthetetal., 2011)。4個基因與擬南芥的CTL2同源(表3)，其中Lkgene13943被確認為核心基因(圖6)，參于纖維素的合成過程及微纖絲和半纖維素的形成過程。CTL2為幾丁質酶家族成員，功能同CTL1/POM1復合體相近，可調控纖維素的組裝過程，并通過與微纖絲結合與半纖維素發生作用(Sanchezetal., 2012)。同時有研究表明LAC17和CTL2與CESA基因在擬南芥中均存在共表達關系(Brownetal., 2005; Perssonetal., 2005)，但在本研究構建的網絡中并未發現CESA基因，可能與構建網絡的閾值設定或物種本身的差異有關。基因高度共表達表明基因之間可能有共同的調控機制或相似功能(Alloccoetal., 2004)，盡管其他14個核心基因是否參與木材形成過程目前尚無報道，但是根據共表達網絡結果，推測其也可能參與木材形成的某些過程，仍需試驗進一步確定其功能。

本研究由于基于三代測序進行二代轉錄組比對，木質部、韌皮部和針葉中的比對率只有36.79%、41.41%和 54.95%，與其他有參物種相比，比對率相對較低。同時由于數據庫注釋信息局限，導致部分基因注釋信息不完整，無法確定其功能。盡管早有報道證實物種間序列有一定的保守性，但不同物種的同源基因也有其自身特點，例如在楊樹中基因PTLF不僅在花絮中表達，同時在幼葉和葉原基處也表達，但在擬南芥中僅發現在花中表達(Rottmannetal., 2000)。因此本研究基于擬南芥為參考的功能推測可能導致對基因功能認識不全面。

4 結論

本研究以三代混樣測序為參考對二代測序結果進行組裝、比對，通過木質部相對于韌皮部和針葉的差異基因比較，篩選出木質部特異上調基因和特異下調基因。通過GO、KEGG和BLASTN對篩選的木質部特異上、下調基因進行基因功能預測。研究結果表明篩選的日本落葉松木質部發育相關基因參與半乳甘露聚糖合成、木葡聚糖合成、纖維素微纖絲形成、細胞壁纖維素合成、次生細胞壁形成過程、纖維伸長過程、調控合成木質素的碳代謝流、木質素生物合成及降解、木質素單體聚合、木質素單體甲基化和細胞程序化死亡等木材形成相關生物學過程。利用WGCNA分析在共表達網絡中鑒定出17個與木材形成相關基因關系緊密的核心基因，預測其在木材形成過程中發揮重要作用。在當前無參考基因組的情況下，本研究為探討針葉樹種木材形成機理、挖掘基因功能提供有力依據。