楸树(Catalpa bungei)是紫葳科(Bignoniaceae)梓属(Catalpa)落叶乔木，是我国所特有的珍贵优质用材树种和观赏树种，素有“木王”之称。其树干通直，花形若钟，花朵盛开时极其优美，自古人们就把楸树作为观赏树种广泛栽培(王新建等，2004)。自然条件下，楸树个体存在雄性不育的现象，小孢子在四分体时期由于绒毡层结构和功能异常而不能正常发育，造成花药败育而不能形成正常的花粉(张博等，2015)。楸树通常采用嫁接和扦插等方式进行无性繁殖，长期的无性繁殖导致品种、类型和无性系单一化，且由于其自花不孕，常常使得楸树只开花不结实，加之对其过度的开发利用，造成楸树的资源缺乏。

近年来，转录组测序技术迅速发展，已被广泛应用到了SSR标记规模化开发(文亚峰等，2015；张振等，2015)、代谢途径的分析(陈昊等，2015；马婧等，2016)、基因挖掘等方面，同时也在植物雄性不育研究上应用并取得进展。根据已报道的雄性不育转录组数据分析表明，葱(Allium fistulosum)细胞质雄性不育与线粒体的氧化磷酸化有关(Liu et al., 2016a)，芝麻(Sesamum indicum)雄性不育的形成是在乙烯和茉莉酸介导的信号通路以及NAC、WRKY 2类转录因子共同作用下形成的(Liu et al., 2016b)，而卷心菜(Brassica oleracea var. capitata)细胞核雄性不育则主要涉及ATP合成酶(Guo et al., 2016)。西瓜(Citrullus lanatus)(Rhee et al., 2015)、辣椒(Capsicum annuum)(Liu et al., 2013)等植物也有报道，但在木本植物中通过转录组数据分析对雄性不育相关基因进行研究尚未见报道。

本研究通过转录组测序技术，对楸树雄性可育与不育花芽进行差异表达基因的分析，在基因转录水平上对楸树雄性不育的机制进行探索，分析楸树雄性不育发育过程中的特定基因的转录表达信息，以期为楸树雄性不育机制的研究提供理论参考，同时也为重要性状相关基因的克隆以及功能分析等奠定基础。

楸树花芽采自河南洛阳虎头山基地的楸树植株，以自然突变的楸树雄性不育花芽为试验组，记为SL，以雄性可育花芽作为对照组，记为FL。试验材料采集完成后立即用锡箔纸包好后做好标记用液氮处理，保存于-80 ℃冰箱中备用，转录组测序样本每一组各取3个生物学重复，分别为SL-1，SL-2，SL-3和FL-1，FL-2，FL-3。

采用多糖多酚植物总RNA提取试剂盒(北京天根生化科技有限公司)提取楸树雄性不育花芽和雄性可育花芽总RNA。利用NanoDrop检测RNA浓度，用Agilent 2100检测28S/18S以及RIN值，同时使用1%琼脂糖凝胶电泳检测提取RNA的质量以及完整性。

提取样品总RNA并使用DNaseⅠ消化DNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA；加入打断试剂在Thermomixer中适温将mRNA打断成短片段，以打断后的mRNA为模板合成一链cDNA，再通过二链合成反应体系合成二链cDNA，并使用试剂盒纯化回收、粘性末端修复、cDNA的3′末端加上碱基“A”并连接接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System质检，合格后使用Illumina HiSeq 4000或其他平台进行测序。

测序所得数据称为raw reads。首先，过滤掉低质量、接头污染以及未知碱基N含量过高的reads，过滤后的数据称为clean reads。对过滤后的数据使用Trinity进行De novo组装得到Unigene，接下来使用Tgicl对转录本进行聚类去冗余得到Unigene。

