马尾松(Pinus massoniana Lamb.)是松科(Pinaceae)松属的双维管束松亚属(subgen. Pinus)多年生常绿树种，它是我国松属树种中分布最为广泛的树种，也是南方荒山绿化、工业用材、纸浆造纸的重要树种，在创造经济价值和发挥生态效应方面占据着重要地位。目前已有马尾松遗传多样性研究、遗传连锁图谱构建等方面的相关报道(张薇等, 2008; 蔡娟娟和季孔庶, 2009)，但马尾松基因组及转录组数据的缺乏，造成马尾松生长发育相关研究、分子标记开发和遗传图谱构建等研究相对滞后。

EST技术是开展转录组研究的有效方法(李滢等, 2010)，已被广泛应用于动、植物和微生物的基因表达谱分析、功能基因预测等(Adams et al., 1991; Boguski et al., 1994; Bouchez and Höfte, 1998; Ewing et al., 1999)。目前公共数据库中已有松属多个树种的EST序列，如火炬松、北美短叶松、海岸松(Allona et al., 1998; Johnson et al., 2002; Temesgen et al., 2001; Chagné et al., 2004)。截止2011年4月，从GenBank查找发现松属EST序列有454 455条，但未发现马尾松的序列。随着测序技术的发展，将高通量测序应用到马尾松的转录组研究中，不仅可降低测序的成本和时间，而且可获得的丰富的数据，有利于松属的转录组和生长发育相关研究。

本研究首次采用Illumina进行松属树种转录组的研究，构建马尾松均一化cDNA文库、测序，对序列再做基因功能注释和分类等，以研究马尾松生长发育过程中重要基因的表达。为进一步开展功能基因组研究提供帮助，同时也可为马尾松转录组SSR位点查找和信息分析，开发分子标记提供丰富的数据。

1结果与分析
1.1高通量测序和EST序列拼接
采用Illumina高通量测序技术对马尾松新梢的转录组进行测序，200个循环得到4~6 G的数据。经过序列拼接，获得83 680个contig，其中300 bp以上的contig数目是53 677，大于1 kb的有17 386。

1.2序列功能注释
通过BLAST序列比对，共有33 772个contig获得基因注释(表1)。33 772个基因被注释为相应的331 669对生物学功能(contig-GO function pairs)；10 647个基因被注释为1 029种酶功能。另外39 261个contig未被注释，认为是可能的新基因。

表1 GO数据库注释结果基本统计 Table 1 Statistics on the GO database annotation results

1.3功能分类研究
通过与蛋白质组序列比对，获得马尾松contig的GO分类信息。其中分别归入分子功能(MF)的有26 686个contig，生物学过程(BP)的有29 762个contig，细胞组分(CC)的有30 250个contig，并对之进行具体归类，该结果说明了马尾松新梢基因表达谱的总体情况(图1; 图2; 图3)。其中，分子功能中具有代表性的功能类别为与蛋白结合、ATP结合、序列特异性等相关；生物学过程中占优势的为与生长相关、盐胁迫响应、抗病防御响应如抗细菌和真菌、叶绿体RNA加工等过程；细胞组分的代表类型为细胞质、细胞膜、叶绿体、线粒体、细胞核等。

图1 Contig GO注释的层次(类别)分布图-MF Figure 1 Distribution bar chats of the levels for contig GO annotation-MF

图2 Contig GO注释的层次(类别)分布图-BP Figure 2 Distribution bar chats of the levels for contig GO annotation-BP

图3 Contig GO注释的层次(类别)分布图-CC Figure 3 Distribution bar chats of the levels for contig GO annotation-CC

1.4生物学通路分析
生物学通路主要包括代谢通路与信号转导通路。根据KEGG pathway数据库，对马尾松转录组的contig进行Pathway生物学通路的注释和预测，共识别出10 647个contig具有对应的1 029种酶(Enzyme, EC)功能，并关联到135条生物学通路(pathway)。

