芝麻显性核不育系可育和不育花芽的比较转录组分析(胡麻属indicuml .)

背景

芝麻(胡麻属indicum油料作物是全球重要的油料作物，含油量高。该作物常表现出较强的杂种优势，可用基因雄性不育(GMS)手段加以利用。我们已经开发了一个优势GMS (DGMS)系W1098A_1.混合动力车。虽然已经对其进行了遗传学和解剖学研究，但其雄性不育的分子机制尚不清楚，因此限制了该GMS系的充分利用。本研究采用基于RNA-seq的转录组分析方法，对两个近同基因的DGMS株系(W1098A和W1098B)进行了雄性不育相关差异表达基因(DEGs)的鉴定。

结果

共检测到1502个显著的DEGs，其中751个在不育花芽中上调，751个在不育花芽中下调。许多DEG参与乙烯和JA的合成和信号通路；其表达在不育芽中分别上调或下调。此外，大多数与DEGs相关的NAC和WRKY转录因子在不育芽中上调。通过查询植物雄性生殖数据库，获得49个芝麻同源基因；其中一些编码转录因子（bHLH089、MYB99和AMS），这些转录因子在不育芽中表达降低，因此暗示可能在指定或决定绒毡层命运和发育中发挥作用。等位基因变体对其相应DEG功能的预测效应强调了几个插入/缺失（INDEL），这可能是DGMS系不育/育性表型的原因。

结论

目前的转录组比较研究表明，激素信号通路和转录因子共同控制着芝麻DGMS雄性不育。结果还表明，位于DEGs中的几个InDels容易导致功能丧失，这可能是导致男性不育的原因之一。这些发现为深入了解DGMS的分子机制提供了宝贵的基因组资源。

芝麻(胡麻属indicumL.）是全球重要的，古代油籽作物，主要用于高质量的石油[1.,2.]．它是油料作物中含油量最高的，富含天然抗氧化剂，如芝麻素和芝麻醇，它们因其独特的降压作用和抗氧化活性而为人所知[3.–5.]．芝麻的种子产量虽然重要，但与油菜籽、花生和大豆相比，芝麻的种子产量不稳定且相对较低。因此，应大力提高芝麻的种子产量。

杂种优势利用是用于提高产量的最有前途的方法，因为非常强的杂交优势（> 15％）已经被这种作物内观察到6.]．细胞质雄性不育(CMS)和基因雄性不育(GMS)都可以有效地利用杂种优势。到目前为止，只有隐性GMS被成功应用于芝麻F的生产_1.混合动力车。然而，该方法可能受到环境敏感性，不完全无菌等某些缺点的限制，并且及时除去50％雄性肥沃的植物，用于杂种种子的两种型线[7.]．近年来，我们通过与野生植物杂交，获得了一种新的优势GMS株系(DGMS)S. mulayanumL. (2n = 26)植物与栽培种美国indicumL. (2n = 26)，具有很大的杂交育种潜力。细胞学研究表明，DGMS株系(W1098A)的花粉败育始于花粉母细胞(PMC)，并在花粉发育过程中持续进行，在小孢子晚期达到顶峰。利用两个紧密相连的SSR标记SBM298和GB50 [8.]．然而，潜在的分子机制仍不清楚。

小二倍体基因组（〜350 MB）使芝麻成为遗传研究的有吸引力的物种[9,10]．最近，组装了芝麻的高质量基因组序列，其含有〜27,148个预测的基因模型，其中91.7％锚定到16个假调子或连杆（LGS）上锚定[11]．利用正向和反向遗传方法，越来越多的基因已被确认在花药发育中有重要作用。因此，植物雄性生殖数据库(PMRDhttp://202.120.45.92/ADDB/)，一个与植物雄性生殖有关的基因和突变体的综合资源已出现[12]．

雄性不育（MS）不仅与缺乏可行的花粉有关，还与花粉释放的失败相关联13]．绒毡层程序性细胞死亡(PCD)对成功形成花粉的重要性已经被一些未能通过正常绒毡层分解的MS突变体所强调[13–15]．原孢细胞数量和绒毡层细胞命运由EXCESS MICROSPOROCYTES1 (EMS1)控制，EMS1是一种富含亮氨酸重复受体样激酶，和一个小分泌蛋白配体，绒毡层DETERMINANT1 (TPD1) [16]．绒毡层的发育是由功能失调的TAPETUM1(DYT1) [17),绒毡层发育和功能有缺陷(TDF1) [18]具有绒毡层成熟、花粉壁形成和绒毡层PCD夭折的小孢子(自动对盘及成交系统) [19),男性STERILITY1(MS1) [20.]．裂解的最后阶段涉及茉莉酸(JA)诱导的基因表达和与内皮次生增厚相关的转录因子[13]．

为了更全面地阐明MS的机制，许多高等植物的转录组已经被测序，包括拟南芥[21],荞麦[22)、棉花(23–25,西瓜26,大豆27]，芸苔栗鸟［28–30.),芸苔属植物oleracea［31]．本研究以长度为2.5 mm的DGMS株系的可育和不育花蕾为样本，进行RNA-seq测序，这是对芝麻DGMS转录组的首次研究。本研究的目的是鉴定与MS相关的差异表达基因(DEGs)，并探索不同的生物过程及其假定的功能。这些结果将有助于阐明DGMS的分子机制，并有助于芝麻杂交种的选育。

