生殖调控如花期、花量与性别的调控，以及控制授粉与避免自交等操作是针叶树种遗传改良与良种生产中的重要环节。在长期的生产过程与实践中，育种工作者已经积累了很多有益的经验，如外施GA₄₊₇可以促进雌球花花量(Ross et al.，1984)，而遮荫可以促进雄球花形成而抑制雌球花形成(Ross et al.，1987)。但是对针叶树雌雄球花发育及调控的分子生物学研究却严重滞后，很重要的一个原因是它们的基因组明显比大多数其他植物基因组大(Morse et al.，2009)，全基因组测序工作进展缓慢，直至2013年5月才公布了第一种针叶树———欧洲云杉(Picea abies)的全基因组草图(Nystedt et al.，2013)。另外，与被子植物显著不同，祼子植物生殖器官多是单性而且缺少花被和心皮(Wang et al.，2010)，以致许多在被子植物生殖遗传学中取得的重要进展对祼子植物中的相关研究推动作用有限(Mathews et al.，2012)。

第2代测序技术的快速发展为这种困境的改变提供了契机(Neale et al.，2008)，转录组测序只关注基因组中转录的部分，可以避开针叶树基因组结构复杂、多重复序列等不利因素，可以较低的费用、快速地获得大量的基因组水平上的遗传信息(Fernandez-Pozo et al.，2011; Lesser et al.，2012)。随着测序成本的降低，公共数据库中基因组信息呈现爆炸式增长，这为基因芯片技术的发展和应用提供了基础，基因芯片已经成为生命科学研究中成熟而有效的工具。目前针对拟南芥(Arabidopsis thaliana)中所有预测基因，已经建成详细的基因时空表达数据库(Schmid et al.，2005)，其中数据的收集依赖于成熟的基因芯片技术与数据分析的有效应用。

转录组测序结合基因表达谱芯片分析为针叶树遗传研究提供了一种可行的、高效的研究手段。采用这种策略，本研究从油松(Pinus tabulaeformis)多个体多组织样本中获得大量的转录本信息，并以此为基础，设计并定制了油松高通量表达谱基因芯片。通过芯片杂交，比较油松雌雄球花发育后期的基因表达谱，筛选出一批稳定的在雌雄球花中表达差异极显著的基因，并对其中2个分别在雌雄球花中优势表达的基因进行Real-Time PCR分析。本研究为油松及近缘树种的雌雄球花发育以及其他生物学发育过程的分子生物学研究提供了基础和新的思路。

油松材料均取自平泉县七沟林场国家重点油松良种基地(41°02' N，118°68' E，海拔480～740 m)。油松转录组测序使用Roche 454 GS-FLX Titanium测序平台。基因芯片是基于获得的转录组信息设计、定制于Agilent公司的60 K油松表达谱芯片。

转录组测序样品: 2011年5月和7月分2次随机选取10株油松，采集针叶、茎干(包含正常木、应拉木及应压木的木质部、形成层、韧皮部组织)和球花(包含不同发育阶段的雌雄球花)混合样品。

芯片杂交样品: 2012年4月20日，随机选取3株成熟油松，分别采集雌雄球花样品。

基因表达分析样品: 将雌雄球花从表型可辨别到发育成熟的过程平均分为4个阶段，由于雌雄球花发育不同步，所以从2012年3月24日到4月20日每隔9天采集1次雄球花样品，从2012年4月2日到4月21日每隔7天采集1次雌球花样品; Realtime PCR分析采集的4个不同发育阶段的球花样品，取3个生物学重复，每个重复所有阶段样品均取自同一株树; 由于前期样品太小，为了提取RNA的需要，所以前2个发育阶段样品为采集于同一枝条上的3个混合花芽。

样品采集后立即放入液氮里带回实验室。在后续试验之前所有样品均保存于-80 ℃超低温冰箱。

油松组织样品在液氮中研磨成粉后用TRIzol(Invitrogen)试剂提取总RNA，RNA 用电泳和紫外分光光度法检验质量，总RNA用QIAGEN RNeasyMini Kit 进行纯化。以纯化后的RNA为模板，用SMART PCR cDNA synthesis Kit(Clonetech)合成cDNA; 用Agilent公司芯片杂交系统试剂合成cRNA，用QIAGEN RNeasy Mini Kit进行纯化。试验步骤参见各试剂盒操作说明书。