马尾松属针叶树种，在新梢的生长过程中会有大量的角质、软木质和蜡质的合成，以保护其体内水分的散失，维持正常的生长发育过程。通过同源性比较，发现47个contig可能参与编码马尾松角质、软木质和蜡质合成途径的3个关建酶(表2; 图4; 图5)，包括长链乙醇-O-脂肪酰基转移酶、二酯甘油-O-酰基转移酶和长链脂肪-酰基辅酶A还原酶。其中长链乙醇-O-脂肪酰基转移酶含有11个contig，二酯甘油-O-酰基转移酶含有14个contig，长链脂肪-酰基辅酶A还原酶含有33个contig。

表2 马尾松角质, 软木质和蜡质的生物学通路注释结果 Table 2 The annotation results of biological pathways for cutin, suberin and wax in P. massoniana

图4 马尾松角质, 软木质合成的生物学通路图 Figure 4 Biological pathway of the cutin and suberin biosynthsis in P. massoniana

图5 马尾松蜡质合成的生物学通路图 Figure 5 Biological pathway of the wax biosynthsis in P. massoniana

此外，细胞色素P450家族在植物防御反应中的作用、对除草剂等外源化学物质的解毒功能以及在一些物质的生物合成(如木质素中间物, 植物激素, 甾醇, 萜类, 黄酮类, 异黄酮和呋喃香豆素等)中起重要作用(戴素明等, 2004; 李明等, 2006; 余小林等, 2004)。从马尾松EST文库注释结果中找到165个contig可能参与编码细胞色素P450基因，其参与编码多种关键酶，包括羰基还原酶(NADPH)、乙醇脱氢酶、11-β-羟基脱氢酶、醛脱氢酶([NADP+])葡萄糖醛酸转移酶、谷胱甘肽转移酶、非特异性单加氧酶等，同时参与细胞色素P450的外源物代谢之中。

1.5 SSR信息分析
从马尾松83 680个contig中查找到889个SSR位点，占contig总数的比例为1.06%。SSR的类型丰富，二至六核苷酸重复类型均存在，所占比例变化较大(表3)。其中，三核苷酸重复所占比例最高，达到48.37%，其次是六核苷酸重复，为19.12%，比例最低的是四核苷酸重复，仅为4.72%，二核苷酸重复和五核苷酸重复基本相同，分别为14.62%和13.16%。在检出的SSR中，共发现349种基元类型，出现频率最高的前10种重复基元为AT/AT (82个)，AGC/CTG (78个)，AAG/CTT (66个)，AAT/ATT (64个)，AGG/CCT (63个)，ATC/ATG (56个)，AAC/GTT (43个)，ACC/GGT (31个)，AG/CT (30个)，AC/GT (18个)。上述SSR的特征分析，有助于开展马尾松及其同属物种的基因组差异分析、通用性标记开发和遗传图谱构建等的研究。

表3 马尾松EST-SSR不同重复基元分布及优势碱基组成 Table 3 Distrbution and compositions of the dominant repeat for the different EST-SSR in P. massoniana

2讨论
Illumina高通量测序的数据量大，效率高，试验成本低，适合于当前马尾松等没有开展全基因组测序的物种开展转录组测序研究。本研究构建了第一个马尾松均一化cDNA文库，首次采用了Illumina高通量测序技术对文库进行了测序，序列拼接后得到83 680个contig，并采用生物信息学分析方法，对参与马尾松相关生长发育过程的基因进行注释，其中33 772个contig得到功能注释，获得39 261个新的contig，所占比例为46.9%，这一结果显示采用高通量测序可以大量挖掘马尾松生长发育过程中表达的重要基因。

本研究构建的马尾松均一化cDNA文库，可有效地使不同丰度基因的拷贝数趋于一致，克服了基因转录水平上的差异给测序效率带来的障碍；也相对提高了低丰度基因的拷贝数，使更容易检测出因时空特异性表达的低拷贝基因，大大增加了低丰度基因的测序几率，避免了相同基因的重复测序。

根据KEGG pathway数据库，对所得马尾松转录组的contig进行Pathway生物学通路注释和预测，共识别出10 647个contig具有对应的1 029种酶功能，他们参与到马尾松体内生化合成和次生代谢过程中，这为进一步开展马尾松基因克隆、功能基因的验证等分子试验提供帮助。

