日本落葉松木質(zhì)部發(fā)育相關(guān)基因篩選及其共表達(dá)網(wǎng)絡(luò)構(gòu)建*

李 慧 代新仁 周再知 李全梓

(1.中國林業(yè)科學(xué)研究院林木遺傳育種國家重點(diǎn)實(shí)驗(yàn)室 北京 100091；2.中國林業(yè)科學(xué)研究院熱帶林業(yè)研究所 廣州510520)

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

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

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

1 材料與方法

1.1 材料收集及RNA提取

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

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

1.2 PacBio文庫構(gòu)建、測序及組裝

將木質(zhì)部、韌皮部和針葉組織樣品的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系統(tǒng)產(chǎn)生1～2 kb、2～3 kb、3～6 kb 和 5～10 kb的 cDNA片段。對于長度大于3 kb轉(zhuǎn)錄本，使用BluPinppin再篩選1次。產(chǎn)生的片段采用Pacific Biosciences SMRTbell template Prep Kit 1.0(part 100-259-100)試劑盒構(gòu)建文庫。用SMRT Cell 8 Pac v3(part100-171-800)在PcBio RSII real-time(RT)測序平臺上測序，每個(gè)樣品7個(gè)SMRT cells：其中1～2 kb和2～3 kb文庫各使用3個(gè)SMRT cells，大于3 kb的片段使用1個(gè)SMRT cell。使用SMRT Analysis Server(v2.3)獲得三代全長轉(zhuǎn)錄本序列：首先將獲得的插入片段進(jìn)行分類(Reads Classify)，該步驟可去掉來自SMRT cell cDNA分子的Reads of Insert cDNA引物和polyA尾，并將Reads of Insert分成全長(full-length)或非全長(non full-length)和嵌合(chimeric)或非嵌合(non-chimeric)類別；接著在LSC軟件中，用二代測序數(shù)據(jù)對全長非嵌合(FLNC)轉(zhuǎn)錄本進(jìn)行校正，具體方法參照Au等(2012)的研究；最后使用CD-HIT-EST(v4.5.4-2011-03-07)對不同片段文庫的FLNC進(jìn)行合并去冗余，最終獲得三代全長轉(zhuǎn)錄本序列(Lietal., 2006)。完成以上步驟后使用Cogent軟件對轉(zhuǎn)錄本進(jìn)行重構(gòu)，獲得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)數(shù)據(jù)庫得到注釋信息。通過NR數(shù)據(jù)庫中的注釋，用Blast2GO(https:∥www.blast2go.com)獲得GO數(shù)據(jù)庫(http:∥geneontology.org/)注釋(Conesaetal., 2005)。

1.3 mRNA文庫構(gòu)建、組裝及差異基因篩選

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

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

1.4 木質(zhì)部特異表達(dá)基因GO功能富集及KEGG通路富集分析

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

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

1.5 木質(zhì)部特異上、下調(diào)基因權(quán)重共表達(dá)網(wǎng)絡(luò)分析

通過R軟件中的WGNCA包對以上篩選的木質(zhì)部特異上、下調(diào)的基因進(jìn)行共表達(dá)網(wǎng)絡(luò)分析。對得到的不同模塊也進(jìn)行GO和KEGG分析，具體方法參照1.3。結(jié)合GO和KEGG 分析結(jié)果，篩選出木材形成相關(guān)的基因富集的模塊。參照以上篩選出的模塊中基因關(guān)系對文件中的權(quán)重系數(shù)，選排名前1 000的基因?qū)?，用Cytoscape中cytoHubba的12種算法：最大團(tuán)中心性(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)、發(fā)散性(Radiality，Rad)、中介性(Betweenness Centrality，BC)、應(yīng)力(Stress，Str)和集聚系數(shù)(Clustering coefficient，Cc)(Aokietal., 2007; Chinetal., 2014)，選出每種算法中排名前30的基因，12種算法中排名前30的基因的交集為本研究篩選出來的核心基因。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)控檢驗(yàn)

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

圖1 不同組織及相同組織不同重復(fù)間相關(guān)性分析及表達(dá)量統(tǒng)計(jì)

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

表1 3個(gè)組織中測序結(jié)果統(tǒng)計(jì)

2.2 木質(zhì)部特異上調(diào)、下調(diào)基因篩選

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

圖2 木質(zhì)部特異上調(diào)基因和下調(diào)基因篩選

2.3 木質(zhì)部特異上、下調(diào)基因GO和KEGG功能富集分析

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

圖3 木質(zhì)部特異基因GO和KEGG分析

2.4 核心基因與木材形成相關(guān)基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建及驗(yàn)證

為了解木質(zhì)部特異上調(diào)和下調(diào)基因內(nèi)部的共表達(dá)關(guān)系，對篩選出的基因進(jìn)行WGCNA分析。經(jīng)過30次迭代計(jì)算(圖4A)，2 596個(gè)木質(zhì)部特異上調(diào)和下調(diào)的基因被分為22個(gè)動(dòng)態(tài)模塊。經(jīng)合并后，最終得到3個(gè)合并模塊：深藍(lán)色模塊、深綠色模塊和綠色模塊(圖4B)。根據(jù)WGCNA相關(guān)理論(Lukensetal., 2012)，相關(guān)性越低，表明基因間獨(dú)立性越高。深藍(lán)色模塊及其中的基因與其他模塊及其中的基因相關(guān)性均偏低(圖4C、D)。結(jié)合模塊內(nèi)基因的KEGG及BLASTN結(jié)果，發(fā)現(xiàn)在深藍(lán)色模塊中木材形成相關(guān)的基因較多，為此選取該模塊內(nèi)的部分基因構(gòu)建共表達(dá)網(wǎng)絡(luò)。

圖4 木質(zhì)部特異上、下調(diào)基因WGCNA分析

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

為驗(yàn)證本研究構(gòu)建的共表達(dá)網(wǎng)絡(luò)中的共表達(dá)關(guān)系，將共表達(dá)網(wǎng)絡(luò)中的所有基因比對到毛果楊(Populustrichocarpa)基因組上。除Lkgene15627在毛果楊中無對應(yīng)同源基因，其余41個(gè)均可找到對應(yīng)同源基因。其中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個(gè)基因與毛果楊中的基因一一對應(yīng)(表2)。將此41個(gè)基因?qū)?yīng)的31個(gè)毛果楊同源基因輸入AspWood網(wǎng)站(http:∥aspwood.popgenie.org)驗(yàn)證共表達(dá)關(guān)系，結(jié)果表明28個(gè)同源基因在木質(zhì)部(≥300 μm)中均共表達(dá)(圖6C)，說明它們參與木質(zhì)部的形成，表明構(gòu)建的共表達(dá)網(wǎng)絡(luò)有效。

圖6 共表達(dá)網(wǎng)絡(luò)內(nèi)基因表達(dá)豐度及共表達(dá)關(guān)系驗(yàn)證

表2 日本落葉松與毛果楊的基因?qū)?yīng)

3 討論

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

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

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

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

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

4 結(jié)論

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