文冠果(Xanthoceras sorbifolia Bunge)属无患子科文冠果属，广泛分布于我国东北、华北、西北地区，种子含油率为30%～40%，籽油适合生产生物柴油，是我国特有的油料树种，具有结实早、易繁殖、适应性强、耐旱、耐寒、耐瘠薄等特点。《全国林业生物质能源发展规划(2011~2020年)》将文冠果作为重点发展树种。大力发展文冠果对改善生态环境、缓解能源短缺局面具有重要意义^[1]。国内外文冠果研究主要集中在栽培和常规育种等方面，对其分子生物学研究相对较少^[2]。基于第二代高通量测序转录组测序(RNA-seq)技术是一种高效、快捷的获得基因序列的研究手段，利用该手段可消除对基因组序列的依赖，从整体水平了解植物在特定阶段的基因功能及基因结构，便于揭示特定生物学过程的分子机制^{[3, 4, 5]}。近年来，已有紫杉、桉树、日本落叶松、橡胶树、茶树等多种植物进行了转录组测序^{[6, 7, 8, 9, 10]}。

本研究利用高通量测序技术平台Illumina Solexa Hi-seq^TM 2000对文冠果花芽进行转录组测序。拼接组装后对得到的unigene进行功能注释和功能分类，还对文冠果转录组序列潜在的SSR位点进行发掘，为功能基因的发掘利用及功能分析等奠定基础。

文冠果材料采集于河北省承德市(E117°50′，N40°58′)，该地区属于低山丘陵阳坡厚层土立地类型。林分为文冠果纯林，林龄50年，平均树高5.11m，平均冠幅19.12m²。于3月末花芽分化临界期采集花芽，样品取自10个单株。样品采集后立即放入液氮中，之后存放于－80℃冰箱备用。

总RNA的提取采用RNeasy Plant Mini Kit(Qiagen,Inc.,Valencia，CA，USA)试剂盒，按照实验手册提取。利用Agilent 的BioAnalyzer 2100 检测RNA完整性，样品RIN≥8。取10μg用于cDNA文库构建。用带有Oligo(dT)的磁珠富集mRNA，然后将得到的mRNA处理成小片段。以mRNA小片段为模板，用随机引物反转录生成dsDNA，末端修复后，在3′端加A，并连接接头，进行片段大小选择，最后进行PCR扩增，构建好的文库用Illumina Hi-seq^TM 2000的paired-end测序方法进行转录组测序。

鉴于Solexa数据错误率对结果的影响，对原始数据进行质量预处理。去除含有adaptor的read；去除N比例大于5%的read；去除低质量read(质量值Q≤10的碱基数占整个read的20%以上)；获得clean read，后续分析都基于clean read进行。然后利用软件Trinity对有效read进行从头拼接^[11]。将具有重叠区域的read组装成长片段(contig)，然后将这些contig连在一起组装成两端不能再延长的序列(unigene)。

使用BLASTx程序将拼接所得unigene与核酸、蛋白质序列数据库比对(比对参数为E值＜0.000 01)，并选取最佳注释。蛋白质数据库包括GenBank非冗余蛋白数据库(Nr)、Swiss-Prot(Swiss-Protein)、京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)、蛋白质直系同源簇数据库(clusters of orthologous groups of proteins，COGs)和Gene Ontology(GO)。其中，unigene通过COG、GO和KEGG数据库的分类参考Wang等^[12]的方法。

利用Microsatellite(MISA)软件在所有unigene中搜索SSR位点，参数设置如下：单核苷酸、二核苷酸、三核苷酸、四核苷酸、五核苷酸和六核苷酸最少重复次数分别为12、6、5、5、4和4。复合SSR两个位点间最大间隔碱基数为100。