与传统方法相比，高通量测序是大量发掘SSR资源的简便方法，效率高，能达到批量开发SSR标记目标。本次试验通过SSR位点查找发现了889个SSR位点，且长度均在18bp以上，利用primer3.0在线引物设计软件共设计出401对SSR引物。下一步可对这些引物进行扩增检测，筛选出目的条带清晰、多态性好的引物，为马尾松群体遗传多样性分析和构建遗传连锁图谱等后续研究和分子育种奠定基础。

3材料与方法
3.1 材料
马尾松材料(Pinus massoniana Lamb.)采自南京林业大学南大山苗圃，剪下7年生马尾松的幼嫩新梢，迅速放入10 mL离心管内，立即放入液氮中，带回实验室保存于超低温冰箱里备用。

3.2 RNA提取和反转录
采用改进的CTAB-LiCl法提取马尾松新梢总RNA，用OligodT磁珠分离纯化mRNA。以mRNA为模板，TruSeq^TM RNA Sample Preparation kit(illumina)反转录合成dscDNA。末端补平，加A，加接头，进行PCR扩增，扩增条件为98℃，30 s；98℃ 10 s；60℃ 30 s；72℃ 30 s，15个循环。然后进行PCR产物的纯化。

3.3均一化cDNA文库构建和测序
参照DSN Normalization kit (illumine)的方法，利用DSN (duplex-specific nuclease)处理纯化后的dscDNA，取稀释的dscDNA进行不同循环数的PCR扩增，扩增条件为：98℃ 30 s；98℃ 10 s；65℃ 30 s；72℃ 30 s，12个循环，然后72℃，5 min，进行检测。应用新一代高通量测序平台Illumina对均一化cDNA样品上机测序。簇生成采用cBot Cluster Generation kit，测序模式为2×100。

3.3.1序列拼接
将测序片段进行k-mer长度的重叠拆分，并进行数据校正(包括通过对具有相同前缀但末端核苷酸不同的k-mer进行校正排除)，当处理全部测序片段成为k-mer并且过滤掉可能有错误的k-mer后，将丰度最大的，出现频率最高的k-mers作为组装拼接转录组草图的seed k-mers。Seed选择算法采用RepeatScout算法(Price et al., 2005)。采用5’端至3’端方向或3’端至5’端的覆盖度指导原则(coverage-guided manner)延伸组装转录组序列。在对样本组装拼接和合并组装拼接中，参数k都选择k=25。

3.3.2序列注释、功能分类和生物学通路分析
将测序得到的contig分别与已知的核酸、蛋白序列数据库进行BLAST比对，先进行功能注释(E值<0.001)，其中包括SwissProt、GenPept、PIR、PDF、PDB和NCBI RefSeq等。再将注释上的序列按照不同的功能划分成不同的类别。最后对KEGG注释的基因功能信息进行生物学通路的注释和预测。

3.3.3 SSR位点查找和分析
从contig序列中查找SSR位点，查找的标准为：二核苷酸重复不少于9次、三核苷酸不少于6次、四核苷酸不少于5次、五核苷酸和六核苷酸不少于4次。同时查找被小于等于100bp打断的复合型SSR。对查找的SSR类型进行特征分析。

作者贡献
季孔庶、王晓锋是本研究的实验设计和实验研究的执行人；王晓锋完成数据分析，论文初稿的写作；何卫龙、蔡卫佳、阮倩倩参与实验；季孔庶是项目的构思者及负责人，指导实验设计，数据分析，论文写作与修改。全体作者都阅读并同意最终的文本。

致谢
本研究由林业公益性行业科研专项经费(NO.201104010)、国家“十二五”科技支撑项目(NO.2012BAD01B02)和江苏高校优势学科建设工程项目资助。

参考文献
Adams M.D., Kelley J.M., Gocayne J.D., Dubnick M., Polymeropoulos M.H., Xiao H., Merril C.R., Wu A., Olde B., Moreno R.F., Kerlavage A.R., McCombie W. R., and Venter J.C., 1991, Complementary DNA sequencing: expressed sequence tags and human genome project, Science, 252(5013): 1651-1656

Allona I., Quinn M., Shoop E., Swope K., St Cyr S., Carlis J., Riedl J., Retzel E., Campbell M.M., Sederoff R., and Whetten R.W., 1998, Analysis of xylem formation in pine by cDNA sequencing, Proc. Natl. Acad. Sci., USA, 95(16): 9693-9698

