米槠雄花发育相关的转录组和MADS-box家族基因表达特性分析

连 辉， 江淑珍，b， 熊远芳，b， 沈 军， 陈世品，①

(福建农林大学： a. 林学院， b. 园林学院 兰科植物保护与利用国家林业和草原局重点实验室， 福建 福州 350002)

在被子植物中，花是重要的生殖器官，也是进化过程中最具多样性的结构[1]。同一个体(雌雄同株)或不同个体(雌雄异株)的单性雄花和雌花可促进异花授粉和基因交流，从而减少近亲繁殖[2]。壳斗科(Fagaceae)含900多种，是重要的经济资源，对森林生态系统的动态平衡有重要作用[3]。单性花是壳斗科植物的主要繁殖策略[4]。在壳斗科植物中，不同属植物有着不同类型花序，如栗属(CastaneaMill.)和锥属(CastanopsisSpach)植物，前者具有雄性和两性的柔荑花序，后者具有雄性和雌性的花序[4]。在栗属植物中，板栗(CastaneamollissimaBlume)雌花中有12枚发育不完全的退化雄蕊，具有短的花丝和可见的花药[5]。在锥栗〔Castaneahenryi(Skan) Rehd. et Wils.〕雌花中也存在退化雄蕊，在雌蕊原基的发育过程中雄蕊原基也发育，随着雌蕊原基的分化完成雄蕊停止发育[6]。而在其他壳斗科种类的研究中，单性花由花分生组织产生，不能同时形成雌性和雄性的性器官原基[4]。因此，在壳斗科植物中，栗属植物由于雄蕊停止生长而形成单性花，而同科其他属种类的单性花不存在雌蕊或雄蕊败育的现象，花发育起始阶段就是单性的。在这些花中，性别决定基因通过控制花分生组织的发生和花器官的形成来调控性别分化[7]。目前，雌雄同株植物的花器官发育已被广泛研究，在大多数情况下，花器官分化主要由MADS-box转录因子控制[8]。

米槠〔Castanopsiscarlesii(Hemsl.) Hayata〕为壳斗科锥属植物，广泛分布于中亚热带海拔1 300 m以下的山地丘陵，为亚热带常绿阔叶林顶级群落之一[9]。米槠具有2种类型花序：雄性柔荑花序在当年生的叶轴上发育，雌花序在枝稍附近发育。目前，在米槠的雄花中尚未发现雌蕊败育的迹象，表明米槠雄花在花发育的起始阶段就是单性的。在壳斗科植物中，花发育起始阶段就是单性的现象较为普遍，如欧洲栓皮栎(QuercussuberLinn.)[10]。然而，关于壳斗科植物从起始阶段就是单性的雌、雄花发育机制的研究还不够深入。本文基于转录组测序技术，对米槠雄花的花芽期、半开期和盛花期3个时期进行转录组分析，对差异表达unigenes进行功能注释。通过系统发育分析，鉴定出米槠MADS-box家族基因，并对其在不同发育阶段雄花的表达量进行分析。以期进一步认识MADS-box家族基因在米槠雄花发育中的作用，为深入探究壳斗科植物单性花发育的分子机制提供参考。

1 材料和方法

1.1 材料

供试米槠雄花来自于福建省福州市鼓山风景区(东经119°23′04″、北纬26°03′13″，海拔402 m)米槠野生群落。树种鉴定参照《中国植物志》[11]和《福建植物志》[12]。于2019年3月至5月，持续观察米槠雄花发育时期，分别在米槠雄花的花芽期(花被片和雄蕊均未分化，F1期)、半开期(花被片和雄蕊分化完成，花被片微微张开，F2期)及盛花期(花被片张开并开始散粉，F3期)采集雄花花序，每个时期分别从3株植株上采集，即为3个生物学重复，共采集9个样品。取样后立即用液氮速冻，保存于-80 ℃冰箱中，用于后续转录组测序分析。

1.2 方法