可育芽和不育芽的转录组分析

我们之前已经证实DGMS系的雄性不育主要发生在PMC期[8.]．因此，我们对这一阶段的可育芽和不育芽进行了取样，并分别制备了cDNA文库。通过Illumina HiSeq 2000平台测序，我们分别从可育花芽和不育花芽中获得了53,126,890和55491,408个高质量的对端reads，然后对其进行清洗，并将其定位到包含27,148个基因模型的芝麻参考基因组序列[11]．从可育芽中83.54%的reads和从不育芽中84.86%的reads被定位到参考基因组中，其中大部分为唯一定位(见表)1.）.通过序列比对发现，从可育芽和不育芽中获得的unique reads分别命中了22,373和22,788个基因，占已知基因模型>的82%。可育芽平均基因长度为1305 bp，不育芽平均基因长度为1297 bp。大部分基因(74%在不育芽中，71%在不育芽中)表现出很高的基因覆盖率(90 ~ 100%)。

为了测量基因在不同组织中的相对表达水平，我们基于唯一的定位Reads计算了RPKM(每Kilobase的外显子模型每百万mapping Reads)值。这些基因的RPKM值为0.012 ~ 16683.020，平均值为40.974。不育芽的最小RPKM为0.008，最大RPKM为33521.52，平均RPKM为40.302。因此，上述所有基因要么在可育芽中表达，要么在不育芽中表达，RPKM阈值≥0.001。不出所料，大多数表达基因(> 95%)在组织间是共同的;然而，我们也发现了少量的独特表达基因(539个在可育芽中表达，954个在不育芽中表达)。

二基因的功能表征