使用Blast对Unigene进行NT、NR、COG、KEGG以及SwissProt注释，使用Blast2GO以及NR注释结果进行GO注释，使用InterProScan 5进行InterPro注释。在NT数据库中比对分析非冗余核酸序列，利用NR、SwissProt数据库分析蛋白序列，然后利用COG数据库根据直系同源家族蛋白进行功能分类以及通过InterPro数据库根据蛋白结构域进行家族分类以提供功能注释。对所有Unigene基因在分子功能、细胞组分和生物过程进行GO功能分析以及在KEGG数据库中进行代谢途径的分析。然后使用RSEM计算各个样品的基因表达水平。然后使用R软件里的princomp函数进行PCA分析。

根据需求，使用NOIseq和PossionDis方法进行差异基因检测。NOIseq(参数：Fold Change ≥ 2.00；Probability ≥ 0.8) 方法基于噪音分布原理，本试验根据Tarazona等(2011)描述的方法进行差异表达基因(DEG)检测。PossionDis(参数：Fold Change ≥ 2.00；FDR ≥ 0.001) 方法基于泊松分布原理，本试验根据Audic等(1997)中描述的方法进行DEG检测。将得到的差异表达基因再进行单独的聚类分析、GO功能分析、Pathway功能分析，并与已知雄性不育同源基因进行比对分析。

使用Illumina Hiseq平台对楸树雄性可育花芽材料FL-1、FL-2、FL-3以及雄性不育材料SL-1、SL-2、SL-3进行转录组测序，共测得27.18 Gb原始数据。过滤掉原始数据中不合格的reads后，过滤后的clean reads Q20都达到了96.84%以上，每组数据质量见表 1。对过滤后的数据使用Trinity进行de novo组装，聚类去冗余得到Unigene，共得到86 076个Unigene。将得到的86 076个Unigene在7大功能数据库进行注释，注释结果分布见表 2，最终总共有64 600(75.05%)个Unigene得到注释。

表 1 过滤后的reads质量统计^□ Tab.1 The filtered reads quality statistics

表 1 过滤后的reads质量统计□ Tab.1 The filtered reads quality statistics

表 2 Unigene功能注释结果统计 Tab.2 Unigene functional annotation statistical results

表 2 Unigene功能注释结果统计 Tab.2 Unigene functional annotation statistical results

采用Audic等(1997)的方法对3组测序的差异表达基因进行筛选，每组差异表达基因分布见图 1。差异表达基因定义为FDR(fail discovery rate)＜0.001并且倍数差异在2倍以上的基因。楸树雄性不育与雄性可育3组转录组测序结果分析获得的差异表达基因数分别为26 149，27 320，23 466，差异上调表达基因分别有13 979，13513，130 55个，差异下调表达基因分别有12 170，13 807，10 411个，3组测序筛选的结果相似。用噪音分布原理对楸树雄性不育与雄性可育进行差异表达基因的筛选，差异表达基因为10 419个，表明这些差异基因与楸树雄性不育的关系密切。在楸树雄性不育花芽中，有3 504个差异基因(33.63%)在花芽组织中下调表达，有6 915个差异基因在花芽组织中上调表达(66.37%)。

图 1 差异表达基因统计结果 Fig.1 The statistical results of differentially expressed genes

对楸树雄性不育以及可育的差异表达基因进行GO功能富集，图 2为差异表达基因GO聚类分析图。发现富集在生物过程中的Unigene有2 329个，且富集较多的是代谢过程和细胞过程亚类；富集在细胞组分中Unigene有1 538个，富集较多的2个亚类是细胞和细胞部分；分子功能富集到的Unigene数为1 361，在该类中结合作用和催化活性2个亚类富集最多。显著富集的(corrected P＜0.05) GO功能亚类9个(表 3)，分子功能中的生长素外排膜转运蛋白活性，生物过程中的种子休眠期的胚胎发育、种子发育、果实发育、繁殖器官结构的发育、繁殖系统的发育、体轴建成、生殖发育过程、繁殖过程。这些类别在楸树雄性不育花芽组织中发挥了重要作用。

