山茶属(Camellia)油茶组(Sect. Oleifera)植物是我国重要的木本油料树种，茶油富含的不饱和脂肪酸对人体具有重要的保健功能(韩宁林等，2009; 王文杰等，2007)，其中油茶(Camellia oleifera)是我国目前第一位的主栽物种，不管其栽培面积还是产量都居主导位置(庄瑞林，2008)。油茶是典型的“抱子怀胎”植物，雌雄同花，花期一般在10月初到翌年2月，盛花期在11月中旬(查小华等，2012)。春梢生长量非常大，一般春梢占三梢(春梢、夏梢和秋梢)98%以上(唐光旭，1984)。油茶花开在当年生的春梢上，花量非常大，一支7张叶片的春梢上常常能开30朵以上，因此春梢是油茶的主要结果枝，与油茶产量密切相关(黎章矩等，1992； 曾燕如等，2009)。油茶是典型的常绿植物，春梢上的叶片一般有15～16个月的寿命，春梢一般3月中下旬开始萌发，5月初开始花芽分化。油茶自交可育性低，如果花期遇到寒冷和阴雨天气会影响授粉受精，同时幼果树上过冬，不同品种及立地条件下的油茶产量差异非常大(林少韩等，1981; 高本年，1981)。因此研究油茶花发育特点，并对花期及花器官形成进行有效调控，对解决油茶产量及提高油茶育种效率具有重要作用。目前对花发育调控方式主要有2种，即通过农事和转基因手段，转基因技术成功运用于油茶的花发育调控，其调控效果较农事手段明显。转基因成功的关键，除建立良好的遗传转化体系外，就是要清楚被调控生物的分子发育机制。

目前对油茶的成花机制研究鲜有报道，多数研究侧重于花芽形态和受精方面的研究(袁德义等，2011; 王湘南等，2011)。全基因组测序是研究生物发育分子机制的重要内容和手段，但是油茶的染色体数为2n=90，其基因组非常庞大，基于目前的测序条件，其成本相当昂贵。从目前相关文献来看，侧重了解生物发育某一方面的基因表达，最有效的办法是通过新一代测序技术对油茶花发育过程的转录组进行高通量测序(Mortazavi et al.，2008)。目前油茶转录组测序虽然也取得一些进展(林萍等，2011; 陈英等，2011)，但是侧重在花发育方面的研究却未见报道。本文对油茶花发育过程的6个阶段的转录本进行高通量测序，并探讨花发育相关基因FLC(Flowering locus C)，FCA(Flowering time control)，FT(Flower locus T)，AP1(Apetala1)，AP2(Apetala2)和PI(Pistillata)在花芽发育过程中的表达变化，以期为油茶成花机制研究提供基础数据。

油茶长林4号(Camellia oleifera‘ChangLin 4’)优良无性系。

根据相关文献报道(庄瑞林，2008; 袁德义等，2011; 王湘南等，2011)，分别在花芽分化前(5月8日)、萼片形成期(5月28日)、花瓣形成期(6月15日)、雌雄蕊形成期(7月3日)、子房和花药形成期(7月25日)及雌雄蕊成熟期(8月25日)采集花芽用于总RNA分离及转录组测序； 在春梢萌芽期(3月19日)分别采集新芽和根尖，而后每隔10天取1次茎尖，最后1次为5月28日(花芽分化完成进入萼片形成期)。样品采集后置于液氮中带回，放在-80 ℃的超低温冰箱中保存备用。

RNA提取采用艾德莱公司EASYspin Plus植物RNA快速提取试剂盒(目录号：RN38)，采用微量分光光度计Nanodrop 2000及琼脂糖凝胶电泳检测RNA质量，保证RNA样品浓度400 ng·μL^-1，28S:18S大于1.8。

利用富含ploy(T)低吸附磁珠对Total RNA中的mRNA进行富集纯化处理，在高温条件下，利用二价阳离子打断mRNA以选取合适大小的目的片段，利用反转录酶和随机引物将打断后mRNA片段反转录形成cDNA第1条链，然后加入缓冲液、dNTPs，在RNase H和DNA polymerase I的作用下合成第2条cDNA，对反转录合成的双链cDNA进行末端修复、3’末端加A、连接接头，通过琼脂糖凝胶电泳，筛选一定范围大小的片段。将前步处理好的样品进行精确定量(qubit)，在芯片(flow cell)表面进行桥式PCR，使DNA片段扩增为单分子DNA簇(此过程在Cluster Station中进行)，单分子DNA簇形成后将Flow cell 移入Hi-Seq中，进行Illumina HiSeq^TM 2000测序。