采用Illumina Hi-seq^TM 2000 测序技术对文冠果花芽进行转录组测序，共获得56 632 980条原始序列。经过预处理，最终得到有效序列52 869 962条。使用软件Trinity对有效read进行de novo组装，得到91 946条序列重叠群(contig)，总长度约为38 742 364bp，平均长度及N50分别为421bp和1 060bp。取每个contig下最长的转录本作为unigene，得到58 311条unigene，总长度约为40 016 117bp，平均长度与N50分别为686bp和1 180bp。其中，11 233条(19.26%)unigene长度大于500bp，6 257条(10.73%)unigene长度大于1 000bp，3 610条(6.19%)unigene长度大于1 500bp，长度大于2 000bp的序列共有3 311条，占unigene总数的5.68%。最长的unigene长度为17 698bp(图 1)。这些结果说明测序和组装效果很好，可以进行基因功能分析。

图 1 unigene和contig的长度及百分比分布 Fig. 1 Length and percent distribution of unigene and contig

图选项

将所获得的58 311条unigene与公共数据Nr、NT、Swiss-Prot、KEGG库进行Blastx比对，通过gene的相似性进行功能注释，共有37 047条unigene获得了基因注释，占unigene总数的63.53%(表 1)。另有21 264条unigene(36.47%)未被注释，认为可能是新基因。

unigene与Nr数据库比对(E值＜0.000 01)，有8 665条(24.58%)与葡萄(Vitis vinifera)的序列同源；7 139条(20.25%)与蓖麻(Ricinus communis)序列同源；6 217条(17.63%)与毛果杨(Populus trichocarpa)序列同源；5 331条(15.12%)与碧桃(Prunus persica)序列同源(图 2)。与Swiss-Prot数据库比对，有22 176条(38.03%)文冠果unigene 找到了同源序列。

图 2 unigene的物种分布 Fig. 2 Species distribution for unigene

图选项

对注释到unigene的E值进行统计，其中E值为0的unigene占总unigene总数的14.63%(5 158条)；E值小于10^－100的unigene共有6 466条，占unigene总数的18.34%(图 3)。

图 3 unigene的E值分布 Fig. 3 E-value distribution for unigene

图选项

为进一步研究文冠果中unigene功能分类，将本研究将所得的58 311条unigene在COG和GO数据库中进行比对及功能注释与分类。

与COG数据库比对(E值≤1e -10)，58 311条文冠果unigene中，11 637条(19.96%)被注释到25个COG类别中(图 4)。其中，“一般功能基因”(general function prediction only)是最大类别，包含3 645条unigene，占被注释unigene总数的31.32%；其次是“转录”(transcription)，包含1 928条(16.57%)unigene，其余依次为“复制、重组和修复”(replication，recombination and repair)(1 846条，15.86%)，“转录后修饰、蛋白质翻转和分子伴侣类”(posttranslational modification，protein turnover，chaperones)(1 587条，13.64%)。而“核酸结构”(nuclear structure)是最小的类别，仅包含1条unigene。

图 4 文冠果unigene的COG分类 Fig. 4 Classification of COG for yellow-horn 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 transport 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

图选项

此外，文冠果unigene的GO注释分类结果表明(图 5)，共有27 514条unigene被注释到生物学过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)三大类中。其中，110 232条unigene归入生物学过程，82 734条unigene归入细胞组分，32 036条unigene归入分子功能。三大类又被划分为55小类。其中细胞、细胞组分、细胞过程、细胞器、代谢过程几类富集程度较高，这可能与花芽分化过程中分生组织不断进行细胞增殖来建成各种花器官原基，导致花芽内代谢活动旺盛有关。生物学过程包括22类，细胞过程(cellular process)是其中最大类别，包含17 164条unigene，占unigene总数的29.44%；细胞组分包括17类，细胞(cell)和细胞部分(cell part)是其中的最大类别，均包含20 883条unigene，占unigene总数的35.81%；分子功能包括16类unigene，催化活性(catalytic activity)是其中最大类别，包含13 519条unigene，占unigene总数的23.18%。