图 2 差异表达基因GO聚类分析 Fig.2 GO clustering analysis of Differentially expressed genes 1-20：生物过程；21-31:细胞组分；32-43：分子功能。1.生物黏附；2.生物调节；3.组织或生物发生细胞组件；4.细胞过程；5.发育过程；6.生长；7.定位；8.细胞活动；9.代谢过程；10.多有机体进程；11.有机体多细胞进程；12.负调控生物过程；13.正调控生物过程；14.生物过程调节；15.生殖；16.生殖过程；17.应激反应；18.节律进程；19.信号传导；20.单有机体进程；21.细胞；22.细胞连接；23.细胞成分；24.细胞外区域；25.复杂大分子；26.膜；27.膜组件；28.附膜腔；29.细胞器；30.细胞器组分；31.共质体；32.抗氧化活性；33.结合活性；34.催化活性；35.电子载体活性；36.酶调节活性；37.分子转导活性；38.核酸结合转录因子活性；39.养分贮液囊活性；40.蛋白结合转录因子的活性；41.受体活性；42.分子结构活性；43.转录活性。 1-20: Biological process; 21-31:Cellular component; 32-43: Molecular function.1.Biological adhesion; 2. Biological regulation; 3. Cellular component organization or biogenesis; 4. Cellular process; 5. Developmental process; 6. Growth; 7. Localization; 8. Locomotion; 9. Metabolic process; 10.Multi-organism process; 11. Multicellular organismal process; 12. Negative regulation of biological process; 13. Positive regulation of biological process; 14. Regulation of biological process; 15. Reproduction; 16. Reproductive process; 17. Response to stimulus; 18. Rhythmic process; 19. Signaling; 20. Single-organism process; 21. Cell; 22. Cell junction; 23. Cell part; 24. Extracellular region; 25. Macromolecular complex; 26. Membrane; 27. Membrane part; 28. Membrane-enclosed lumen; 29. Organelle; 30. Organelle part; 31. Symplast; 32. Antioxidant activity; 33. Binding; 34. Catalytic activity; 35. Electron carrier activity; 36. Enzyme regulator activity; 37. Molecular transducer activity; 38. Nucleic acid binding transcription factor activity; 39. Nutrient reservoir activity; 40. Protein binding transcription factor activity; 41. Receptor activity; 42. Structural molecule activity; 43. Transporter activity.

表 3 差异表达基因的GO功能显著性富集可视化分析 Tab.3 Visual analysis of differentially expressed genes which are significant enrichment by GO function

表 3 差异表达基因的GO功能显著性富集可视化分析 Tab.3 Visual analysis of differentially expressed genes which are significant enrichment by GO function

以楸树雄性不育与雄性可育花芽组织的差异表达基因KEGG pathway聚类分析，差异表达基因映射到127个不同的生物途径，显著富集的Pathway (Q＜0.05) 有42条(表 4)，主要为：代谢途径、生物合成的次生代谢产物、剪接体、RNA转运、内吞作用、甘油磷酸酯代谢、淀粉和蔗糖代谢、醚脂类代谢、苯丙素类生物合成等。这些途径在楸树雄性不育花芽组织中起重要调节作用。

表 4 差异表达基因在KEGG pathway显著性富集可视化分析 Tab.4 Visual analysis of differentially expressed genes which are significant enrichment by KEGG pathway

表 4 差异表达基因在KEGG pathway显著性富集可视化分析 Tab.4 Visual analysis of differentially expressed genes which are significant enrichment by KEGG pathway

根据同源性比对分析，最终共比对出246个Unigene。将得到的Unigene进行COG功能分类(图 3)，大多数都能够在RNA加工和修饰，细胞周期控制、细胞分裂、染色体分区，转录中得到注释，在信号转导机制与复制、重组和修复中也有Unigene得到注释。将比较分析后筛选得到的Unigene映射到代谢通路当中，最终多集中在丙酮酸代谢通路(图 4)以及与油菜素内酯(图 5)相关的代谢通路当中。在丙酮酸的代谢途径中有丙酮酸转化为草酰乙酸的途径中(图 4加粗线条标示)富集结果最为显著；而在油菜素内酯的信号转导途径中，合成的油菜素内酯在对油菜素内酯不敏感相关受体激酶及油菜素内酯不敏感蛋白复合物(BRI1) 的作用(图 5加粗框标示)过程中富集结果显著。因此根据富集结果可以分析出丙酮酸代谢途径和油菜素内酯信号转导途径对楸树雄性不育的形成起着重要的作用。