用GS FLX TitaniumGeneral Library Preparation Kit进行转录组测序文库构建; 用Roche 454高通量测序平台测序。转录组测序原始数据已经上传NCBI SRA(Sequence ReadArchive)公共数据库(登录号: SRA056887)。所分析的2个基因，MCE1在本数据库中的编号为lw_isotig04844，FCE1为lw_isotig08739。

测序原始数据经SeqClean和LUCY2软件去除接头和低质量部分(Li et al.，2004)，用Newbler2.5.3软件拼接(Kumar et al.，2010)，拼接好重叠群(isotig)和单一序列(singlet)用CD-HIT软件聚类分析(Huang et al.，2010)，相似度大于95%的序列聚为一类，其中最长的1条序列代表1个单基因(unigene)用于芯片设计。所有单基因用GeneOntology(GO)数据库注释(Garzon-Martinez et al.，2012)。

将序列信息提交到Agilen公司网站，利用该平台进行芯片设计与定制(https://earray.chem.agilent.com/earray)。提取油松雌雄球花总RNA并合成cRNA，每组样品取3个生物学重复。纯化后的cRNA用Agilent公司芯片杂交系统试剂片段化后加入杂交液，上芯片，65 ℃，10 r·min^-1滚动杂交17 h。杂交后的芯片经洗涤后用Agilent扫描仪进行扫描，分辨率为5 μm，以100%和10% PMT各扫描1次，2次结果用Agilent软件合并。

以油松4个不同发育阶段的雌雄球花为材料，提取总RNA并合成cDNA，用SuperReal PreMix(SYBR Green)(TIANGEN)进行Real-Time PCR分析。程序为: 95 ℃15 min预变性; 95 ℃ 10 s，60 ℃ 30 s，重复40个循环; 每个阶段取3个生物学重复，每个样品做3个技术重复。基因表达分析所用引物序列信息如下: 18SrRNA-F: 5'-cggctaccacatccaaggaa-3'，18S rRNA-R:5'-gctggaattaccgcggct-3'; MCE1-F: 5'-gctgcctctctacttcg-3'，MCE1-R: 5'-tgctgtcttctctcccac-3'; FCE1-F: 5'-catgaaattgttgagaaggca-3'，FCE1-R:5'-ggactatgacccacagga-3'。

基于Roche 454/GS FLX高通量测序平台，利用1个测序反应(a run)，共获得了822 891条EST序列，平均长度为358 bp，高质量的EST序列进一步拼接出46 584个单基因簇(unigenes)。

将获得的单基因簇进行GO功能注释，由于裸子植物基因功能研究不足以及与被子植物功能基因的差异，只有约1/2的基因可以预测出功能，其中21%参与生物学过程，12%属于细胞组分，17%属于分子功能基因(图 1)。

	图 1 油松转录组单基因簇GO注释分类 Fig. 1 The unigenes of P. tabulaeformis transcriptome functionally categorized into GO categories

利用Agilent软件对杂交芯片扫描图进行质量分析表明，芯片样本背景信号值与噪点值都在适当范围内，芯片数据质量较高，可以用于后续分析。对雌雄球花芯片数据进行主成分分析(principal component analysis，PCA)，结果显示2组样本内3个生物学重复间的重复性很好，而雌雄球花2组样本聚类差别非常明显，表明它们在各自的发育过程中基因表达差异极大(图 2A)。

	图 2 油松雌雄球花芯片数据聚类 Fig. 2 The cluster of the microarray data between female cones and male cones A: 芯片数据主成分分析( PCA) ; B: 芯片数据的GO 注释分类。 A: The principal component analysis (PCA) of the microarray data; B: The GO annotation categories of microarray data.

比较雌雄球花组织中基因表达差异，参照FC(fold change)≥5，3个生物学重复样品间P≤0.01的标准，共筛选到2 460个基因。将这些基因进行GO功能注释，有85%的基因有预测的功能，其中29%参与生物学过程，26%属于细胞组分，30%属于分子功能基因(图 2B)。

从这2 460个基因的散点分布图(图 3A)和频度统计图(图 3B)可以看出，在取样期，雄球花组织中优势表达的基因相对较多，而且雌球花中优势表达的基因更多地集中在相对差异较低区域。有1 309个基因在雄球花中优势表达，其中有20%的基因的表达量比雌球花中高64倍(log₂64=6)以上;有1 151个基因在雌球花中优势表达，但只有7%的基因的表达量比雄球花中高64倍以上。如果参照更加严格的重复样品间P≤0.001的标准，表达差异大于64倍以上的基因占总数的比例分别为26%(雄球花)和10%(雌球花)，共计322个。

	图 3 雌球花与雄球花差异表达基因的频度分布 Fig. 3 The distribution of differentially expressed genes of female cones vs. male cones 图A上方点集与图B负值统计区域代表在雄球花中优势表达的基因，图A下方点集与图B正值统计区域代表在雌球花中优势表达的基因。 The genes predominantly expressed in male cones fall into the upper left of A and the negative part of B, and the genes predominantly expressed in female cones fall into the lower right of A and the positive part of B.

从筛选到的2 460个基因中随机选取100个基因(图 4A)，另外选取雌雄球花中表达差异最大的各50个基因(图 4B)进行可视化绘图。从热图中可以看出所得到的雌雄球花中差异表达基因的数量和表达量的差异是稳定的。推测这些基因与雌雄球花发育的进程密切相关，同时为寻找关键调控基因提供了思路。

图 4 雌雄球花中差异表达基因的热图显示 Fig. 4 Heat map of differentially expressed genes in the male and female cones A: 随机选取的100个基因; B: 在雌雄球花中分别差异最大的50个基因。图中单元格颜色越深代表表达量越高。F1，F2，F3指雌球花3个生物学重复，M1，M2，M3指雄球花3个生物学重复。 A: The 100 genes randomly selected. B: The 50 genes most predominantly expressed in male cones and female cones, respectively. Colors indicate relative expression: from dark (highly abundant) to white (weakly abundant). F1 , F2, F3: The three biologicalreplicates of female cones; M1, M2, M3: The three biological replicates of male cones.

从雌雄球花差异表达基因中，选取1个在雄球花中优势表达基因，命名为MCE1(Male Cone Expressed 1)，芯片数据显示MCE1在雄球花中的表达量是雌球花中表达量的5 844倍(P=1.73×10^-6)，BLAST 结果显示MCE1在拟南芥、杨树及其他松树中存在同源基因，与现有数据库中白云杉(Picea glauca)的未知功能序列高度同源［NCBI登录号: DR555469; 相似性为95%(499/528)］; 选取1个在雌球花中优势表达基因，命名为FCE1(Female Cone Expressed 1)，芯片数据显示FCE1在雌球花中的表达量是雄球花中表达量的849倍(P = 5.78 × 10^-6)，通过与公共数据库中已知序列BLAST 结果显示FCE1与现有数据库中海岸松(Pinus pinaster)的未知功能序列高度同源［NCBI登录号: CT581928; 相似性为96%(633/657)］。

为了研究MCE1和FCE1在雌雄球花发育过程中可能的生物学功能，进一步分析了它们在雌雄球花不同发育阶段的表达规律。将雌雄球花的发育过程分为4个阶段，分别标记为F1-4和M1-4，提取每一阶段不同样本的总RNA，反转录成cDNA后进行Real-time PCR分析。

分析结果表明，MCE1在不同发育阶段的雄球花中表达规律呈单峰型，在M3阶段表达量最高，比M1阶段表达量高16倍(图 5A)。MCE1在雌花中也呈单峰型表达规律，在F2阶段表达量最高，不过MCE1在任何阶段的雌球花中始终都比雄球花中的表达量低1 000倍(log₂1 000=10)以上(图 5A)。上述结果与芯片分析结果相一致，同时说明维持MCE1高水平表达对雄球花的发育至关重要。

图 5 MCE1和FCE1在不同发育阶段雌雄球花中的表达 Fig. 5 The expression analysis of MCE1 and FCE1 gene of male and female cones in the different developmental stages M1-4与F1-4分别指按时间顺序不同发育阶段的雄球花与雌球花样本;图中表达差异倍数均以F1阶段中的样本1的数据为参照(设为1)。 The M1-4 and F1-4 indicate male cone and female cone samples from the early to the late developmental stages; The relative expression levelsare normalized against the sample one of F1 (set at 1).

FCE1在不同阶段的雌球花中一直维持在较高的表达水平，不过其中F1阶段表达量稍低，约为其他阶段表达量的1/4(图 5B)。有趣的是FCE1在雄球花M1阶段表达量比任一阶段雌球花中的表达量都高16倍以上，但在发育后期，其mRNA含量分别以29倍、221倍的速度急剧降到非常低的水平，到芯片分析阶段(F4与M4)，它在雄球花中的表达水平已经只有雌球花中的1 /1 000左右(图 5B)。这些结果表明FCE1并不只在雌球花的发育过程中起作用，而对雄球花早期的发育也十分重要。推测其在雌雄球花不同发育阶段可能发挥着不同的作用。

油松是我国华北、西北及东北部分地区特有的重要针叶用材树种，分布于中国北方的14个省市自治区，生态适生区达300万km²(Li et al.，2011)。本研究基于Roche 454/GS FLX高通量测序平台，对油松多个个体、多种组织样本转录组进行了均一化混合测序，获得了822 891条EST序列，极大地丰富了公共数据库中油松基因组遗传信息(此前只有288条EST序列信息)。为油松以及近缘树种比较基因组学、基因表达谱分析及分子标记开发等研究提供了重要资料。

基于获得的转录本信息，设计并定制了60 K油松表达谱芯片。通过芯片杂交，筛选到了2 460个在雌雄球花中差异表达的基因，其中有90个差异在500倍以上。这些基因的筛选为雌雄球花发育过程的分子生物学调控机制研究提供了有益的参考。

从中选取了2个分别在雌雄球花中优势表达的基因，由于功能未知，分别命名为MCE1与FCE1;基因表达分析表明MCE1可能在雄球花发育过程中起着重要作用，而FCE1可能在雌雄球花不同发育阶段起着不用作用，它的表达对雄花的早期发育可能具有调控作用。另外由于MCE1特异地在雄花各个阶段高效表达，它的启动子可以成为以后油松及近缘树种雄球花遗传操作中重要的工具。

本研究结果表明油松雌雄球花的发育是2个极其复杂的生物过程，涉及到众多的基因，这些基因还可能在不同发育阶段起着不同的作用。对油松雌雄球花发育调控分子生物学基础解析需要针对这些基因开展更加详细和深入的研究。关于油松的生殖调控研究可以为以针叶树为代表的祼子植物与被子植物生殖调控的比较生物学研究提供基础，有助于更加深入地理解花的起源与进化进程; 生殖调控分子生物学基础的解析可以为种子园中父母本的性别调控提供理论基础和新的调控策略，为提高种子产量与品质提供技术手段; 生殖调控在用材林中的潜在应用，可以实现大幅度地降低成熟林在生殖过程中的营养消耗，提高木材增长量。因此，关于油松的生殖调控的研究具有重要的理论与实践意义，需要开展更加详细和深入的研究。

目前生殖调控的分子生物学基础研究在被子植物中，特别是在拟南芥(Drouaud et al.，2007; Giraut et al.，2011)、黄瓜(Cucumis sativus)(Kater et al.，2001)、玉米(Zea mays)(Banks，2008; Dellaporta et al.，1994; Hartwig et al.，2011)等模式植物中，取得了较大的突破。但是在祼子植物中，相关研究却相当缺乏。

在第二代测序技术快速发展的背景下，本研究中采用了转录组测序和基因芯片相组合的方法，为祼子植物生殖调控的转录组水平分子生物学基础研究提供了一种有效的可操作方案。本研究的结果与方法为油松及近缘树种生殖调控的深入研究提供了基础。