花药是雄蕊重要的组成部分，含有与花粉粒形成及释放有关的生殖和营养组织^[1]。花药开裂作为花药发育过程的最后一个阶段，能否适时开裂直接关系到花粉能否及时到达柱头、植物能否完成授粉受精^[2]。植物功能雄性不育是雄性不育的一种类型，表现为开花时花药不开裂，而花粉具有活力。茄子Solanum melongena L.是世界各地广泛栽培的蔬菜作物之一，具有明显的杂种优势，并已在杂交种生产中得到广泛应用，而生产上主要采用人工去雄授粉的方法获得F1代种子。由于人工去雄和授粉制种过程繁琐，生产成本较高，加上去雄时选择花蕾不当或去雄不彻底，因而很容易形成假杂种，降低种子的纯度和质量，给生产带来一定的损失^[3]。利用雄性不育系生产茄子杂交一代种子不仅可以简化育种程序，提高育种效率，还可以提高种子质量、降低生产成本，具有显著的经济效益和社会效益。植物功能雄性不育系因自身有可育花粉，无需再选育保持系，在育种上转育省事^[4]，因此开展功能雄性不育的研究具有广阔的应用前景。

从Jasmin^[5]首次认定茄子功能雄性不育材料开始，育种家们陆续开展了茄子功能雄性不育系的选育及遗传机制等研究，获得了一些不育材料^[6-9]，但与大白菜^[10]、甘蓝^[11]、番茄^[12]、辣椒^[13]等其他蔬菜作物雄性不育的研究相比还较滞后，尤其对其不育机理的研究较少。转录组测序(RNA- sequencing)具有灵敏度高、信号数字化、检测范围广等优点，被广泛用于生物体的多种功能研究^[14]。随着高通量转录组测序技术的发展与成熟，已采用高通量测序技术对水稻Oryza sativa、小麦Triticum aestivum、玉米Zea mays等20多个物种的细胞核雄性不育、细胞质雄性不育、光敏雄性不育等多种败育类型进行了雄蕊败育发生特征的研究，这些研究对于深入解析不同类型雄性不育的发生机制具有重要意义^[15]。

本研究以茄子功能雄性不育系和可育系为材料，通过组织切片和转录组分析，以期找到茄子功能雄性不育花药发育的组织结构特点和差异表达转录组信息，从而为进一步解析茄子花药开裂调控机制，加速茄子雄性不育新品种选育奠定一定的理论基础。

茄子可育系F142和功能雄性不育系S13由重庆市农业科学院蔬菜花卉研究所提供。取开花当天的花蕾，TTC法检测花粉活力。分别取F142和S13开花前8、7、6、5、4、3、2、1天和开花当天9个时期的花蕾，迅速剥离花药，FAA固定，按常规石蜡切片程序脱水、透明、浸蜡、包埋、切片(厚度8 μm)、脱蜡、染色(1%硫酸铁铵媒染22 min，0.1%的苏木精染色30 min)、分色(0.5%的硫酸铁铵分色1 min)，加拿大树胶封固，自然风干。Leica荧光显微镜Y35白光通道下观察并照相。

取F142和S13开花前8天、5天及开花当天的花蕾分别编号、锡箔纸包装，投入液氮保存，送至北京百迈克生物科技有限公司进行RNA提取和质量检测，RNA样品采用NanoDrop和Agilent Technologies 2100 Bioanalyzer检测。使用带有Oligo (dT)的磁珠富集mRNA；加入片段化缓冲液将mRNA打断成短片段，以打断后的mRNA为模板合成一链cDNA，然后加入缓冲液、dNTPs、RNase H和DNA聚合酶Ⅰ配制二链合成反应体系合成二链cDNA，经过PCR扩增，建立测序文库；构建好的文库经Agilent 2100 Bioanalyzer和ABIStepOnePlus Real-Time PCR System质检合格后，使用Illumina HiSeq^TM 2000进行测序。