以至少2倍变化和错误发现率(FDR)<0.001为标准，通过比较可育芽和不育芽之间的基因表达水平，我们获得了1,502个显著的DEGs，其中不育芽中有751个表达上调，751个表达下调(附加文件)1.:表S2)。通过锚定之前释放的16个假分子(或lg)的基因序列，分析了芝麻基因组中所有deg的分布，这些假分子包含了85.3%的芝麻基因组组装[11]．通过整合公共领域的基因组信息，我们可以将DEGs分配到每个LG上。结果表明，LG4的deg最少(4.47%)，LG11次之(4.76%)。LG7的deg比例最高(6.83%)。在LG16、LG8和LG15中，上调基因的比例是下调基因的近2倍。而LG3、LG4、LG5、LG9、LG11、LG12则呈相反趋势。此外，在其余四种LGs中，向上和向下调节的基因数量几乎相同(图。1.）.

然后用GO (Gene Ontology)和KEGG (Kyoto Encyclopedia of Genes and Genomes)数据库对每个DEG假定的功能进行了表征。由于氧化石墨烯类的数量多且分支结构复杂，因此仅以最丰富的三个官能团为例，即“细胞组分”、“分子功能”和“生物过程”(图1)。2.）.In the sub-category of ‘Cellular Component’, the largest numbers of genes were found to be associated with ‘cell part’, which can be further sub-divided into cascades of ‘intracellular’, ‘cytoplasmic vesicle’ and ‘intrinsic to membrane’. In the next main sub-category of ‘Molecular Function’, ‘ion binding’ and ‘catalytic’ were the most abundant cascades that have a respective of 71 and 19 genes. Moreover, ‘hydrolase activity acting on glycosyl bonds’ and ‘iron ion binding’ were the two dominant groups in the cascade of ‘catalytic’. Within the last sub-category ‘Biological Process’, ‘cellular process’ and ‘metabolic process’ were the two most prevalent cascades that can represent the typical activities of biological processes. Specifically, the most intriguing GO terms in ‘cellular process’ were found to be ‘meiosis I’ and ‘pollen wall assembly’, suggesting their active roles in MS. It was noted that ‘DNA recombination’ was highlighted in the cascade ‘metabolic process’.

在KEGG分析中，共富集了34条通路，其中13条来自于上调和下调基因，其余均来自于上调和下调基因(见表)2.）.研究表明，大多数基因参与“代谢途径”和“次生代谢产物的生物合成”。有趣的是，在不育芽中下调的基因列表中，至少有6个基因(SIN_1006103、SIN_1017099、SIN_1014074、SIN_1023392、SIN_1015497和SIN_1014349)注释为“减数分裂-酵母”或“卵母细胞减数分裂”，与GO注释结果一致。在“次生代谢产物的生物合成”途径中，上调基因的数量是下调基因的近3倍。此外，更多上调的基因被注释为“多环芳烃降解”和“α -亚麻酸代谢”。相比之下，不育芽中许多下调的基因在“抗坏血酸和醛达酸代谢”和“甘油磷脂代谢”中富集。此外，在类黄酮生物合成途径中有14个基因上调(见表)2.）.

通过使用MapMan对代谢途径进行更具体的比较，进一步支持了这些发现[32]．所有鉴定出的1502个不育芽与可育芽之间的DEGs均在TAIR数据库中进行了注释(http://www.arabidopsis.org）.结果发现，1445个DEGs与1240个拟南芥基因同源2.:表S3)。解剖1445个DEGS是可能与MS表型相关的假定功能，我们完全可视化与MapMan拟南芥同源基因和推断的候选途径网络（图3.）.

在该网络中，基因转录丰度的最显著变化与“蛋白质”、“靶向”、“激素”和“DNA”有关。此外，“乙烯和JA合成”相关基因的表达在不育芽中上调，而“信号通路”相关基因在DGMS不育芽中下调。此外，参与“脂质（FA合成）”、“氧化还原（抗坏血酸和谷胱甘肽）”和“能量（转运p-和v-ATP酶）”的DEG均下调，而在无菌芽中参与“第二代谢（类黄酮）”、“细胞壁（修饰）”和“能量（发酵）”的DEG上调，如果与那些处于可育芽中的相比。在“RNA-TF”组内差异表达的转录因子中，所有NAC、三螺旋和WRKY（除了一个WRKY）均上调，而C2C2（Zn）DOF、CCAAT和SET均下调。此外，在“信号”类别中，两个MAP激酶编码基因在不育芽中下调（图。3.；额外的文件2.:表S3)。

雄性不育/雄性生殖相关基因的鉴定

为了更深入地了解MS的分子机制，我们对PMRD中含有548个拟南芥雄性不育/雄性繁殖相关基因的芝麻DEGs进行了查询。检索到49个与植物雄性生殖有关的同源基因;其中一些基因编码转录因子(如bHLH089、MYB99和AMS)。转录因子编码基因在不育芽中表达减少，表明它们在确定绒毡层命运和发育中发挥重要作用3.）.

deg的等位变异

为了更好地了解DEGs，我们进一步利用SnpEff预测器预测等位变异对其靶基因功能的影响。在可育芽中表达的982个基因中共检测到1057个插入/缺失(InDels)，其中52个位于48个DEGs内(部分基因有两个InDels)3.：表S4）。类似地，在无菌芽中表达的1,354个基因中检测到1,432个诱导，86个吲哚美含量在83点内（附加档案4.:表S5)。我们共同鉴定了131个基因中的138个InDels，这些基因要么在可育芽中表达，要么在不育芽中表达。在已鉴定的138个InDels中，62个位于不育芽中上调的57个基因中，76个位于不育芽中下调的68个基因中5.:表S6和6.:表S7)。