图 5 文冠果unigene的GO分类 Fig. 5 GO classification for yellow-horn B1：Biological adhesion；B2：Biological regulation；B3：Cellular component organization or biogenesis；B4：Cellular process；B5：Developmental process；B6：Establishment of localization；B7：Growth；B8：Immune system process；B9：Localization；B10：Locomotion；B11：Metabolic process；B12：Multi-organism process；B13：Multicellular organismal process；B14：Negative regulation of biological process；B15：Positive regulation of biological process；B16：Regulation of biological process；B17：Reproduction；B18：Reproductive process；B19：Response to stimulus；B20：Rhythmic process；B21：Signaling；B22：Single-organism process；C1：Cell；C2：Cell junction；C3：Cell part；C4：Extracellular matrix；C5：Extracellular matrix part；C6：Extracellular region；C7：Extracellular region part；C8：Macromolecular complex；C9：Membrane；C10：Membrane part；C11：Membrane-enclosed lumen；C12：Nucleoid；C13：Organelle；C14：Organelle part； C15：Symplast；C16：Virion；C17：Virion part；M1：Antioxidant activity；M2：Binding；M3：Catalytic activity；M4：Channel regulation activity；M5：Electron carrier activity；M6：Enzyme regulator activity；M7：Metallochaperone；M8：Molecular transducer activity；M9：Nucleic acid binding transcription factor activity；M10：Nutrient reservoir activity；M11：Protein binding transcription factor activity；M12：Protein tag；M13：Receptor activity；M14：Structural molecule activity；M15：Translation regulator activity；M16：Transporter activity

图选项

将文冠果unigene序列映射到KEGG数据库的参考代谢通路(pathway)中，共有19 556条unigene被注释，涉及128个代谢通路。其中包含unigene最多的代谢通路是代谢途径(ko01100，metabolic pathways)，共有4 208条(21.52%)unigene。其次是次生代谢产物的生物合成(ko01110，biosynthesis of secondary metabolites)，包含2 150条(10.99%)unigene(表 2)。Pathway富集性分析的前10个pathway数据列于表 2。

利用Perl操作平台下的MISA软件，在58 311条文冠果花unigene序列中筛选SSRs位点，其中含不同重复基元SSRs的序列有9 794条，共12 213个SSRs，包含SSR的一致序列出现频率为16.66%，平均每3.28kb出现1个SSR，而包含有2个及2个以上SSR的unigene共有1 866条。在12 213个SSR中共有369种重复基元(motif)。其中，1、2、3、4、5和6碱基重复基元出现频率进行统计，最多的重复基元分别是(A/T)n、(AG/CT)n、(AAG/CTT)n、(AAAG/CTTT)n、(AAAT/ATTT)n、(AAAAG/CTTTT)n、(AAAAAC/GTTTTT)n、(AGAGGG/CCCTCT)n。它们在各自重复基元类型中的比例分别是99.88%、75.59%、31.99%、16.05%、16.05%、17.88%、4.26%和4.26%。在所有类型的重复基元中，单核苷酸重复基元出现的频率最高，为34.95%，其次分别为二核苷酸(32.74%)、三核苷酸(28.64%)、六核苷酸(1.54%)、五核苷酸(1.47%)和四核苷酸(0.66%)重复基元(表 3)。除此之外，不同核苷酸的重复次数也有很大变化。对这些SSR的鉴定，将为进行文冠果的基因组差异性分析、遗传图谱构建等研究提供帮助。

植物成花诱导主要有4种途径：光周期途径、春化途径、赤霉素途径和自主途径^[13]。光周期途径相关基因中，本研究发掘出29条光敏色素基因PHYA同源序列，35条PHYB基因同源序列，隐花色素基因CRY1有17条同源序列，CRY2有7条同源序列。春化途径相关基因中，LEO1基因同源序列31条，ACT1基因同源序列84条。赤霉素途径(GA 途径)相关基因GAI同源序列277条，RGL1同源序列39条，RGA同源序列180条，DELLA1同源序列38条。自主途径相关基因，如FLD同源序列13条、FVE同源序列3条。此外，从转录组数据中还获得了与年龄调节相关的SPL基因同源序列3 936条。