图 3 同源基因COG功能分类 Fig.3 COG functional classification of homologous genes A. RNA加工和修饰；D.细胞周期控制、细胞分裂、染色体分区；K.转录；L.复制、重组和修复；R.一般功能预测；T.信号传导机制。 A. RNA processing and modification; D. Cell cycle control, cell division, chromosome partitioning; K. Transcription; L. Replication, recombination and repair; R. General function prediction only; T. Signal transduction mechanisms.

图 4 丙酮酸参与相关代谢途径 Fig.4 Pyruvate involved in the metabolic pathways 粗线条标示丙酮酸转化为草酰乙酸的代谢途径，丙酮酸→磷酸烯醇丙酮酸→草酰乙酸。 The thick lines indicate the metabolic pathways of pyruvate conversion to oxaloacetate: Pyruvate→Phosphoenol-pyruvate→Oxaloacetate.

图 5 油菜素内酯相关信号转导途径 Fig.5 Brassinosteroid related signal transduction pathways 加粗框标示的BRI1为油菜素内酯不敏感蛋白1。BAK1:油菜素内酯不敏感1受体激酶1; BKI1:油菜素内酯不敏感蛋白1激酶抑制剂1；BSK:油菜素内酯信号激酶；BSU1:苏氨酸蛋白磷酸酶激酶；BIN2:油菜素内酯不敏感蛋白2；BZR1/2:油菜素内酯耐受蛋白1/2；TCH4:低聚木糖葡糖基转移酶；CYCD3:细胞周期蛋白D3。 The BRI1 marked with the bold box is brassinosteroid insensitive protein 1. BAK1: Brassinosteroid insensitive 1 associated receptor kinase 1; BKI1:Brassinosteroid insensitive protein 1 kinase inhibitor 1; BSK: Brassinosteroid-signaling kinase; BSU1: Threonine-protein phosphatase; BIN2: Brassinosteroid insensitive protein 2; BZR1/2: Brassinosteroid resistant protein1/2; TCH4: Xyloglucosyl transferase; CYCD3: Cyclin D3.

植物雄性不育的性状是在雄蕊发育过程中由一系列极其复杂的与调控育性有关的基因在时间和空间上共同作用调控的。引起植物雄性不育主要与花粉形成过程中的调节基因有关，这些基因在植物花粉发育的各个过程，包括减数分裂异常(Zhou et al., 2011；Nonomura et al., 2004)、胼胝质代谢异常(Wan et al., 2011)、绒毡层发育异常(Jung et al., 2005；Li et al., 2006)、花粉壁发育异常(Shi et al., 2011)以及花药开裂异常(Stintzi，2000；Steiner-Lange et al., 2003)等其他过程中都起着极其重要的作用。

众所周知，可育花粉的形成与绒毡层的正常发育密切相关。绒毡层细胞发育异常将直接或间接导致花粉败育。张博等(2015)发现，楸树雄性不育的形成，是由于绒毡层出现高度液泡化，导致其提前解体。本研究结合差异表达基因的GO分析和Pathway分析以及同源基因的分析，寻找出与楸树雄性不育有关的差异基因显著富集通路。通过对差异表达基因的GO功能分析，楸树的生殖过程、生殖发育过程、生殖系统发育、生殖结构发育均有显著性富集，生长素流出跨膜转运蛋白活性也显著性富集。差异表达基因的Pathway显著性分析表明，代谢途径、次生代谢产物的生物合成、剪接体、RNA转运及甘油磷脂、淀粉和蔗糖的代谢富集的差异基因数目较多。在同源基因的比对分析基础上，将与已知相似度高的差异基因再次进行功能富集分析。最终显著性富集结果主要映射在丙酮酸代谢以及植物激素信号转导途径中，这些分析为阐明楸树雄性不育的调控途径提供了线索。