其中，在上调基因列表中发现了SIN_1002610(乙烯响应转录因子ERF106)、SIN_1024026 (NAC2)、SIN_1019334 (WRKY 28)和SIN_1011023 (WRKY 33)等转录因子编码基因。一些编码“brassinosterid -regulated protein BRU1”(SIN_1022411)、“COP9 signalosome complex subunit 2”(SIN_1015172)和“Defensin J1-2”(SIN_1021298)的基因也被突出显示(附加文件)5.：表S6）。在下调的基因，SIN_1008339（E3泛素 - 蛋白连接酶MARCH1），SIN_1010740（L-抗坏血酸氧化酶同系物），SIN_1026145（花粉特异性蛋白SF3），SIN_1005014（蛋白质二硫键异构酶5-3）和SIN_1010051列表（糖转运蛋白8）有兴趣的，他们可能与花粉发育（附加文件有关6.: S7)。

我们选择21个含有预测导致功能丧失(LOF)和/或密码子改变(CC)的InDels的基因子集进行进一步分析(表)4.）.其中，在6个不育等位基因和6个可育等位基因中检测到可能导致CC的InDels(称为“CC型”)。此外，在6个可育等位基因和7个不育等位基因中也检测到lof型InDels，它们分别在可育芽和不育芽中表达量较高(用星号标记;表格4.）.因此，LOF-Type InDel可能导致其所在的转录本丰度增加。在可育芽中上调的11个基因中，大多数InDels(9 / 11)在可育等位基因中检测到，进一步证实了这一观察结果。同样，在不育芽中上调的其他10个基因中，大多数InDels(80%)在不育等位基因中检测到。

特别是SIN_1025190 (SCP18，丝氨酸羧肽酶)，SIN_1017245 (F3PH，类黄酮3'-单加氧酶)和SIN_1018350 (综合产品组在不育芽中，腺苷酸异戊烯基转移酶(Adenylate ispentenyltransferase)在不育等位基因中均表达上调。可育性等位基因中含有LOF-和CC-型InDels的激酶(SIN_1004626)基因在可育芽中表达上调(不育芽中表达下调)。有趣的是，在另一个基因SIN_1005818 (HMGB9, High mobility group B protein 9)中，在两个等位基因中均检测到InDel，在不育等位基因中可能存在干扰性inframe_deletion，在可育等位基因中可能存在LOF。该基因在不育芽中表达下调，在可育芽中表达上调4.额外的文件7.：表S8）。综上所述，在这些DEG中检测到大量序列变异，它们对转录物丰度的影响不是决定性的。

实时定量PCR验证

为了验证RNA-SEQ结果，我们为调节和下调的DEG选择了另一种策略。通过使用用于RNA-SEQ的相同的RNA样品，随机选择二十个基因进行验证以验证。底漆设计用于跨越外显子连接（附加文件8.:表S1)。结果表明，虽然qRT-PCR检测到的基因表达倍数变化在大多数情况下高于RNA-Seq，但两者的变化趋势相似，从而证实了RNA-Seq检测的准确性和可靠性。以随机选取的12个雄性不育/雄性生殖基因表达谱为例，列于表中5.，表明qRT-PCR与RNA-Seq的表达水平高度相关(r= 0.762,P< 0.01, n = 12)。

据我们所知，我们首次在转录组水平对芝麻DGMS进行了研究。利用Illumina测序平台，RNA-Seq获得了可育芽和不育芽的转录本丰度。然后绘制了高质量的转录组读到芝麻参考基因组和识别超过22数以千计表达基因,其中只有1502个基因(~ 6.6%)以不同的方式表达了无菌或肥沃的味蕾,表明有限数量的关键基因足以改变特征显著,花药的发育是一个复杂的多基因过程。

我们在芝麻中鉴定了49个与花药发育相关的基因，这些基因在拟南芥中具有同源性，其中一些基因编码转录因子(bHLH089、MYB99和AMS)，可能与绒毡层的命运和发育有关(见表)3.）.这些中，32个下调和17，其余均上调。此外，MS基因（克隆）的同源物占了将近一半每个控制类别中的基因，和基因的其余部分被注释为MR相关（雄性再现有关的基因，与GO证据），从而证明所有这些基因可能负责MS（表很好的候选人3.）.这可以解释为这里提到的芝麻MS是从PMC开始的，PMC是花药和花粉发育途径的第二阶段[13]因此导致口腔发育的失败，如在雄性无菌芽中所观察到的[8.]．

具体来说，我们发现了DYT1和TPD1是在32个下调的deg名单中(附表3.）.先前的研究表明DYT1可能通过表达调控花药发育自动对盘及成交系统和许多绒毡层优先的基因，从而间接地影响花粉壁形成[17]．TPD1是一种小肽，主要在小孢子细胞中表达，并可能分泌到绒毡层与雄性生殖细胞的界面中，与富含亮氨酸的重复受体样激酶EMS1相互作用形成受体复合物，从而决定绒毡层细胞的命运[16,33]．因此，很有可能下调的调控DYT1和TPD1芝麻可能通过决定绒毡层的细胞命运来影响花粉的释放。

另一个感兴趣的基因是RBOHE(RESPIRATORY BURST OXIDASE HOMOLOGUE E)RBOHE(At1g19230)是花药优先或绒毡层富集基因，功能缺失RBOHE导致绒毡层延迟变性，因此表达RBOHE减少在DYT1和tdf1［33]．与此相一致的是，我们发现RBOHE芝麻同源基因SIN_1024646和SIN_1007549在不育芽中的表达也显著降低(log_2.S/F =−1.7和−0.9)，如果与可育芽相比(附加文件1.:表S2)。因此,RBOHE在芝麻DGMS中可能具有类似的功能。

除了DYT1上面所提到的,QRT2(QUARTET2)也在MS基因(克隆)列表中(表)3.）.三个QRT基因包括QRT2当吉麦孢子从它们的四纹释放时，花粉母细胞壁的降解是所必需的[12]．此外,QRT2所需要的花药开裂。在花脱落的过程中由JA，乙烯和脱落酸（ABA），其共同被调节，QRT2由乙烯和ABA调控[34]．此外，与花药开裂相关的聚半乳糖醛酸酶活性可能受到JA、乙烯和ABA的调控[13]．在本研究中，表达减少QRT2与乙烯合成相关基因的上调有关。

拟南芥中有17个上调的芝麻基因，同源物（8与MS基因和9至Mr基因，表3.）.其中，SIN_1007695 (亚精胺羟基肉桂酰转移酶,SHT)，不育芽>增加200倍，使人联想到SHT在拟南芥花药绒毡层中表达的[35]．而且，SHT被FlowerNet的在线工具[分为“集群81”36]，包括KCS10、GH31和ATA7等几个基因;它们在芝麻中的同源基因SIN_1007525、SIN_1025709和SIN_1002500在不育芽中共同上调(附加文件)1.:表S2)，暗示它们可能参与了MS.这个‘cluster 81’也包含TSM1(绒毡层特异性METHYLTRANSFERASE1)，它编码一种阳离子依赖的ccoaomt样蛋白，参与苯丙多胺结合生物合成，并在拟南芥的雄蕊/花粉发育中发挥作用[37];具有未知功能的其余基因可能在花粉外来和脂质生物合成中发挥作用，基于他们在Atensembl中的描述[36]．因此，我们有必要对该基因簇内的其他基因进行研究，以明确它们的功能。