Illumina高通量测序不需要知道物种遗传背景，且具有高通量、成本低、灵敏度高、重复性好等优点，已成为转录组研究的主要手段^{[14, 15, 16]}，适于文冠果等没有开展全基因组测序的物种开展转录组测序研究。本研究采用Illumina高通量测序技术对文冠果花芽形成期转录组进行研究。共拼接得到58 311条unigene，有37 047条unigene获得基因注释，另有21 264条unigene(36.47%)未被注释。文冠果为单属单种，推测这些未注释到的unigene可能是文冠果组织特异性新基因，但也可能是短的非编码序列。这表明高通量测序技术是批量发现文冠果功能基因的有效手段。

SSR标记已被广泛应用于分子标记辅助选择育种，通过关联分析，发掘与目的性状连锁紧密的相关基因，可大大缩短多年生植物育种年限^[17]。本研究发掘到12 213个SSR位点，SSR出现频率为0.305 SSR/kb，平均每3.28kb出现1个SSR。微卫星出现频率与桉树^{[18, 19]}、杏^[20]等相近。并发现文冠果以小于20bp的短微卫星重复序列最多，大于20bp的长微卫星占微卫星总数的13.75%。据此推断文冠果所含的微卫星可能受到较强趋同选择的压力，使其富集在较短的序列范围内。下一步可对这些引物进行扩增检测，筛选出目的条带清晰、多态性好的引物，为群体遗传多样性分析和构建遗传连锁图谱等后续研究及分子育种奠定基础。本研究中二核苷酸出现的频率为32.74%，略高于三核苷酸(28.64%)。与麻疯树(Jatropha curcas)^[21]、橡胶树(Hevea brasiliensis)^[22]、银杏(Ginkgo biloba)^[23]中均是以二核苷酸重复单元较多的现象一致。为保证SSR位点潜在多态性，我们在筛选过程中对于四核苷酸、五核苷酸和六核苷酸的最小重复次数设置为5、4、4，一定程度上影响到了这3类核苷酸重复在总SSR位点中所占的比例。

植物成花诱导的光周期途径中，光敏色素主要吸收红光和远红光，隐花色素主要吸收蓝光和紫外光^[24]。在拟南芥中共发现了4类光敏色素PHYA、PHYB/D、PHYC/F、PHYE，以及3种隐花色素CRY1、CRY2和CRY3，它们感受昼夜长短和光的强弱，产生昼夜节律，进而启动或抑制开花进程^[25]。本研究发掘出29条PHYA同源基因，35条PHYB同源序列，CRY1和CRY2也都发现了同源序列，可见这些基因的保守性。

春化途径相关基因中本研究发掘出LEO1基因同源序列31条，ACT1基因同源序列84条。但没有发掘到VIN1、VIN2、VIN3基因，推测文冠果春化机制与拟南芥存在差异^{[26, 27, 28]}。

赤霉素途径最终作用于GAI、RGA和RGL1这3个基因，从而影响开花。GAI在所有组织中广泛表达，RGL1仅在萌发的种子、花芽和果荚中高度表达^[29]。这些基因对花时起主要调节作用。DELLA蛋白在GA调控花时和花器官发育中是重要参与者^{[30, 31]}。本研究在文冠果转录组中发现了277条GAI同源序列，RGL1同源序列39条，DELLA1同源序列38条。

此外，植物内部还存在控制开花的自主途径。自主途径的一系列基因则对FLC的表达起负调控作用。在拟南芥中克隆到了FCA、FY、FLD、FPA、FVE、LD和FLK共7个基因^{[32, 33]}。本研究发掘出与自主途径相关基因，如FLD同源序列13条、FVE同源序列3条。以上分析表明，文冠果的成花调控途径比较复杂，其调控机制还有待于进一步研究。本研究样品的采集时期是成花转变的关键时期，结果将为进一步研究花发育分子机制提供基础和可利用的基因资源。