1.2.1 转录组测序与组装 使用PureLinkTMPlant RNA试剂盒(美国Invitrogen公司)分别提取米槠3个时期9个样品的总RNA，再使用RNase-FreeDNase (美国Promega公司)去除DNA污染。采用质量体积分数1.25%琼脂糖凝胶电泳检测总RNA的完整性，使用NanoDrop 1000超微量分光光度计(美国Thermo Fisher Scientific公司)测定总RNA的浓度和纯度。使用寡聚体(dT)贴附磁珠对总RNA进行纯化，将纯化后的总RNA用片段缓冲液破碎成小块，利用随机的六聚体引物反转录产生第1链cDNA，接着合成第2链cDNA。双链DNA末端补平，5′端磷酸化，3′端进行接头。使用PCR扩增双链cDNA片段，用AMPure XP磁珠(美国Thermo Fisher Scientific公司)纯化扩增产物，纯化后的双链PCR产物经热变性，溶解于EB溶液中，经安捷伦2100生物分析仪(美国Agilent公司)质量控制验证后，溶解于EB溶液中的双链PCR产物经95 ℃热变性成单链，经寡核苷酸序列循环，最终得到单链环状DNA作为文库。采用第二代测序技术，并基于Illumina HiSeq测序平台，对这些文库进行双末端测序。使用Trinity软件[13]进行转录组从头组装，为确保RNA测序数据的准确性和可靠性，从raw reads中剔除低质量、接头污染以及未知碱基N含量过高的reads，仅对剩下的高质量测序数据(clean reads)进行统计分析。

1.2.2 转录组差异表达unigenes注释 采用FPKM值表示转录组差异表达unigenes的表达丰度。采用Benjamini-Hochberg校正法对原有假设检验得到的P值进行校正，最终采用FDR(错误发现率)作为差异表达unigenes筛选的关键指标。利用DESeq软件，以FDR<0.05、FC(差异倍数)≥2作为差异表达unigenes的筛选标准，差异表达unigenes的功能注释信息通过KEGG数据库(http：∥www.genome.jp/kegg/)和GO数据库(http：∥www. geneontology.org/)获得。使用BLAST软件将差异表达unigenes序列与GO数据库比对，进行注释，在获得每个差异表达unigene的GO注释后，使用WEGO软件(http：∥wego.genomics.org.cn/cgi-bin/wego/index.pl)对所有差异表达unigenes进行GO功能分类统计，通过KEGG数据库的注释信息得到基因的path-way注释。

1.2.3 MADS-box家族基因鉴定与表达量分析 在TAIR在线数据库(https：∥www.arabidopsis.org/index.jsp)下载拟南芥〔Arabidopsisthaliana(Linn.) Heynh.〕MADS-box家族蛋白序列。以拟南芥MADS-box家庭蛋白序列作为索引在米槠基因组中进行本地BLASTp检索。将获得的候选序列提交到NCBI-CDD网站(https：∥www.ncbi.nlm.nih.gov/cdd)和HMMER网站(https：∥www.ebi.ac.uk/Tools/hmmer/search/hmmsearch)，根据Pfam(PF00319)在米槠基因组数据库中筛选出具有完整保守结构域的成员作为米槠MADS-box家族基因。由于米槠type-Ⅰ型MADS-box家族基因在花发育过程中几乎不表达，且该家族基因在花发育过程中调控作用的研究较少，因此，未开展该家族基因的系统进化树及表达量分析。采用MEGA 7.0软件对米槠和拟南芥的type-Ⅱ型MADS-box家族基因的蛋白序列进行多序列比对后，采用邻接法(neighbor-joining method)构建系统进化树，对建成的系统进化树进行自展法系数(bootstrap)校正，设置1 000次重复检验，最后优化进化树。从转录组数据中获取米槠雄花3个时期MADS-box家族基因的表达丰度值(FPKM)，通过log2(FPKM+1)计算表达量，并使用TBtools软件中的HeatMap程序绘制基因表达量热图。

2 结果和分析

2.1 测序质量分析和差异表达unigenes筛选

2.1.1 测序质量分析 测序结果(表1)显示：米槠雄花3个时期样品的clean reads为628 800 000～643 700 000，高质量总碱基量为6.29～6.44 Gb，Q30碱基百分比均不小于90.67%，GC含量为45.09%～45.39%。