茉莉酸是花药发育过程中花药开裂所必需的[38]．参与JA生物合成和感知的基因突变导致花药开裂失败或延迟，花粉不能生存，从而导致雄性不育[39]．这类基因的例子包括花药开裂缺陷1(DAD1)，编码磷脂酶A1，催化JA生物合成的起始步骤;AOS，编码氧化烯合成酶的基因;裂开1(DDE1）/OPR3，编码JA合成途径中的OPR蛋白12-氧-植物二烯酸还原酶[40]．JA通路的所有阶段的缺陷似乎都导致了相似的表型，即灯丝伸长减少和缺乏开裂。在JA生物合成酶缺陷的突变体中观察到延迟开裂或不开裂表型[13]．在本研究中，SIN_1016850(同系物PLA15,磷脂酶A1-Igamma1)在不育芽中显著上调，而氧化丙烯合成酶编码基因的同源基因没有差异(数据未显示)。而SIN_1022877和SIN_1022878是OPR1(12-oxophytodienoate还原酶1）在拟南芥中，显示在无菌芽明显下调（附加文件1.:表S2)。这些数据强烈表明，参与JA通路的基因也与芝麻的MS有关。

植物基因表达调控是一个复杂的网络。转录因子通过与顺式作用的靶标元件的特定相互作用，调控一系列相关的下游靶标，这些靶标在植物的发育和对环境胁迫的响应中发挥着重要作用。拟南芥花药盲目因素(如果）一种类似NAC的基因，作为阻遏物，其通过调节茉莉酸盐生物合成中的基因来控制花药的裂缝[38]．事实上，在Swissprot注释的NACs中，9个芝麻同源基因在不育芽中的表达均上调，这加强了NACs在MS调控中的作用(图2)。3.额外的文件1.:表S2)。此外，在不育芽中显著上调的12个wrky中，有11个被注释为WRKY33的同源基因(图。3.额外的文件1.:表S2)。WRKY33蛋白是进化保守的，在广泛的植物胁迫反应中发挥着关键作用，而拟南芥WRKY33是激素和代谢反应的关键转录调控因子[41]．此外，参与氧化还原平衡的基因，水杨酸（SA）信号，乙烯 - JA介导的交叉通信和camalexin的生物合成被鉴定为WRKY33 [的直接目标42]．此外，ja相关反应的下调似乎涉及多个茉莉酸锌域基因(编码ja反应途径的抑制因子)的直接激活，包括WRKY33功能的丧失和额外的sa依赖的WRKY因子。在本研究中，NACs和wrky的共表达行为表明它们在调控芝麻MS中起着关键作用(图2)。3.额外的文件1.:表S2)。