测序得到的原始数据去除只含有测序接头序列，或N含量过高以及过短的序列数据。利用Trinity(http://trinityrnaseq.sourceforge.net)软件通过序列之间的Overlap将序列延伸成Contig，再根据序列的Paired-end信息，将 Contig连接成转录本序列； 对于Trinity组装结果，通过组装的Component从潜在的可变剪接转录本选取组装最长的转录本作为该样品的Unigene序列，最后将Unigene序列与Nr，Nt，SwissProt，TrEMBL，GO，COG，KEGG数据库比对，获得Unigenes 的注释信息。如果不同库之间的比对结果有矛盾，则按Nr，Swiss-Prot，KEGG和COG的优先级确定Unigene的序列方向，比对不上的Unigenes 则用软件ESTScan预测其编码区并确定序列的方向。

根据相关文献及前期试验在转录组测序数据中，利用同源比对，筛选出FT，FCA，FLC，AP1，AP2和PI基因序列，利用Primer 5进行引物设计，引物序列见表 1。将1.3中提取的RNA，通过TAKARA公司PrimeScript^@Reverse Transcriptase试剂盒(商品编号：D2680A)方法(参考说明书)合成cDNA第1链，取其中4 μL合成产物稀释10倍，用于试验； 采用CESA内参(Zhou et al.，2013)，根据TAKARA公司SYBR^@Premix Ex Taq^TM(Tli RNaseh Plus)(商品编号为DRR420S)说明书进行荧光定量PCR反应(ABI 7300)，每个反应重复3次。

表 1 引物序列 Tab.1 Primer sequences

表 1 引物序列 Tab.1 Primer sequences 基因名称 Gene name 引物序列 Primer sequences 基因长度 Gene length/bp 退火温度 Annealing temperature(Tm)/℃ FT FTF：5′-GATGCTCCCAGCCCAAGTGAT-3′ 231 61.9 FTR：5′ -CTCCGCAAAGTCCCGAGTGTT-3′ FCA FCAF：5′-TCATCAGAGGACCAAGACGACG-3′ 138 61.9 FCAR：5′-TGGTGTGGGCAGATGGATAGAG-3′ FLC FLCF：5′-GCAAATCTGGACTCTGCTTCAC-3′ 205 60.1 FLCR：5′-GAGGTGGATTGGTCAAGAAAGC-3′ AP1 AP1F：5′-GGAGATGCGTGCTTGTTGAGGA-3′ 247 61.9 AP1R：5′-ACCCTAAGCCTCCCAACAGCAT-3′ AP2 AP2F：5′-GCTGCTTCACCACTCTTACCTCT-3′ 202 62.0 AP2R：5′-CACTGAGTCTGAGGTGGAAGGTA-3′ PI PIF：5′-TGGAAAGATGCTGGAGGAGGATA-3′ 255 60.2 PIR：5′-AGGCTTCTCTCCCACACAAAACA-3′

将花芽发育的6个时期的RNA按照等量的方式混合成1个样，测序总片段数(total reads)为28 448 847条，总碱基数为5.74 Gbp，GC含量为46.52%。CycleQ是测序质量评估的一个重要标准，其值越大测序误差越小，从图 1可以看出本次测序CycleQ20 percentage为100.00%，CycleQ30 percentage绝大部分也接近100.00%，说明测序结果符合要求，可以进行后续分析。

	图 1 CycleQ分布 Fig. 1 Distribution map of the CycleQ

使用 Trinity软件对样品数据进行组装，获得对应的Transcripts(转录本)序列。Trinity是一个专门对高通量RNA-seq数据进行全长转录本构建的软件，不同于其他的组装软件，最终构建出来的全长转录本不包含Gap，是一个比较好的无基因组全长转录本的组装软件。 N50长度是衡量组装效果的一个重要指标，从表 2中Contig和Transcripts的N50 长度的数值大小可以看出，组装效果基本符合要求。

表 2 测序数据组装和统计 Tab.2 Sequencing data assembling and statistics

表 2 测序数据组装和统计 Tab.2 Sequencing data assembling and statistics 重叠群Contigs 转录本Transcripts 重叠群长度 Contig length/bp 总数 Total number 占总数百分比 Percentage(%) 转录本长度 Transcript length/bp 总数 Total number 占总数百分比 Percentage(%) 0～300 5 150 211 98.71 200～300 48 993 30.96 300～500 35 277 0.68 300～500 37 781 23.87 500～1 000 18 361 0.35 500～1 000 30 411 19.22 1 000～2 000 9 602 0.18 1 000～2 000 27 357 17.29 2 000+ 4 180 0.08 2 000+ 13 721 8.67 总数Total number 5 217 631 158 263 总长度Total length 262 993 383 126 640 737 N50长度N50 length 49 1 381 平均长度Mean length 50.404 75 800.191 7

由于后续的表达谱分析是建立在同一个参考基因的基础上，故采用相似度的方式构建 Unigenes数据库，即当序列的比对长度达到较长序列的80%或者较短序列的90%时，认为2条序列同源； 取较长的序列作为代表，对样品先前聚类得到的基因序列作进一步的聚类，整合得到该项目物种Unigenes数据库。从表 3看出共得到 Unigenes 94 476 条，其中长度在 1 kbp以上的Unigenes 有12 643 条，占Unigene库总数的13.38%，N50长度为806 bp，说明组装效果较好。

表 3 Unigenes长度统计结果 Tab.3 The statistical results of unigenes length

表 3 Unigenes长度统计结果 Tab.3 The statistical results of unigenes length Unigenes 长度 Unigenes length/bp 总数 Total number 占百分比 Percentage(%) 200～300 40 923 43.32 300～500 27 155 28.74 500～1 000 13 755 14.56 1 000～2 000 8 386 8.88 2 000+ 4 257 4.51 总数Total number 94 476 总长度 Total length 53 644 108 N50长度 N50 length 806 平均长度 Mean length 567.806 723 400 652

从图 2可以看出，在目前的Nr数据库中，油茶转录组测序组装的Unigenes与葡萄(Vitis vinifera)的Unigenes相似数量最多，达到了53.44%，其次是其他物种，毛果杨(Populus trichocarpa)为第三，而与油茶近亲的种即茶树(Camellia sinensis)相似的数量仅为0.58%。

	图 2 在Nr数据库中Unigenes与其他物种的比较 Fig. 2 Unigenes similarity comparison with other species in Nr database

使用 BLAST软件将 Unigenes 序列与 Nr，Nt，SwissProt，TrEMBL，GO，COG，KEGG数据库比对，获得 Unigenes 的注释信息。从表 4可以看出，最终获得注释信息的Unigenes有41 240条，在COG和KEGG数据库中注释的Unigenes数量都在10 000条以下，其他数据库则都超过20 000条，大部分集中在30 000条左右。

表 4 Unigenes的功能注释在各数据库分布情况 Tab.4 The database distribution of unigenes functional annotation

表 4 Unigenes的功能注释在各数据库分布情况 Tab.4 The database distribution of unigenes functional annotation 注释数据库名称 Database name 被注释数量 Annotated number 300 bp≤长度< 1 000 bp 300 bp≤Length< 1 000 bp 长度≥ 1 000 bp Length≥ 1 000 bp COG 9 095 2 981 5 104 GO 27 201 11 406 10 515 KEGG 6 431 2 372 2 904 SwissProt 24 534 10 064 10 376 TrEMBL 36 393 16 403 12 113 Nr 36 400 16 412 12 117 Nt 30 858 13 165 11 460 被注释Unigenes总量 Annotated total number 41 240 18 743 12 220

Unigenes通过与KEGG数据库比对，对13 647条 Unigenes 进行了Pathway 注释，涉及的 Pathway有286个。其中有200条Unigenes以上的Pathway有：剪接体(346条)占总数的2.54%，分子伴侣和可折叠催化(247条)占总数1.81%，染色体(360条)占总数2.64%，内质网中的蛋白质加工(213条)占总数的1.56%，泛素系统(314条)占总数2.3%，植物激素信号转导(236条)占总数1.73％，转录因子(267条)占总数1.96%，糖基转移酶(224条)占总数1.64%，蛋白激酶(201条)占总数1.47%，DNA复制和重组蛋白(281条)占总数2.06%，氧化磷酸化(237条)占总数1.74%，核糖体(294条)占总数2.15%，其他Pathway的Unigenes数量都在200条以下，大部分集中在100条以内。涉及与植物花发育相关的Pathway有：植物激素信号转导，共有236条Unigenes，占总数1.73%； 植物昼夜调节，共13条Unigenes，占总数0.1%； 淀粉和蔗糖代谢，共有168条Unigenes，占总数1.23％。

在GO数据库中(图 3)，执行细胞组成功能(cellular component)Unigenes有42 106条，其中达到1 000条以上Unigenes，分别为：核5 102条、质膜4 259条、线粒体3 414条、叶绿体2 589条、不可或缺的膜2 116条、胞质溶胶1 999条、细胞膜1 548条、胞浆膜界囊泡1 413条、细胞质1 330条、质体1 194条、胞间连丝1 140条、胞外区1 053条、高尔基体1 047条，其他都在700条以下。执行分子功能(molecular function)有47 192条Unigenes，其中大于1 000条有：ATP结合(3 069条)、锌离子结合(1 680条)、核酸结合(1 658条)、核苷酸结合(1 470条)、金属离子结合(1 442条)、DNA结合(1 353条)、蛋白结合(1 229条)、丝氨酸/苏氨酸激酶活性(1 037条)、绑定(1 019条)。涉及油茶生物过程功能(biological process)的Unigenes有80 231条，其中最多是氧化还原过程1 833条，其次是蛋白磷酸化1 161条，其他全部都在1 000条以下。

图 3 Unigenes在GO中功能注释统计结果 Fig. 3 Unigenes functional annotation in the GO 细胞组成：A：细胞； B：细胞连接； C：细胞要素； D：细胞外基质； E：胞外区； F：胞外区要素； G：大分子复合物； H：膜； I：膜要素； J：膜封闭腔； K：核仁； L：细胞器； M：细胞器要素； N：共质体。分子功能：A：抗氧化活性； B：绑定； C：催化活性； D：电子载体活性； E：酶调节活性； F：分子转导活性； G：核酸结合的转录因子的活性； H：营养库活性； I：蛋白结合转录因子活性； J：受体活性； K：结构分子活性； L：转运活性。生物过程：A：生物附着； B：生物调控； C：细胞增殖； D：细胞成分和生物合成； E：细胞过程； F：凋亡； G：发育过程； H：定位系统建立； I：生长； J：免疫系统过程； K：定位； L：代谢过程； M：多机体过程； N：多细胞组织过程； O：色素沉着； P：再生； Q：再生过程； R：刺激应答； S：律动过程； T：信号传导。 Cellular component: A: Cell; B: Cell junction; C: Cell part; D: Extracellular matrix; E: Extracellular region; F: Extracellular region part; G: Macromoleular complex; H: Membrane; I: Membrane part; J: Membrane-enclosed lumen; K: Nucleoid; L: Organelle; M: Organelle part; N: Symplast. Molecular function: A: Antioxidant activity; B: Binding; C: Catalytic activity; D: Electron carrier activity; E: Enzyme regulator activity; F: Molecular transducer activity; G: Nucleic acid binding transcription factor activity; H: Nutrient reservoir activity; I: Protein binding transcription factor activity; J: Receptor activity; K: Structural molecule activity; L: Transporter activity. Biological process: A: Biological adhesion; B: Biological regulation; C: Cell proliferation; D: Cellular component organization or biogenesis; E: Cellular process; F: Death; G: Developmental process; H: Establishment of localization; I: Growth; J: Immune system process; K: Localization; L: Metabolic process; M: Multi-organism process; N: Multicellular organism process; O: Pigmentation; P: Reproduction; Q: Reproductive process; R: Response to stimulus; S: Rhythmic process; T: Signaling.

图 4可以看出，在COG数据库中Unigenes执行的功能分成：翻译，核糖体结构和生物合成，RNA加工和修饰，转录、复制、重组和修复，染色质结构和动力学，细胞周期调控，细胞分裂，染色体分区，核结构，防御机制，信号转导机制，细胞壁、膜和包膜合成，细胞运动，骨架，胞外结构，细胞内运输，分泌和膜泡运输，翻译后修饰，蛋白质，伴侣，能源生产和转换，碳水化合物的运输和代谢，氨基酸转运和代谢，核苷酸转运和代谢，辅酶运输和代谢，脂质运输和代谢，无机离子转运和代谢，次生代谢产物的生物合成、运输和代谢等。其中只有一般功能的Unigenes数量最多，共有2 489条，其次是具复制、重组和修复功能的Unigenes，数量共有1 420条，转录功能Unigenes数量有1 214条，信号转导机制Unigenes数量有1 065条，其他功能的Unigenes都在1 000条以下，其中较多的为翻译后的修饰功能，蛋白质周转和伴侣功能Unigenes有837条，糖类转运与代谢Unigenes有690条。

图 4 COG数据库中功能注释的unigenes分类 Fig. 4 COG function classification of unigenes A: 翻译、核糖体结构和生物发生; B: RNA加工和修饰; C: 转录; D: 复制、重组和修饰; E: 染色体结构和动力学; F: 细胞周期调控、细胞分裂、 染色体分离; G: 核结构; H: 防御机制; I: 信号传递机制; J: 细胞壁/细胞膜生物发生; K: 细胞活性; L: 细胞骨架; M: 胞外结构; N: 细胞间运输、分泌物和囊泡运动; O: 翻译后修饰、蛋白翻转、伴侣; P: 能源产生与转化; Q: 碳水化合物转运和代谢; R: 氨基酸转运和代谢; S: 核酸转运和代谢; T: 辅酶转运和代谢; U: 脂类转运和代谢; V: 无机离子转运和代谢; W: 次生代谢物生物合成、转运和代谢; X: 只有一般功能预测; Y: 未知功能。 A: Translation, ribosomal structure and biogenesis; B: RNA processing and modification; C: Transcription; D: Replication, recombination and repair; E: Chromatin structure and dynamics; F: Cell cycle control, cell devision, chromosome partitioning; G: Nuclear structure; H: Defense mechanisms; I: Signal transduction mechanisms; J: Cell wall/ membrane/envelope biogenesis; K: Cell motility; L: Cytoskeleton; M: Extracellular structures; N: Intracellular trafficking, secretion, and vesicular transport; O: Posttranslational modification, protein turnover, chaperones; P: Energy production and conversion; Q: Carbohydrate transport and metabolism; R: Amino acid transport and metabolism; S: Nucleotide transport and metabolism; T: Coenzyme transport and metabolism; U: Lipid transport and metabolism; V: Inorganic ion transport and metabolism; W: Secondary metabolites biosynthesis, transport and catabolism; X: General function prediction only; Y: Function unknown.

从图 5可以看出FLC，FT和FAC较AP1，AP2和PI在茎尖表达量相对较低。FLC在根系中表达量最高，5月18日表达量最低； FT在3月19日表达量最高，4月18日开始到5月18日表达量变动不大； FAC在3月份的表达量较其他时期高，其与FLC基因表达量存在相互消长的关系。PI在5月中旬表达量最高，AP2则在根部和5月中上旬表达量较高，AP1在3月份表达量都较高。

	图 5 成花相关基因在茎尖的表达 Fig. 5 The expression of flowering gene in shoot tips 除03-19含根尖部位，其他时期只指春梢茎尖部位。 In addition to 03-19 containing root tip, other times only refer to the spring shoot tip

从图 6可以看出，在花器官形成过程中AP2和PI基因相对表达水平较其他基因高，AP2基因在雌雄蕊成熟期(8月25日)相对表达量最高，其次是在雌雄蕊形成期(7月3日)，在萼片形成期(5月28日)相对表达量最低。PI基因在雌雄蕊形成期(7月3日)相对表达量最高，其次为子房与花药形成期(7月25日)，再次为花瓣形成期(6月15日)。AP1基因在花发育的5个时期表达差异不明显，FT和FAC在6月15日和7月3日表达量较其他时期高，但从总体上看在花器官发育过程中FCA和FLC基因的表达量较在成花转变期低。

	图 6 花器官形成过程中成花相关基因的表达 Fig. 6 The expression of floral organ gene during the process of floral organ formation

转录组学(transcriptomics)是功能基因组学研究的一个重要内容，它是从整体水平上研究细胞中基因转录的情况及其转录调控规律。转录组的概念由Velculescu等(1997)首先提出。转录组是指特定细胞在某一功能状态下全部表达的基因总和，一般指的是mRNA总和，基因表达具有特定的空间性和时间性特征(Lander et al.，2001)，因此知道mRNA序列是转录组研究的重要基础； 随着新一代测序技术的发展，其高通量性、低成本、快速、准确的特点可大大提高测序的效率，降低试验成本，为从RNA水平研究基因表达提供有效方法(Nagalakshmi et al.，2008; Rosenkranz et al.，2008; Hoen et al.，2008; Wilhelm et al.，2009; Jack，2004)。本研究从花发育的6个时期获得了近6 Gbp的测序数据，测序数据的拼接基本符合后续生物信息学分析，并获得了近12 643条1 kbp以上Unigenes，其中大部分具有完整的编码区，为后续准确获得目的基因提供便利，同时为有效运用转基因技术奠定基础。

根据现有的研究文献表明：CO(Constant)(Cerdan et al.，2003)是植物光周期途径的关键基因； LFY(LEAFY)(Weigel et al.，1992)和FT(Huang et al.，2005)是成花调控网络的关键基因。本文研究根据转录组测序的基础数据库，通过对CO和LFY基因的同源和相似性比对，结果只发现成花相关基因FT。Zeevaart(2008)研究发现，在拟南芥(Arabidopsis thaliana)的花发育不同阶段FT基因都有一定的表达，其主要在植物叶片中表达，在其他生物体部分(根、茎和花芽)表达量相对较低，这与本文分析结果一致。在花温度诱导自主成花途径中(Blazquez，2001)，发现大量编码热激蛋白的基因HSP(heat shock proteins)、促花基因FCA及抑花基因FLC。本研究通过对测序数据中HSP基因的表达量进行分析发现，随着环境温度升高，油茶花芽中HSP基因表达量急剧增加，表达量极高，通常是其他基因表达量的10³倍以上，同时对油茶采样株系的观察可以发现，随着春季环境气温回升，油茶的春梢生长在40天里可以生长50 cm。从本文研究中还可以看出FCA基因在春梢发育初期及花芽分化期表达量明显增加，在花芽分化完成后的花器官形成期表达量相对较小，FLC基因其表达量与FCA存在相互消长的关系。AP1和AP2都属于MADS-box基因，是真核生物中一类重要的转录调控因子，属于开花ABC模型的A类基因(Mandel et al.，1992; Bowman et al.，1993)； 从本文成花过程表达分析来看，A类基因不仅与花的萼片形成有重要关系，而在前期的花芽转变过程也可能具有调控作用。PI属于B类基因，对花的雄蕊和花瓣形成具有重要作用(Howell，1998)； 从本文花芽分化过程各个时期的表达量来看，PI基因还可能影响油茶的花芽分化。因此对FT及PI基因的改造及遗传改良可以获得早花和雄性不育油茶品种，为油茶高效栽培和育种打下良好基础。

目前植物成花转变主要有光周期途径、春化途径、赤霉素途径、糖代谢途径和自主途径(潘瑞炽，2004; 王忠，2002)。在本研究中没有发现与植物光周期关键基因CO的同源基因，也没有发现与春化相关的基因及控制花发端由营养生长转变成生殖生长的关键基因LFY，而发现随着环境温度的升高(从18 ℃到37 ℃)合成热激蛋白的相关基因(HSP)表达量急剧攀升，春梢在良好水肥条件下快速生长，并在春梢生长停滞期完成了花芽分化，开始进入花器官形成阶段； 鉴于前期试验发现，长林4号油茶经赤霉素处理对其花芽分化存在一定的抑制作用，而在室温(25 ℃)进行不同光周期和光强照射则不能完成花芽分化。结合上面分析可以初步判定，油茶成花途径是基于一定温度作用(持续高温)自身遗传决定的自主过程。首先，油茶春梢在持续升高的环境温度作用下编码热激蛋白的HSP基因表达量迅速升高，油茶在良好的水热及养分条件作用下春梢快速生长，春梢生长过程中促进成花基因FCA表达，并抑制FLC基因表达，同时上调成花网络关键基因FT，FT基因表达促进AP1和PI基因表达，而AP1和PI基因的表达引起花器官一系列相关基因的表达，从而完成了花芽分化。