2.1.2 差异表达unigenes筛选 对花芽期、半开期和盛花期3个时期米槠雄花的转录组结果进行比对分析，结果(表2)显示：花芽期与半开期间有15 086个差异表达unigenes，其中，与花芽期相比，半开期表达量上调的unigenes有5 283个，表达量下调的unigenes有9 803个。半开期与盛花期间有6 770个差异表达unigenes，其中，与半开期相比，盛花期表达量上调的unigenes有3 536个，表达量下调的unigenes有3 234个。花芽期与盛花期间有16 150个差异表达unigenes，其中，与花芽期相比，盛花期表达量上调的unigenes有6 588个，表达量下调的unigenes有9 562个。

表1 米槠雄花转录组数据统计

表2 不同时期间米槠雄花转录组差异表达unigenes数

2.2 功能注释分析

2.2.1 GO功能注释分析 在GO数据库对米槠雄花的花芽期与半开期、半开期与盛花期以及花芽期与盛花期3组的差异表达unigenes进行比对和功能注释，结果见图1。

由图1可以看出：米槠雄花的花芽期与半开期、半开期与盛花期以及花芽期与盛花期3组分别有31 888、13 228和35 203个差异表达unigenes注释在分子功能、细胞组分和生物过程3大类47个亚类中。在分子功能的13个亚类中，3组差异表达unigenes主要注释在催化活性和连接2个亚类中，其中，注释在催化活性亚类的差异表达unigenes分别有4 411、2 001和4 876个；注释在连接亚类的差异表达unigenes分别有4 437、1 911和4 896个。在细胞组分的14个亚类中，3组差异表达unigenes主要注释在膜、膜要素、细胞和细胞器4个亚类中，其中，注释在膜亚类的差异表达unigenes分别有3 215、1 493和3 464个；注释在膜要素亚类的差异表达unigenes分别有2 997、1 399和3 206个；注释在细胞亚类的差异表达unigenes分别有2 594、997和2 915个；注释在细胞器亚类的差异表达unigenes分别有1 862、674和2 111个。在生物过程的20个亚类中，3组的差异表达unigenes主要注释在细胞过程、代谢过程、生物调节和应激反应4个亚类中，其中，注释在细胞过程亚类的差异表达unigenes分别有2 738、983和3 071个；注释在代谢过程亚类的差异表达unigenes分别有2 347、690和2 403个；注释在生物调节亚类的差异表达unigenes分别有1 033、398和1 121个；注释在应激反应亚类的差异表达unigenes分别有981、393和1 037个。

2.2.2 KEGG代谢通路分析 在KEGG数据库中对米槠雄花的花芽期与半开期、半开期与盛花期以及花芽期与盛花期3组的差异表达unigenes进行比对和功能注释，结果见图2。

由图2可以看出：花芽期与半开期、半开期与盛花期以及花芽期与盛花期3组分别有10 826、5 241和11 896个差异表达unigenes注释在KEGG代谢通路，共6大类21亚类，且主要注释在代谢和遗传信息处理2大类中。3组注释在代谢相关通路的差异表达unigenes最多，分别有6 580、3 414和7 130个，包括全局及概要图、碳水化合物代谢、其他次生代谢产物的生物合成、氨基酸代谢和脂质代谢等11个亚类，其中，注释在全局及概要图亚类的差异表达unigenes最多，分别有2 473、1 266和2 673个。3组注释在遗传信息处理相关通路的差异表达unigenes分别有2 060、843和2 379个，包括翻译、转录、复制和修复以及折叠、分类和降解4个亚类。3组注释在环境信息处理相关通路的差异表达unigenes分别有846、384和943个，包括信号转录和膜运输2个亚类。在生物系统相关通路中，3组分别有846、384和943个差异表达unigenes注释在环境适应性亚类中。在细胞过程相关通路中，3组分别有372、171和443个差异表达unigenes注释在运输和分解代谢亚类中。在人类疾病相关通路中，3组共有32、8和40个差异表达unigenes注释在内分泌和代谢性疾病以及耐药性：抗菌2个亚类中。