为了了解序列变异对基因表达的影响，利用SnpEff预测等位变异对其靶基因功能的影响。有趣的是，在可育性等位基因中发现了6个InDels，这些InDels在可育芽中表达上调(而野生型不育等位基因在不育芽中表达水平较低);不育等位基因中有7个InDels在不育芽中表达上调(表1)4.）.这一观察结果表明，序列变异对转录本丰度的因果影响并非如此直接，而是令人困惑的。这可以解释为，大多数InDels是在编码区域而不是启动子区域检测到的，在启动子区域，它可以直接影响转录本的丰度。偶尔，我们也发现了具有转录调节功能的InDels，其中，致病InDels的存在降低了转录丰度。例如，不育等位基因中含有InDels的两个基因(SIN_1025700和SIN_1005818)导致不育芽转录丰度降低，而可育等位基因中含有InDels的两个基因(SIN_1004703和SIN_1019529)导致可育芽中基因表达下调，表现为顺式作用。

正如Rutley和Twell所建议的[43]，对雄配子体的转录组研究不仅增加了我们的知识和理解，而且通过告知实验设计（如通过反向遗传学的基因选择）以及基于查询和共表达分析，提高了实验策略的有效性。本研究提供了许多DEG和一些候选基因，可用于通过转基因验证阐明芝麻DGMS的分子机制。

