李 慧 代新仁 周再知 李全梓

(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个与木材形成相关基因关系紧密的核心基因，预测其在木材形成过程中发挥重要作用。在当前无参考基因组的情况下，本研究为探讨针叶树种木材形成机理、挖掘基因功能提供有力依据。