2.3 type-Ⅱ型MADS-box家族基因鉴定及表达量分析

利用MEGA 7.0对米槠和拟南芥type-Ⅱ型MADS-box蛋白构建系统进化树(图3)。结果显示：米槠type-Ⅱ型MADS-box家族基因有53个成员，可分为13类：C/D类(6个成员)、ANR1类(4个成员)、AGL15类(4个成员)、AP3/PI类(1个成员)、Bs类(1个成员)、MIKC*类(7个成员)、SVP类(14个成员)、FLC类(2个成员)、SOC1类(3个成员)、AP1/FUL类(3个成员)、AGL12类(2个成员)、AGL6类(2个成员)和SEP类(4个成员)。

： 花芽期与半开期的差异表达unigenes的代谢通路Metabolic pathway of differentially expressed unigenes between flower bud stage and half flowering stage； ： 半开期与盛花期的差异表达unigenes的代谢通路Metabolic pathway of differentially expressed unigenes between half flowering stage and flowering stage； ： 花芽期与盛花期的差异表达unigenes的代谢通路Metabolic pathway of differentially expressed unigenes between flower bud stage and flowering stage.

基于花芽期、半开期和盛花期3个时期米槠雄花的转录组数据对米槠type-Ⅱ型MADS-box家族基因进行表达分析，发现其在雄花不同发育阶段的表达具有差异性。3个时期米槠雄花type-Ⅱ型MADS-box家族基因的表达量热图(图4)显示：在C/D类同源基因中，SHP基因(Maker00006492)和AG基因(Maker00069044)表达趋势类似，在花芽期表达量较低，半开期和盛花期表达量升高；其余4个C/D类同源基因在雄花发育过程中几乎不表达。在AP3/PI类同源基因中，AP3基因(Maker00011188)随雄花发育表达量逐渐升高。在MIKC*类同源基因中， Maker00016147和Maker00054870在雄花发育过程中表达量较高，其余基因表达量较低或不表达；其中，Maker00016147在花芽期的表达量较低，在半开期的表达量升至最高，在盛花期的表达量降至半开期的一半；Maker00054870的表达趋势类似，在半开期的表达量最高，在盛花期的表达量略下降。在AP1/FUL类同源基因中，Maker00031743和Maker00060929为AP1基因，Maker00031743在花芽期的表达量最高，在半开期的表达量降至最低，在盛花期的表达量较半开期升高；Maker00060929在花芽期、半开期和盛开期的表达量均较高，且几乎没有变化；FUL基因(Maker00039576)在花芽期和盛花期的表达量较高，但在半开期的表达量却较低。在AGL6类同源基因中，AG6基因(Maker00067599)在雄花发育3个时期的表达量均较高。在SEP类同源基因中， Maker00050389、Maker00039517、Maker00060937和Maker00031685均在雄花发育过程中高表达；其中，SEP3基因(Maker00050389)的表达量最高，且在雄花发育3个时期的波动不大；Maker00039517和Maker00031685的表达趋势相似，在花芽期的表达量较高，在半开期和盛花期的表达量逐渐升高；Maker00060937在花芽期的表达量较低，在半开期和盛花期的表达量逐渐升高。其余6类同源基因在米槠雄花发育过程中几乎不表达。

： 米槠Castanopsis carlesii (Hemsl.) Hayata； ： 拟南芥Arabidopsis thaliana (Linn.) Heynh.

F1： 花芽期Flower bud stage； F2： 半开期Half flowering stage； F3： 盛花期Flowering stage.

3 讨 论