本研究提供了一套基于转录组分析的芝麻DGMS株系可育芽和不育芽中差异表达的1502个基因。这些基因中有一半在不育芽中上调，表现出复杂的表达模式。乙烯和JA合成及信号转导相关基因在不育芽中分别表达上调和下调。此外，大部分NAC和WRKY转录因子在不育芽中上调。

此外，49个基因芝麻与雄性不育/男性再现相关的拟南芥同源物在无菌花显示出降低的表达。这些基因编码转录因子（bHLH089，MYB99，和AMS）的是可能在指定或确定绒毡层的命运和发展中的作用。此外，等位基因变体对靶基因的功能预测效果强调了几个个InDel，这可能导致生育能力的决心。

植物材料及RNA制备

本研究使用的芝麻植物材料包括新开发的DGMS株系W1098A和其可育性对应的W1098B，两者仅在花粉育性方面存在差异[8.]．这两个品系均在中国农业科学院油料作物研究所(中国湖北省武汉市)试验田栽培。从5株雄性不育株和可育株中分别剥离出2.5 mm长的芽并进行膨大转录组分析。育芽体和无菌芽体立即在液氮中快速冷冻，然后在−80°C冰箱中保存直到使用。用TRIzol试剂(Gibco-BRL)从无菌芽和可育芽中按说明书提取总RNA。然后从不育芽和可育芽中构建两个cDNA文库，如前面在芝麻[44]．简而言之，将大约5mg mRNA分段，转化为cDNA，并根据Illumina RNA-SEQ协议（Illumina，Inc.San Diego，CA）扩增PCR。使用Illumina Genome Analyzerii（Sandiego，CA）和Illumina Hiseq 2000平台（San Diego，CA）在北京基因组学院（中国沉珍）生成序列读数。

差异表达基因的鉴定