GO分析结果表明楸树雄性不育的形成与花发育过程中花繁殖器官的形成过程相关，才会形成不育的花粉。这个结果与本实验室在细胞学上观察的楸树不育花粉的形成是绒毡层的发育异常相一致。Pathway通路分析表明楸树雄性不育的形成与生化过程有着较为紧密的联系。在已经报道的雄性不育与绒毡层发育相关的基因中，有一类基因属于bHLH转录因子。这一类转录因子主要参与调控植物发育进程以及生化过程(Andrade et al., 2014；Andriankaja et al., 2014)；而且这类转录因子表达异常后绒毡层也会出现高度液泡化(Zhang et al., 2006；Jung et al., 2005；Moon et al., 2013；Liu et al.，2014)。因此，在对差异基因Pathway分析结果做进一步的验证时，应将潜在的bHLH转录因子当作进一步分析的对象。

在与雄性不育有关的同源基因的分析中，得到Unigene集中映射在丙酮酸代谢和油菜素内酯参与的信号转导途径中。丙酮酸在生物体的3大物质代谢过程中起着极其重要的作用。在丙酮酸的代谢途径中，它可以通过乙酰CoA和重要的代谢途径三羧酸循环实现3大物质的相互转换。而一旦丙酮酸代谢异常就有可能对与生殖过程有关的包括减数分裂、绒毡层发育等途径造成影响。丙酮酸参与的代谢网络途径复杂，但由丙酮酸到草酰乙酸的转化这一过程得到显著的富集。该过程也是生物体内糖代谢的关键步骤。植物激素转导过程中，油菜素内酯的合成至关重要。油菜素内酯广泛存在于包括花粉在内的许多花器官当中。在油菜素内酯参与调节的生物过程中，通过合成的油菜素内酯对油菜素内酯不敏感相关受体激酶和油菜素内酯不敏感蛋白复合物的作用，最终对DNA的转录表达产生影响，使得细胞伸长和细胞分裂出现异常从而影响花粉的正常形成。

结合本实验室对楸树雄性不育的细胞学观察(张博等，2015)以及对转录组数据的分析，对楸树雄性不育得出一个初步的结果：在楸树花药发育过程中，由于丙酮酸参与的代谢途径出现异常，导致油菜素内酯合成途径受到影响。合成的油菜素内酯是通过作用于细胞表面由BAK1和BRI1形成的复合物来进行信号转导的。油菜素内酯合成异常，使得由BAK1和BRI1形成的复合物活性受到影响，导致参与调控细胞伸长和分裂的基因表达异常，从而使绒毡层细胞出现液泡化提前解体，楸树小孢子在减数分裂时期出现异常，最终形成不育的花粉。楸树的花粉败育形成的雄性不育，在本试验中已经有了初步的结论，但绒毡层异常形成的雄性不育，其形成过程十分复杂，只有通过更加深入的探索和研究才能对楸树雄性不育的机制有更加清晰的了解。

为了研究楸树雄性不育形成的机制，对1株雄性不育的楸树花芽和另1株正常发育的楸树花芽进行转录组测序并对差异表达基因进行分析，结果显示显著性富集的差异表达基因多在与花的繁殖器官相关的途径、重要的代谢途径和与油菜素内酯相关的信号转导途径中。据此推断可能是在花粉发育过程中，由于重要代谢途径异常使得油菜素内酯的合成途径出现变化最终导致绒毡层的不正常发育，使得小孢子减数分裂异常形成花粉败育。本研究从转录水平初步分析了楸树雄性不育的相关代谢途径，获得的有关结果对楸树雄性不育的研究有着重要的参考价值。