对测序获得的Raw Data进行数据过滤，去除其中的接头序列及低质量Reads，获得高质量的Clean Data。之后，将其与参考基因组(Solanum melongena L.，ftp://ftp.solgenomics.net/genomes/ Solanum_melongena/)进行序列比对，获取在参考基因组或基因上的位置信息，以及测序样品特有的序列特征信息。使用DEGseq软件进行样品组间的差异表达分析，以FDR < 0.01且|log₂ (FC)|≥2作为标准，筛选出差异表达基因(Differentially expressed genes，DEGs)。将获得的DEGs用Blastx软件在8大功能数据库(COG库、GO注释库、KEGG库、KOG库、NR库、Pfam库、SwissProt库和eggNOG库)中进行比对，获得DEGs的主要代谢通路、功能注释以及生物学过程。

利用R软件中的WGCNA软件包^[16]完成加权基因共表达网络构建。首先，以每对基因的相关系数构建基因共表达矩阵，共18个转录组样本(两种材料，3个时期，各3次重复)。根据无尺度网络的原则来确定加权系数β。经过阈值筛选，选择β=18对原本基因共表达相关矩阵进行幂处理，得到无尺度化邻接矩阵^[17]。之后结合某个基因和分析中的其他所有基因之间的关系，将邻接矩阵转化为拓扑矩阵，并利用节点的相异程度进行拓扑重叠构建基因联系网络。按照动态剪切算法进行基因聚类及模块划分。设置每个模块内最少基因数目为30^[17]，再依次计算每个模块的特征向量值ME (module eigengene)，选择的相似模块合并为新的模块，相似模块合并阈值为0.25 (cutHeight=0.34)。

通过cor函数量化模块特征向量ME与不同样本之间的相关系数。相关系数越大，模块内基因与样本特征关联度越高，基因相对表达量越高，而相关系数值越小，模块内基因与样本特征关联度越低，基因相对表达量越低。选择与样本特征具备较高相关系数的模块作为备选特异性模块^[18]。

按照Trizol试剂说明书提取总RNA，然后按照反转录试剂盒(天根生化科技北京有限公司)说明书合成cDNA。以cDNA为模板，对筛选出的差异基因的转录水平进行实时荧光定量PCR分析。GAPDH作为内参基因，每个样品进行3次技术重复。利用2^–ΔΔC_t方法计算基因的相对表达量。

观察发现，茄子可育系F142的叶片为长卵圆形，盛开的花瓣为浅紫色，花柱被5–7枚花药包围，花药呈黄色，开花后花药顶端的小孔张开，花粉散出(图 1，A1、A2)；不育系S13的叶片为长卵圆形，边缘波状，盛开的花瓣为浅紫和白色相间，花柱周围花药6–9枚，黄色花药的边缘纵贯一条深褐色线条，开花后花药顶孔紧闭，不开裂，花粉囊饱满(图 1，B1、B2)。花粉活力测定显示，可育系F142的成熟花粉活力为87% (图 1，A3)，不育系S13的花粉活力为83% (图 1，B3)，表明茄子功能雄性不育系的花粉活力与可育系相近。

图 1 茄子可育系与不育系的花、花药顶端及花粉活力 Fig. 1 The flower, anther tips and pollen vitality of fertile eggplant and male sterile eggplant. A1, A2, A3: flower, anther and pollens of F142; B1, B2, B3: flower, anther and pollens of S13.

图选项

石蜡切片观察发现(图 2)：在开花前8天，茄子不育系S13 (图 2，B1)与可育系F142 (图 2，A1)花药的组织结构没有明显区别，花药的横切面由外至内可见表皮、药室内壁、中层、绒毡层以及裂口组织凹口处的环形细胞簇，孢原细胞分布于药室中。开花前7天–开花前5天，F142药隔组织细胞与花药壁相连，可育系药隔横向扩增(图 2，A2)，药室呈马蹄状分布，裂口组织细胞呈纵向紧密排列(图 2，A3)，凹口处的环形细胞簇几乎完全解体，同侧药室已连通(图 2，A4)；而S13的药隔晚1天开始横向扩增(图 2，B3)，药室呈不明显的马蹄状，药隔与花药壁仍处于连接状态，环形细胞簇结构完整(图 2，B4)；花前4天–花前1天，F142花药中的环形细胞簇完全解体(图 2，A5)，药隔增大(图 2，A6)，裂口组织从缢缩(图 2，A7)到出现轻微解体(图 2，A8)；而S13花药中的环形细胞簇开始出现解体(图 2，B5)，与F142花前6天(图 2，A4)相似，但药隔较可育系大，药室裂口组织出现缢缩(图 2，B8)；开花当天，F142裂口组织断裂层明显，内绒毡层全部解体，花粉囊腔增大，由于药孔开裂散粉，囊腔中只存有少量成熟的花粉粒(图 2，A9)，而S13变化不明显(图 2，B9)，未见内绒毡层解体和裂口组织断裂，花粉囊腔内有大量成熟花粉粒。从花药发育过程可以发现，环形细胞簇的延迟解体及开花当天裂口细胞的不解体与花药的不开裂有关。

图 2 茄子花蕾发育各时期的花药横切图 Fig. 2 The cross-sectional diagram of each period in bud development. A1–A9: the nine stages of flower development of F142; B1–B9: the nine stages of flower development of S13; Ep: epiderm; En: endothecium; ML: middle layer; T: tapetum; ST: stomiun tissue; CCC: circular cell clusters; IT: inner tapetum; bar=100 μm.

图选项

通过对可育系F142和功能雄性不育系S13花蕾发育3个时期的转录本进行比较分析，得到了7 585个差异表达基因(DEGs)，花前8天、5天和开花当天分别有4 105个(1 762个上调，2 343个下调)、5 765个(2 739个上调，3 026个下调)、2 736个(1 340个上调，1 396个下调)差异表达基因(图 3A)。有1 436个基因为两种材料的3个花蕾发育阶段共有的DEGs (图 3B)，其中651个差异表达上调基因(图 3C)，785个差异表达下调基因(图 3D)。

图 3 茄子可育系F142和不育系S13花蕾不同时期的差异表达基因分析 Fig. 3 Analysis of DEGs between the fertile line F142 and male sterile line S13 in different stages of flower buds. (A) The total number of up-regulated and down-regulated DEGs. (B) Venn diagram of all DEGs. (C) Venn diagram of up-regulated genes. (D) Venn diagram of down-regulated genes.

图选项

基于基因表达的相似性，对在F142和S13的3个花蕾发育阶段都表达的1 436个DEGs进行层次聚类分析。结果表明，大多数DEGs在F142和S13的表达模式完全相反。大部分在F142中表达水平较高的基因在S13中表达水平较低，反之亦然(图 4A)。这1 436个DEGs被分为7类基因表达趋势(图 4B)。在Subcluster 1、Subcluster 2、Subcluster 5和Subcluster 7中的DEGs在F142中有较高的转录水平，在Subcluster 3、Subcluster 4和Subcluster 6中的DEGs在S13中有较高的转录水平。

图 4 茄子可育系F142和不育系S13差异表达基因聚类分析 Fig. 4 Clustering analysis of DEGs between the fertile line F142 and male sterile line S13. (A) Hierarchical clustering of the 1 436 DEGs. (B) Expression patterns of DEGs in the seven clusters. F142-1, F142-2 and F142-3 represent F142 at 8 days before flowering, 5 days before flowering and on flowering day, respectively. S13-1, S13-2 and S13-3 represent S13 at 8 days before flowering, 5 days before flowering and on flowering day.

图选项

对样本中的1 436个差异表达基因进行GO显著性富集分析，按照基因的分子功能(Molecular function)、参与的生物过程(Biological process)和所处的细胞位置(Cellular component)对注释的基因进行分类，明确其在各GO项中的分布情况。这三大功能分类可分为30个亚类，其中生物学过程有13个亚类，细胞组分有10个亚类，分子功能有7个亚类(图 5A)。其中在生物学过程中涉及单细胞生物过程(Single-organism process)、代谢过程(Metabolic process)和细胞过程(Cellular process)的Unigene较多，分别有385条、557条和561条；细胞组分中涉及较多的是细胞部分(Cell part)和细胞(Cell)，分别有600条和603条；分子功能中涉及较多的有催化活性(Catalytic activity)和结合(Binding)，分别有570条和593条。KEGG注释表明，差异基因主要富集在次生代谢产物的生物合成(Biosynthesis of secondary metabolites)、代谢途径(Metabolic pathways)、内质网蛋白质加工(Protein processing in endoplasmic reticulum)、氨基酸的生物合成(Biosynthesis of amino acids)、碳代谢(Carbon metabolism)和植物激素信号转导(Plant hormone signal transduction) (图 5B)。

图 5 茄子花蕾不同时期差异表达基因GO和KEGG分析 Fig. 5 Analysis of GO enrichment and KEGG pathway enrichment. (A) Analysis of GO enrichment for the common DEGs. (B) Analysis of KEGG pathway enrichment for the common DEGs.

图选项

过滤去除表达矩阵中变异较低的基因，最终选取16 465个基因进行加权基因共表达网络构建分析，确定加权系数β为18，利用所有基因间的相异系数为基本元素构建分层聚类树，并使用动态剪切算法对基因进行模块划分，再通过计算各个模块的模块特征向量，对模块划分进行优化合并，最终共得到15个基因的共表达模块(图 6A)。构建每个模块与样本特征的相关性热图和模块内基因数量的条形图，热图中颜色越深表示模块与特征值的相关性越高(图 6B)。经筛选发现，Plum2模块与开花前5天的S13存在高度正相关(r=0.97，P=2E-11)，Royalbiue模块与开花前8天的S13存在高度正相关(r=0.81，P=4E-05)，Bisque4模块与开花当天的S13存在高度正相关(r=0.99，P=4E-15)。

图 6 功能雄性不育茄子相关基因的权重基因共表达网络分析 Fig. 6 Weighted correlation network analysis of functional male sterile-related genes. (A) Gene co-expression network gene clustering number and modular cutting. (B) Heat map of the correlation between modules and traits.

图选项

将与S13三个花蕾发育时期相关的Plum2、Royalblue和Bisque4模块的基因表达量进行聚类构建热图(图 7A、B、C)。图中绿色越深表明基因的相对表达量越低，红色越深则表达量越高，可以看出这3个模块在S13中上调表达。将3个模块内的基因进行KEGG代谢通路富集(图 7D、E、F)。Plum2模块主要富集到α-亚麻酸代谢(Alpha-linolenic acid metabolism)、多糖的生物合成和代谢(Glycerophospholipid metabolism)以及脂肪酸降解(Fatty acid degradation)；Royalblue模块主要富集到苯丙烷类生物合成(Phenylpropanoid biosynthesis)、光合作用(Photosynthesis)以及卟啉和叶绿素代谢(Porphyrin and chlorophyll metabolism)；Bisque4主要富集到戊糖和葡萄糖醛酸的相互转化(Pentose and glucuronate interconversions)。

图 7 模块的聚类分析和KEGG富集分析 Fig. 7 Cluster analysis and KEGG pathway enrichment analysis of modules. (A–C) Cluster analysis of the genes in three modules related to S13 in 18 samples of expression profile. (D–F) KEGG enrichment analysis of modules Plum2, Royalblue and Bisque4.

图选项

基于差异表达基因的GO和KEGG分析及鉴定得到的与S13花蕾发育高度相关的特异性模块，从中挑选9个基因进行实时荧光定量PCR (qRT-PCR)验证。从图 8可以看出，丝氨酸羟甲可以看出这3个模块在S13中上调表达。将3个模块内的基因进行KEGG代谢通路富集(图 7D、E、F)。Plum2模块主要富集到α-亚麻酸代谢(Alpha-linolenic acid metabolism)、多糖的生物合成和代谢(Glycerophospholipid metabolism)以及脂肪酸降解(Fatty acid degradation)；Royalblue模块主要富集到苯丙烷类生物合成(Phenylpropanoid biosynthesis)、光合作用(Photosynthesis)以及卟啉和叶绿素代谢(Porphyrin and chlorophyll metabolism)；Bisque4主要富集到戊糖和葡萄糖醛酸的相互转化(Pentose and glucuronate interconversions)。

图 8 差异表达基因的qRT-PCR验证 Fig. 8 qRT-PCR validation of DEGs.

图选项

对于开花植物，花药开裂是植物有性繁殖的先决条件。植物花药开裂释放出成熟花粉是由一系列的细胞解体过程组成，包括绒毡层、隔膜和最后的裂口组织^[19]。研究者通过对大麦^[20]、水稻^[21]、拟南芥^[22]和番茄^[23]等不同植物花药开裂过程中花药的解剖学和分子生物学研究发现，植物花药开裂是一个复杂过程，并且不同物种花药的开裂存在不同的调控机制。Keijzer等^[24]指出裂口组织是花药开裂的位置，是由特殊细胞构成的一层单细胞层；在番茄花药开裂过程中，发现绒毡层细胞解体后，药室外壁细胞加强，在裂口处产生了最终引起花药释放的力量。茄科中还有一种特殊的组织――环形细胞簇，位于裂口组织的凹口处，Beals等^[25]发现环形细胞簇与花药的开裂有关。本研究发现，在茄子花药发育过程中可育系花药中环形细胞簇的解体发生在开花前第6天，而功能雄性不育系花药环形细胞簇的解体时间比可育系推迟2 d，解体程度也较可育系弱；在开花当天，可育系的裂口细胞及内绒毡层细胞死亡并解体，而不育系中却没有发生此现象，且内绒毡层细胞依然存在。因此推测花药发育过程中环形细胞簇的轻微解体以及开花当天裂口细胞的不凋亡可能是导致茄子功能雄性不育系花药不开裂的原因。

花药开裂涉及到花药壁脱水、隔膜的酶解、植物激素的调控等^[1]。拟南芥中黄酮转运因子FFT可通过调节花药脱水来控制花药开裂，在fft-1突变体中花药不开裂^[26]。Ogawa等^[27]研究发现拟南芥开裂区域的多聚半乳糖醛酸酶ADPG1、ADPG2和QRT2均与花药的开裂有关。茉莉酸(Jasmonic acid，JA)和生长素协调控制拟南芥雄蕊的晚期发育^[28]，JA生物合成或信号转导突变体的花药通常不能正常开裂^[29]。JA信号转导途径关键基因SmCOI1在茄子功能雄性不育系中的表达量比可育系低^[30]。除JA外，生长素在花药开裂中也起着重要作用^{[21, 31]}，ARF6和ARF8的缺失影响茉莉酸的形成，导致花药延迟开裂^[32]，本研究发现，ARF5只在功能雄性不育系S13中表达，而其在茄子花药开裂中的分子功能还有待进一步研究。MYB转录因子在花药的发育中起着重要的调控作用^[33]。Yang等发现，AtMYB26蛋白存在于拟南芥的花药内，通过MYB26、NST1、NST2之间的一系列调控，可使次生增厚在花药内形成，推测AtMYB26蛋白在花药开裂中起关键作用^[34]。此外，MYB转录因子还可通过茉莉酸、生长素、赤霉素等激素途径影响花药发育。本研究也发现MYB26在茄子可育系与不育系的花药发育过程中都有不同程度的表达，而对于MYB转录因子如何参与植物激素信号转导途径调控茄子花药开裂将是进一步研究的重点。

通过转录组分析，即使是无参考基因组的物种，也可以获得在某一状态下基因的表达丰度、剪接方式等信息^[15]。本研究通过对开花前8天、5天及开花当天的可育系F142和功能雄性不育系S13花蕾进行转录组测序，筛选出了7 585个差异表达基因，其中有1 436个基因在两种材料的3个花蕾发育阶段都差异表达。利用GO分析发现相关差异表达基因富集在单细胞生物过程、细胞过程和代谢过程，差异基因的分子功能主要富集于催化活性和结合；KEGG通路富集分析表明，差异表达基因主要富集于次生代谢物的生物合成、代谢途径、内质网蛋白质加工、氨基酸的生物合成、碳代谢和植物激素信号转导等代谢途径，这也暗示了茄子功能雄性不育形成机理的复杂性。Yang等^[35]对茄子雄性不育系EP26A与其保持系EP26进行转录组比较分析发现在花粉发育的3个阶段共有1 716个共同差异表达基因，两种材料的DEGs主要参与氧化还原、碳水化合物和氨基酸代谢及转录调控，推测相关基因和通路可能与雄性不育有关。此外，茄子反向温敏核不育系05 ms与其同源可育系S63的差异表达基因也富集到次生代谢物的生物合成、代谢途径、碳代谢、苯丙烷类生物合成和植物激素信号转导等代谢途径^[36]，这些通路都参与能量和营养物质的供应和分配，从而影响细胞和细胞壁的发育。在棉花中发现，高温胁迫会导致棉花花粉活性降低及花药不开裂，高温耐受型和高温敏感型两个棉花品系不同花药发育时期受高温诱导后表达趋势不同，差异表达基因主要富集在碳水化合物代谢、植物激素信号和内质网蛋白质路径^[37]。在本研究中功能雄性不育系S13差异表达基因没有富集到淀粉和蔗糖代谢途径，而是富集到内质网蛋白质加工途径，推测S13的花药不开裂与内质网蛋白质加工相关。对这些内质网蛋白质加工相关基因的进一步研究将有助于揭示茄子花药开裂的分子调控机制。

利用WGCNA构建基因共表达网络，通过基因的层次聚类，共划分出15个共表达模块。每个模块代表了具有相同表达模式的基因集。根据模块与特征性状的相关性关联分析，筛选出与功能雄性不育系S13三个花蕾发育时期高度相关的3个模块：Plum2、Royalblue和Bisque4模块。并且3个模块的基因在S13中都上调表达。对模块内的差异表达基因进行KEGG代谢通路富集发现S13花蕾发育与苯丙烷类生物合成、光合作用、卟啉与叶绿素代谢、α-亚麻酸代谢、多糖的生物合成和代谢、脂肪酸降解以及戊糖与葡萄糖醛酸的相互转化等代谢过程相关。这些差异表达基因和筛选出的关键代谢途径为进一步研究功能雄性不育系茄子的育性转换以及花药开裂分子调控机理提供了理论依据和参考。