这些干净的reads被映射到参考基因组序列美国indicum(http://ocri-genomics.org/Sinbase/) [11]使用SOAP对准/ SOAP2（简称读对准的改进的超高速工具）45]．RPKM用于测定每个基因的相对转录本丰度。利用DEGseq程序，在可育芽库和不育芽库之间鉴定出显著差异的基因表达[46]．FDR用于确定阈值p-在这项研究中，一个严格的FDR ≤ 0.001和│log2（不育/可育的倍数变化率）│ ≥ 1.00作为阈值，选择一个显著不同的表达基因。

遗传变异特征

使用SnpEff 4.1版本对序列变异(如InDels)进行了表征[47]，参考从Sinbase下载的芝麻基因组注释(http://ocri-genomics.org/sinbase.)根据Wang等人[48]．根据Saeed等人描述的方法，预测了可能对转录本/蛋白质有高度影响的序列变异(InDels，帧移位，停止获得，停止丢失和非同义编码)[49]．

GO和KEGG通路富集分析

DEGs用于氧化石墨烯和途径富集分析。一个纠正P选取≤0.05为显著性阈值，确定基因集富集[50]．根据集成数据库提供的GO映射来分配从deg中推断出来的函数类。Blast2GO计划(https://www.blast2go.com/)用于获取所有deg的GO注释[51]。然后，结果提交给WEGO(http://wego.genomics.org.cn)生成所有deg的GO分类图[52]．

KEGG通路分析是基于我们的定位基因与当前KEGG数据库的比较结果[53]．MapMan (version 3.5.1 R2)也用于将DEGs注释到代谢途径上。

通过qRT-PCR确认候选deg

为了验证由RNA-SEQ检测DEGS，20度的视角，随机从52共同差异表达的基因在两个文库中选择，然后进行qRT-PCR分析，根据Qi等。［54]．利用在线工具Primer3设计基因特异性引物[55]，基于选定的unigenes序列(附加文件8.:表S1)。在Bio-Rad CFX96仪器中使用SYBR Green Real time PCR Master Mix (TOYOBO, Japan)进行反应。对于每个样本，每个基因在96孔板上重复3次。各基因的相对表达量由2^-ΔΔC_T方法(56]．所有数据均用均数±标准差表示。

阿坝:

乙烯和脱落酸

如果:

花药盲目因素

AMS：

夭折的小孢子(AMS)

AOS:

丙二烯氧化氮合酶

CMS：

细胞质雄性不育

DAD1:

花药开裂缺陷1

DDE1:

裂开1（也称为OPR3)

度:

差异表达基因

DGMS：

显性基因雄性不育系

DYT1:

功能失调的TAPETUM1

EMS1:

过量的小孢子细胞1(也称为体外产孢细胞)

罗斯福:

假发现率

罗斯福:

假发现率

总经理:

细胞核雄性不育

走:

基因本体论

走:

基因啊ntology.

HMGB9:

高迁移率B组蛋白9

InDel:

插入/缺失

IPT：

腺苷酸isopentenyltransferase

JA：

茉莉酸

ke:

京都enc基因和基因组百科全书

格林:

连锁群

格林:

连杆组

LOF:

损失函数

女士:

雄性不育

女士:

雄性不育

MS1:

男性STERILITY1

OPR1:

12-oxophytodienoate还原酶1

纤毛运动:

程序性细胞死亡

PLA15:

磷脂酶A1-Igamma1

PMC:

花粉母细胞

PMRD:

植物男性生殖数据库

QRT:

四方

QRT2:

QUARTET2

存在:

实时定量PCR

RBOH：

呼吸爆发氧化酶同系物

RBOHE:

呼吸爆发氧化酶同源物e

RPKM:

读取每千碱基外显子模型每百万映射读取

RPKM:

读取每千碱基外显子模型每百万映射读取

SA：

水杨酸

SHT:

亚精胺hydroxycinnamoyl转移酶

SHT:

亚精胺hydroxycinnamoyl转移酶

TDF1:

绒毡层发育和功能有缺陷

TDF1:

绒毡层的发育和功能1

TPD1:

绒毡DETERMINANT1

我们感谢华大基因(中国北京基因组研究所)的同事就RNA-seq取样进行了有价值的讨论，并帮助解释转录组数据。

资金

农业部油料作物生物学与遗传改良重点实验室开放基金项目(2016003)资助感谢国家自然科学基金(31101180)和国家农业科研体系(CARS-15)的资助。

数据和材料的可用性

本研究使用的原始RNA-Seq数据已保存在美国国家生物技术信息中心(NCBI)序列读取档案(SRA)数据库中，登录号为SRP076254 (http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP076254）.

作者的贡献

HL、MT和YZ设计了本研究;HL和YZ准备了用于测序的植物材料。LL和MT对NGS数据进行生物信息学分析;HY、FZ、MY、TZ进行qPCR实验和统计分析;HL和MT解释数据并撰写手稿。所有作者阅读并批准了最终的手稿。

同意出版

不适用。

竞争利益

两位作者宣称他们没有相互竞争的利益。

伦理批准和同意参与

不适用。

作者笔记

隶属关系

1. 石油作物，农作物，中国农业科学院石油作物，湖北武汉石油作物研究所的生物学和遗传改善重点实验室，湖北武汉430062

   刘红艳，周方，杨敏敏，周婷，赵英忠

2. 南京农业大学生命科学学院，江苏南京210095

   谭明普，于海娟，李亮

相应的作者

对应于映赵。

开放获取本文遵循知识共享署名4.0国际许可协议(http://creativecommons.org/licenses/by/4.0/）如果您向原始作者和源给出适当的信用，则允许在任何介质中进行不受限制的使用，分发和再现，提供指向Creative Commons许可证的链接，并指示是否进行了更改。Creative Commons公共领域奉献豁免（http://creativecommons.org/publicdomain/zero/1.0/)适用于本条提供的数据，除非另有说明。

再版和权限