Boguski M.S., Tolstoshev C.M., and Bassett D.E., 1994, Gene discovery in dbEST, Science, 265(5181): 1993-1994

Bouchez D., and Höfte H., 1998, Functional genomics in plants, Plant Physiol., 118(3): 725-732

Cai J.J., and Ji K.S., 2009, Diversity of SSR markers and strategy of genetic linkage map making with masson's pine Megagametophytes, Fenzi Zhiwu Yuzhong (Molecular Plant Breeding), 7(5): 934-940 (蔡娟娟, 季孔庶, 2009, 马尾松大配子体的SSR多样性及其遗传作图策略研究, 分子植物育种, 7(5): 934-940)

Chagné D., Chaumeil P., Ramboer A., Collada C., Guevara A., Cervera M.T., Vendramin G.G., Garcia V., Frigerio J.M., Echt C., Richardson T., and Plomion C., 2004, Cross-species transferability and mapping of genomic and cDNA SSRs in pines, Theor. Appl. Genet., 109(6): 1204-1214

Dai S.M., Zhou C.A., Xie B.Y., Feng D.X., and Xiao Q.M., 2004, Role of cytochromes P450 expression in plant defence responding to pathogens, Shihezi Daxue Xuebao (Journal of Shihezi University (Natural Science)), 22(S): 184-187 (戴素明, 周程爱, 谢丙炎, 冯东昕, 肖启明, 2004, 细胞色素P450表达在植物防御反应中的作用, 石河子大学学报(自然科学版), 22(S): 184-187)

Ewing R.M., Ben Kahla A., Poirot O., Lopez F., Audic S., and Claverie J.M., 1999, Large-scale statistical analyses of rice ESTs reveal correlated patterns of gene expression, Genome Res., 9(10): 950-959

Johnson A., Kinlaw C., Loopstra C., MacKay J., Neale D., O'Malley D., Peter G., Retzel E., VanZyl L., and Sederoff R., 2002, A genomic approach to wood formation in loblolly pine, In: Scherago International Inc (ed.), Plant, Animal and Microbe Genome X Conference, Town and Country Convention Center, San Diego CA, pp.12-16

Li M., Zeng R.S., and Luo S.M., 2006, Cytochrome P450s and plant chemical defense, Zhongguo Zhibao Daokan (China Plant Protection), 26(11): 13-15 (李明, 曾任森, 骆世明, 2006, 细胞色素P450酶系与植物的化学防御, 中国植保导刊, 26(11): 13-15)

Li Y., Sun C., Luo H.M., Li X.W., Niu Y.Y., and Chen S.L., 2010, Transcriptome characterization for Salvia miltiorrhiza using 454 GS FLX, Yaoxue Xuebao (Acta Pharmaceutica Sinica), 45(4): 524-529 (李滢, 孙超, 罗红梅, 李西文, 牛云云, 陈士林, 2010, 基于高通量测序454 GS FLX的丹参转录组学研究, 药学学报, 45(4): 524-529)

Price A.L., Jones N.C., and Pevzner P.A., 2005, De novo identification of repeat families in large genomes, Bioinformatics, 21(S): 351-358

Temesgen B., Brown G.R., Harry D.E., Kinlaw C.S., Sewell M.M., and Neale D.B., 2001, Genetic mapping of expressed sequence tag polymorphism (ESTP) markers in loblolly pine (Pinus taeda L.), Theor. Appl. Genet., 102(5): 664- 675

Yu X.L., Cao J.S., Cui H.M., and Ye W.Z., 2004, Plant cytochrome P450, Xibao Shengwuxue Zazhi (Chinese Journal of Cell Biology), 26: 561-566 (余小林, 曹家树, 崔辉梅, 叶纨芝, 2004, 植物细胞色素P450, 细胞生物学杂志, 26: 561-566)

Zhang W., Gong J., and Ji K.S., 2008, Genetic diversity for seeding orchard of masson's pine, Fenzi Zhiwu Yuzhong (Molecular Plant Breeding), 6(4): 717-723 (张薇, 龚佳, 季孔庶, 2008, 马尾松实生种子园遗传多样性分析, 分子植物育种, 6(4): 717-723)