功能注释结果表明：通过与GO数据库比对，花芽期与半开期、半开期与盛花期以及花芽期与盛花期3组差异表达unigenes注释在分子功能、细胞组分和生物过程的47个亚类中。在分子功能中，注释在催化活性和连接亚类的差异表达unigenes最多，可见催化酶在米槠雄花发育过程中发挥着重要作用；在细胞组分中，差异表达unigenes主要注释在膜、膜要素、细胞和细胞器4个亚类中，在雄花发育过程中，花粉粒的成熟伴随着新的细胞膜、细胞壁以及细胞器的生成，推测差异表达unigenes与花粉粒的成熟相关。通过KEGG数据库比对，3组差异表达unigenes主要注释在代谢和遗传信息处理2大类中。注释在代谢相关通路的差异表达unigenes最多，主要注释在全局及概要图、碳水化合物代谢和其他次生代谢产物的生物合成3个亚类中。相关研究表明：植物成花过程中需要消耗大量的碳水化合物和蛋白质等[14]，因此，注释在代谢相关通路的差异表达unigenes对米槠雄花成花有重要作用。此外，KEGG代谢通路注释结果还显示：注释在遗传信息处理相关通路的转录亚类和环境信息处理相关通路的信号转录亚类的差异表达unigenes较多，表明转录因子对米槠雄花的基因表达有重要作用。

相关研究表明：转录因子家族在植物的生殖过程中发挥着重要作用，其中MADS-box家族是控制花器官形成的重要转录因子家族[15]。相关研究根据系统发育学和基因结构的特点，将MADS-box家族分为type-Ⅰ型和type-Ⅱ型[16，17]。通过对米槠雄花的花芽期、半开期和盛花期的基因表达量分析，发现type-Ⅰ型MADS-box家族基因在米槠成花过程中几乎不表达，印证了大多数植物中以type-Ⅱ型MADS-box家族基因为主的推断[18]。通过对米槠type-Ⅱ型MADS-box家族基因进行鉴定及表达量分析，发现其在米槠雄花不同发育阶段的表达具有差异性，其中，鉴定到AP1/FUL类同源基因3个，AP1基因2个，其中1个在花芽期和盛花期表达量较高，另外1个在3个时期的表达量均较高；FUL基因1个，在花芽期和盛花期表达量较高。研究发现，FUL基因在花序分生组织中表达，对花分生组织的起始分化有重要作用，可能参与调控开花时间，而AP1基因在花分生组织和萼片的形成中有重要作用[19]。因此，米槠FUL基因在花芽期和盛花期表达量较高，推测对米槠雄花花序分化和开花起重要作用，AP1基因在半开期和盛花期高表达，调控花被片的形成。鉴定到AP3/PI类同源基因1个，为AP3基因，表达量随着雄花的发育逐渐升高。AP3/PI类基因通常在花瓣和雄蕊中表达，对控制花瓣和雄蕊的形成有重要作用[19]。因此，在米槠雄花发育过程中，AP3/PI类同源基因对雄蕊和花被片的发育有重要作用。鉴定到C/D类同源基因6个，其中米槠AG基因表达量随着雄花的发育逐渐升高，主要在雄花发育后期表达。在欧洲栗(CastaneasativaMill.)和欧洲栓皮栎的研究中，AG基因主要在雄花发育后期表达，与花粉发育密切相关[4，10]。在米槠雄花发育过程中，AG基因的表达趋势与欧洲栗和欧洲栓皮栎相似，因此，AG类同源基因可能参与调控米槠的花粉发育。米槠SHP基因在时间上的表达模式与AG基因相似，主要在雄花发育后期表达。在欧洲栗和欧洲栓皮栎的研究中，SHP基因在雄花发育后期高表达，同时还在雌花中表达，表明壳斗科的SHP类同源基因可能是功能性的C类基因[4，10]。而STK基因主要参与调控胚珠的发育[20]，在米槠雄花发育的3个时期均未表达，印证了米槠雄花中不存在雌蕊的结构。鉴定到SEP类同源基因4个，在米槠雄花发育的3个时期有较高表达，表明SEP类同源基因在米槠雄花发育过程中起协同作用。

综上所述，在米槠雄花发育过程中，AP1/FUL类、AP3/PI类、C/D类和SEP类同源基因有较高的表达，且在不同时期的表达具有差异性，而这些基因在其他物种中已被证实参与调控各个花器官的发育。米槠是雌雄同株的单性花植物，AP3/PI类和C/D类同源基因对雄蕊和雌蕊的形成有着重要的作用。因此，下一步可重点分析MADS-box家族基因分别在米槠雄花和雌花不同花器官中的表达情况，以期阐明这些基因在米槠性别分化和花器官发育中的功能和作用机制。