比较转录组分析揭示花园芦笋(芦笋officinalis）

背景

花园芦笋(芦笋officinalis)是一种非常有价值的蔬菜作物，具有商业和营养价值。它也常用于研究植物性别决定和分化的机制。然而，对芦笋的性别表达机制仍知之甚少。

结果

通过Illumina配对端测序进行从头转录组测序，揭示了来自雄性和雌性芦笋花蕾的超过260亿个高质量序列数据。共组装了72,626个unigenes，平均长度为979 bp。通过比较转录组分析，鉴定出4876个差异表达基因(deg)存在于雌性和雄性/超雄性花蕾可能的性别决定阶段。其中433个基因被标记为花相关基因，其中285个为雄性/超雄性偏倚基因，149个为雌性偏倚基因。在雄性/超雄性偏向的花相关基因中，102个可能与花药发育有关。此外，还发现了43个与激素反应和生物合成有关的DEGs，这些DEGs被认为与性表达和生殖有关。此外，还发现了来自不同家族的128个转录因子(TF)相关基因的差异表达，这一发现暗示了TF在芦笋性别决定或分化中的重要作用。相关分析表明，miRNA-DEG对也与芦笋的性发育有关。

结论

我们的研究发现了大量参与芦笋性别表达和生殖的deg，包括参与植物生殖的已知基因、植物激素信号、TF编码基因和功能不明确的基因。我们还发现mirna可能参与了性别分化过程。本研究为进一步研究芦笋性别决定和分化的调控网络提供了有价值的基础，并为进一步开展这一雌雄异株物种的遗传和基因组研究提供了依据。

雌雄异株系统是动物中最常见的规则，但大多数开花植物是共性的，例如，单株植物由具有两性功能的两性花组成。大约6%的被子植物是雌雄异株植物，在不同的个体上有单性花，即雄花在雄蕊上，雌蕊花在雌蕊上[1]。植物的性别决定机制多种多样，从XY性染色体到完全的常染色体决定系统[2]。此外，植物的雌雄异株已经独立进化了很多次。然而，雌雄异株植物通常来自于雌雄同体植物的谱系。不同性别的决定和分化是由一对携带决定性别基因的性染色体在基因上控制的[3.]。至少有两个显性基因，即一个促进雄蕊的基因和一个抑制雌蕊的基因，在性染色体中相连[4]。尽管近年来人们对性别决定调控机制进行了大量的研究，但目前国内只报道了一种在雌雄异株植物中被证实的性别决定基因Diospyros莲花(5]。在各种雌雄异株植物中，甚至在基因组测序的物种中，负责性别决定的机制和基因仍不清楚，主要是因为性别决定基因所在的Y/W染色体的非重组区域通常被重复序列占据[6]。

与动物的生殖细胞系在早期发育阶段分化不同，植物，包括雌雄异株物种，不具有独特的生殖细胞系[7]。在植物中，全能分生组织细胞通常经历一个漫长的营养期，然后进入生殖阶段形成花，这是复杂的性器官[8]。性别特异性表型，特别是不同个体的雄花和雌花的发育和维持是在一系列的代谢途径和调控遗传网络下进行的，其中涉及各种基因、转录因子(TFs)和其他调节因子，如microRNAs (miRNAs) [qh]9，10，11]。与哺乳动物相似[12]，植物性别分化所必需的下游代谢途径和遗传网络可能受到上游性别决定基因的控制[5]。然而，植物性别决定和分化的机制尚不清楚。

花园芦笋(芦笋officinalis百合科植物，是世界上重要的雌雄异株蔬菜作物。它包含一个单倍体基因组大小为1308 Mb和2n= 2×= 20条染色体[13]。由于芦笋具有严格的基因型性别决定方式、稳定的性别表型和可行的YY基因型个体(超级雄性)，因此芦笋也是性别决定机制研究的好物种[14]。芦笋的两性二态性是由卵黄决定的米在同态性染色体L5上的位点。特别是,毫米，毫米,毫米分别代表超级男性、男性和女性[14，15]。虽然各种分子标记与米已确定基因座[16，17，18),米基因还有待分离和鉴定。最近，Murase等人发现了一个男性特有的基因MSE1，在芦笋的性别决定中起重要作用[19]。然而，这个基因是否是确切的性别决定基因米应进一步调查。各种数据集和大量分析，包括EST、转录组和miRNA，为研究涉及性别决定和分化的复杂过程提供了有用的资源[6，20.，21]，但与芦笋性别决定和分化有关的基因仍不清楚。

随着下一代测序技术的进步，转录组测序已成为基因发现和相关途径富集的宝贵工具。这项技术已被应用于鉴定植物中潜在特征的候选基因[22，23]。更重要的是，已经进行了几项研究来检测不同性别类型的deg，例如在柳树suchowensis(24]，番木瓜(25]，Pelteobagrus fulvidraco(26等等。最近，Harkess等(2015)对雄性和雌性芦笋矛尖进行了转录组测序，发现了许多参与雄性和雌性配子体发育的基因。虽然他们提供了关于花药和雌蕊发育途径的有用信息，但他们分析了芦笋矛尖的数据，其中包括营养器官和生殖器官。他们还提出，研究花蕾和花蕾原基中的性别偏见表达比研究矛尖中的这种现象更有意义[6]。在我们的研究中，我们使用Illumina Highseq 2500测序平台对雄性和雌性芦笋的花蕾进行了比较高通量的转录组测序，以确定它们对应的deg。然而，超级男性系列并没有出现。为了更好地阐明雌雄芦笋的性别决定和分化机制，我们利用芦笋的超级雌雄花蕾将我们的测序数据与先前发表的数据联系起来[6]。我们还进行了miRNA-DEG相关分析，以确定芦笋中参与性别决定或性别分化的关键调控mirna靶向基因。

芦笋花蕾的测序和转录组组装

我们对雄性和雌性花蕾样本进行Illumina测序，每个样本分两个重复制备。从芦笋雌雄花蕾中共获得210,849,964个原始测序reads。在除去适配器序列、不明确的核苷酸、低质量序列和所有可能的污染后，还剩下205,739,016个干净reads。其中，107,013,940个来自雌性花蕾，98,725,076个来自雄性花蕾1)。为了提高组装的质量，我们将我们的测序数据与之前报道的雌性和超级雄性花蕾的数据结合起来[6]进行序列组装。这些高质量reads的组装得到72,626个unigenes, N50长度为1569 bp，平均长度为979 bp。最短和最长的单基因分别为201和15,615 bp。此外，34,995个(48.2%)unigenes大于500 bp, 19,967个(27.5%)unigenes大于1000 bp1)。BUSCO评估表明，83.2%的植物同源物被组装的单一基因所覆盖(附加文件2)。Transrate评估的程序集得分为0.242，高于Transrate得分0.22的阈值。这些结果表明，我们获得了一个高质量的芦笋花蕾转录组。

芦笋转录组的功能注释与分类

在4个公共数据库(Nr、Swiss-Prot、COG和KEGG)中进行BLAST检索，对组装的unique进行注释。结果显示，这些数据库中的单基因注释分别约为40.96%、29.39%、24.75%和16.07%(表1)2)。总共有29967个(41.26%)unigenes在至少一个公共蛋白数据库中被成功注释(表1)2、附加文件3.)。对于COG注释，17973个unigenes被分配到COG分类，并分为25个类别(图2)。1)。“仅一般功能预测”(6564,36.52%)是最大的群体，其次是“翻译后修饰、蛋白质周转、伴侣”(3344,18.61%)、“信号转导机制”(2727,15.17%)和“RNA加工和修饰”(1725,9.60%)。很少的unigenes被分配到“细胞运动”(12.0.07%)。

物种分布分析表明，芦笋的独特基因与各种植物的相同(图2)。1 b)。在众多的植物物种中，Elaeis guineensis(6736)，然后是凤凰dactylifera(5515),穆萨acuminata(2695)和栽培稻(1853)，包含了与芦笋的独特基因匹配度最高的序列。因此，大部分被注释的单基因序列与单子叶植物的序列高度吻合，这与芦笋是单子叶植物的事实是一致的。

芦笋雌雄花蕾转录组的比较

数字表达谱分析显示，22,335个基因在雄性和雌性花蕾中差异表达(图2)。2)，而18550个基因在超雄性和雌性花蕾中差异表达(图2)。2 b)。一个负责性别决定或性别分化的特定基因在雄性和超雄性植物中表现出与雌性植物相似的表达模式。因此，我们确定了雄性/超雄性和雌性植物之间具有相似表达模式的基因，以便进一步分析。雄/超雄与雌花蕾之间共有4876个基因差异表达(图2)。2摄氏度)。在这些基因中，2401在雄性/超雄性花蕾中的表达明显高于雌性花蕾，2475在雌性花蕾中表现出偏表达(图2)。二维、附加文件4)。

在这些deg中，通过检测性别特异性基因来寻找在不同性别类型中存在或缺失的潜在基因。使用以下标准:deg在一种性别类型中含有绝对零计数的估计丰度，但在另一种性别类型中显示一定的表达(每百万碱基读取量，RPKM >0.03)。最终鉴定出166个性别特异性基因。其中120个为雄性/超雄性特异性，46个为雌性特异性5)。因此，大多数基因表现出雄性/超雄性特异性表达模式。

为了证实RNA-Seq数据的可靠性，我们通过实时定量PCR检测了半随机选择的24个单基因的转录水平(图2)。3.)。我们随机选择了18个男性偏倚基因和6个女性偏倚基因，因为男性偏倚基因被认为比女性偏倚基因更具信息性。虽然RT-PCR与RNAseq的fold-change值不同，但其unigenes的表达模式与RNA-Seq数据一致(R²= 0.7944)。

然后将雄性/超雄性和雌性花蕾之间的deg作为输入进行基因本体(GO)注释，并将GO术语分为三类:细胞成分、分子功能和生物过程。共有1619个、1060个和2193个雄性/超雄性偏倚基因被分配到这些类别，而1118个、642个和1402个雌性偏倚基因被分配到这些类别6)。

为了更好地注释deg，我们收集了所有拟南芥与植物繁殖和花发育有关的基因。BLAST分析显示5809个单基因推定同源拟南芥基因产品。在这些单基因中，有433个是雄性/超级雄性和雌性花蕾之间的基因变异，其中包括284个雄性/超级雄性偏倚基因和149个雌性偏倚基因。4、附加文件7)。在标记为生殖和花发育相关的序列中，102个单基因被有效地标记为花药相关基因或雄性不育突变相关基因(图2)。4)。芦笋花蕾雄性/超雄性偏表达基因包括同源的夭折的小孢子s (自动对盘及成交系统)，聚半乳糖醛酸酶2（PG2)，小孢子特异性启动子2（MSP2)，以及其他花药相关基因。其中，11个表达偏雄/超雄的基因参与了花粉壁形成过程，如花粉外壁形成、孢粉素生物合成、花粉亚精胺形成和花粉管形成/生长。一些基因，包括同源盒基因（ATH1)、F box基因、植物激素相关基因和编码tf的基因也被鉴定出来，这些基因的表达具有雌性偏向性，与雌性配子体发育有关。

芦笋花蕾KEGG通路分析

对雄/超雄和雌花蕾之间的deg进行KEGG分析，发现300个unigenes被分配到KEGG生化途径的前18位p价值(无花果。5)。这些途径可分为五类，即“氨基酸代谢”、“其他次生代谢物的生物合成”、“碳水化合物代谢”、“能量代谢”和“脂质代谢”。苯丙素和类黄酮生物合成途径中的许多基因在雄性/超雄性和雌性花蕾中表达差异，这一观察结果与先前在鱼尖上的发现一致[6]。基于通路的分析进一步表明，鉴定出43个deg，并与“植物激素信号转导”通路有关。在这些基因中，28个表现为男性/超男性偏倚表达，15个表现为女性偏倚表达。这些基因包括加勒比海盗，蚂蚁，阿富汗二月,小块土地，参与调节多种激素稳态和生殖过程。富集分析表明，大多数这些deg参与生长素、乙烯、细胞分裂素和JA信号通路(图2)。6、附加文件8)。在生长素信号传导中，鉴定了8个雄性/超雄性偏倚转录本和10个雌性偏倚转录本，编码生长素合成酶、生长素诱导蛋白、生长素响应蛋白和生长素转运蛋白。此外，还检测到1个编码生长素内流转运蛋白的雄性/超雄性偏倚转录本和2个编码生长素抑制蛋白和生长素外排蛋白的雌性偏倚转录本。在乙烯信号通路中，3个雄性/超雄性偏倚的单基因和1个雌性偏倚的单基因编码乙烯应答转录因子(ERF)，另一个雌性偏倚的单基因为ap2样蚂蚁基因。在细胞分裂素信号通路的4个雄性/超雄性偏倚转录本中，3个是加勒比海盗基因和1是编码含组氨酸的磷酸转移蛋白的基因。此外，在雌性花蕾中，负责细胞周期蛋白d3 - 2样蛋白的1个单基因被下调。在茉莉酸(jasmonic acid, JA)信号通路中，3个tify样基因和1个myc4样基因在雄/超雄花蕾中表达上调。还检测到一些与其他植物激素信号通路相关的deg。

TFs在雌雄花蕾中的差异表达

共有128个TF基因在雄性/超雄性和雌性花蕾中差异表达(附加文件)9)。其中70个基因在雄性/超级雄性中高表达，58个基因在雌性中优先表达。热图显示了TF基因在雄性/超雄性和雌性花蕾之间的差异表达谱。7。差异表达的TF包括C2H2、C3H、bZIP、ERF、bHLH、MADS、MYB、WRKY、NAC等小TF家族成员。各家族中以bHLH和MYB家族居多。值得注意的是，17个bHLH基因在雌雄花蕾中差异表达，其中5个在雄/超雄中偏表达，12个在雌/超雄中偏表达。MYB/MYB相关TF家族成员也表现出明显的差异表达模式，12个在雄性/超级雄性中高表达，5个在雌性中优先表达。其他TF家族成员，包括C2H2、LBD、C3H、DRF、NAC和WRKY，也具有男性/超男性偏倚的表达特征。在女性偏倚基因中，C2H2家族基因最多，其次是MADS和TALE家族基因。

miRNA-mRNA相关性分析

从雄性和雌性花蕾的混合物中构建了两个小RNA文库。从雄性和雌性花蕾中分别产生了6,837,349和6,571,668个原始reads。裁剪适配器后，去除低质量序列，过滤rnas、tRNAs、anRNA和snoRNAs，剩余的2,949,390和2,902,544个清洁reads，分别代表1,832,639和1,777,464个独特reads，分别用于鉴定雄性和雌性花蕾中的miRNA。在这些序列中，有72个已知的保守mirna，属于53个家族，在两个样本中鉴定出81个新的mirna(附加文件)10)。

差异表达miRNA (differential expression miRNA, DEM)分析显示，雄性花蕾中有15个DEM上调(已知10个，新发现5个)，雌性花蕾中有32个DEM上调(已知12个，新发现20个)。结合DEM分析结果和RNA-seq数据，根据miRNA与其靶标之间的负调控关系，从转录组数据中检测出DEM- deg靶对。共有90个靶点相对于30个相应的mirna表现出相反的表达模式(附加文件10)。在雄性花芽中，12个目标deg上调，57个目标deg下调，包括与花发育有关的unigene0002102 (squamosa promoter-binding like protein)和unigene0031841 (scareclike protein)。在雌性花蕾中，共鉴定出18个上调的dem和33个下调的靶deg，如unigene0034398 (MADS-box转录因子AP3)和unigene0001585(转录因子gamyb样蛋白)。

在本研究中，我们从转录组水平分析了芦笋的花蕾，以加深我们对芦笋性别决定和分化机制的理解。因此，我们对从芦笋雌雄花蕾中收集的cDNA进行了测序，这可能代表了芦笋性别决定和分化的关键阶段。芦笋包含两种类型的雄性个体:普通XY型雄性和YY型超级雄性。这两种类型的雄花发育过程相似。我们假设调节雄性和超级雄性性别决定和分化的主要因素是相同或相似的。因此，我们将测序项目中获得的数据与先前报道的减数分裂前雌性和超雄性花蕾的发现相结合，以确定雄性/超雄性和雌性花蕾之间的deg。这种组合和比较可以解决我们测序项目中仅使用两个重复的局限性。RT-qPCR验证了基于比较筛选的deg在不同个体和不同实验条件下的准确性和可靠性。这些男性/超男性偏倚和女性偏倚的基因对于阐明调节性别决定和分化的潜在分子变异是必要的。

雄性/超雄性与雌性花芽表达模式

我们观察到许多基因在雄性/超雄性花蕾或雌性花蕾中过表达或特异性表达。这些基因的功能注释揭示了花器官身份的维持、雄蕊和雌蕊的发育以及配子体发育等与花发育有关的生物过程的分子机制。

在FlowerNet数据库的基础上，我们将284个同源物标注为花相关基因。在这些同源物中，102个被鉴定为花药相关基因或雄性不育突变相关基因。一些与花药有关的基因，如galacturonase（PG),自动对盘及成交系统，是根据先前的一项研究[6]。这一发现表明，这些基因在花发育的不同阶段对花药发育至关重要。花药相关聚类分析进一步表明，有11个基因与花粉壁形成相关。阳离子/H+交换剂（CHX19)，DUF064，SULTR，SYTB,TSPO在花粉母细胞减数分裂中表达上调，并参与花粉外壁形成[27]。酰基辅酶:一种合成酶（ACOS5),细胞色素P450（CYP704B1)与孢粉素的生物合成有关。淘汰ACOS5和CYP704B1不要因为花粉壁受损或存在不规则的外胚层而产生花粉[28，29]。CYP86C3和花药7（A7)是雄蕊专一的拟南芥并且可能负责花粉亚精胺的形成。我们的研究结果进一步显示了两个雄蕊特异性基因，即，拟南芥LAT59同源（AT59),糖转运体（STP11)，在花粉管的形成和生长中起作用[27]。总的来说，这些表达偏向雄性或超雄性的芦笋同源物可能参与了花粉的形成和发育。然而，负责花粉发育的基因参与下游途径，不应该涉及性别决定。

花发育的同源物共有149个拟南芥在雌性中比在雄性或超级雄性中表达得更明显。这些基因与雌蕊和配子体发育有关。在不同性别类型的差异表达基因中，大量未知基因未产生显著的BLAST命中。这些基因可能包括那些芦笋特有的基因，那些在其他植物中功能不明确的基因，那些与其他已知基因序列差异很大的基因，以及其他基因。随着大量基因功能的阐明，可能参与芦笋性别决定或性别分化的基因逐渐被了解。因此，与花发育相关的deg应该进一步研究。

植物激素信号在芦笋性别分化中的作用

植物激素是具有多种信号功能的内源性调节剂，几乎影响植物生长发育的各个方面[30.]。多种植物激素可能参与调节单性花的性别决定基因和发育途径[11]。在雌雄同株植物中，如黄瓜(Cucumis巨大成功)和甜瓜(Cucumis梅洛)，单性花的发育可以通过影响性别表型的植物激素来调节。例如，乙烯抑制雄花发育[31]，而赤霉素则促进雌花的雄蕊发育[32]。在黄瓜和甜瓜中，经典鉴定米(生育)和F(雌性化)基因座是调节激素信号和促进器官流产过程的因素[33，34]。雌雄异株物种的性发育可能受遗传变异控制，对植物激素的敏感性较弱。然而，激素可能刺激从雌雄同株植物进化而来的雌雄异株植物的性别分化[11]。在麻(大麻)，几种植物激素的应用可以影响它们的性表型。例如，GA诱导雄性化，而IAA、乙烯和激动素对大麻的性别产生雌性化效应[35]。因此，应该检查激素调节途径，以获得鉴定性别分化和性别决定基因的有用信息。在我们的研究中，几种激素的信号通路，包括生长素、乙烯、细胞分裂素和JA，通过通路分析富集(图2)。5)。生长素途径可以与其他途径相互作用，在花卉发育中发挥重要作用[j]。36]。许多与生长素代谢和信号传导相关的基因参与了生长素介导的开花调控过程[36]。鉴定了芦笋雄性/超雄性和雌性花蕾中生长素诱导蛋白、生长素负责蛋白和生长素运输蛋白的独特转录本，并对其进行了差异调控。结果表明，一个生长素内流转运基因在雄性/超级雄性花蕾中优先表达，而一个生长素抑制蛋白编码基因和一个生长素外排蛋白编码基因在雌性花蕾中比在雄性/超级雄性花蕾中更多地表达。这些结果表明，雄性/超级雄性芦笋花蕾中的生长素水平可能高于雌性芦笋花蕾，并且在雄性/超级雄性芦笋花蕾中生长素的作用可能比雌性芦笋花蕾中更重要。这一观察结果与先前的结果相反，该结果表明生长素增加了雌雄同株植物中雌花的百分比[37]。然而，产生高浓度游离生长素的花器官会抑制或延缓邻近器官的发育[38]。我们预测，雄性/超雄性花蕾中高水平的生长素抑制了雌蕊等邻近器官的发育，促进了单性花的形成。事实上，在不同的研究中，生长素在植物性别调节中的作用和可能的作用显著不同[39，40，41]。因此，生长素的功能及其对性发育的影响有待进一步研究。其他基因，如小块土地，加勒比海盗,TIFY，也被发现与乙烯、细胞分裂素和JA信号有关。综上所述，这些结果表明植物激素介导的转录调控可能对芦笋的性别分化至关重要。对植物激素产生和信号传导关键基因的鉴定和进一步分析，将极大地促进对芦笋性别分化复杂遗传网络的研究。

芦笋性别决定与分化的TF基因

TFs在植物生殖发育中起着重要作用。许多基因，如AP1，AP2，东盟地区论坛,APL /出处，它们参与生殖发育，编码一个TF [10，42]。的确,对不起是一种众所周知的人类性别决定因素，它编码一个TF [43]。动物中一些与性别决定和性别分化相关的基因也编码tf [44，45]。最重要的是，唯一被鉴定和确认的性别决定基因编码HP-Zip型TF [5]。越来越多的研究表明，tf在植物和动物的性别决定和分化中起着重要作用。在我们的研究中，许多TF基因在雄性/超雄性和雌性芦笋花蕾中表现出显著的差异表达。在各TF家族中，锌指、MYB、NAC、bZIP和bHLH显著富集，表明它们在芦笋性别决定或性别分化中起着重要作用。这些科的成员已被证明参与其他雌雄异株或雌雄同株植物的性别分化[5，33]。在柿子中，同源结构域TF基因政府鼓励属于HD-Zip家族，以剂量依赖的模式调节花药生育力。奥吉是一种男性特有的基因编码一种小RNA政府鼓励，以及两者的相互作用奥吉和政府鼓励控制柿子的性别[5]。在我们的研究中，HD-ZIP家族的两个TF基因(Unigene0023095和Unigene0041406)在雌花芽中表达上调。先前的一项研究表明，HD-Zip转录本在雄性和雌性芦笋矛尖中不表现出性别偏见的表达[6这些观察结果表明，这两个HD-Zip基因可能在花蕾中而不是在花原基中启动偏转录。尽管雌雄异株在不同的开花植物谱系中独立进化了几次，但对这两个TF基因的进一步表达、遗传和功能研究应该为理解芦笋的性别决定或分化提供有用的信息。

锌指TF家族由大量含有一个或多个锌指结构域的蛋白质组成，并进一步分为不同的亚家族，如C2H2型、C3H型和C2HC型[46，47]。属于锌指家族的蛋白质参与各种调节过程[48，49]。在甜瓜,CmWIP1该基因编码C2H2锌指TF，通过间接抑制乙烯刺激促进雄蕊发育CmACS-7基因。两者之间的相互作用CmWIP1和CmACS-7控制瓜中雄花、雌花和雌雄同体花的发育[33]。在我们的研究中，6个c2h2型和5个C3H型TF基因表现出男性偏倚的表达模式，而5个c2h2型和1个C3H型TF基因表现出女性偏倚的表达模式。值得注意的是，一个清晰的SUPERMAN-LIKE同源基因在雌花芽中的表达量是雄花芽的4倍。在拟南芥,超人基因在细胞增殖过程中具有负向活性，在花轮和胚珠发育中起着重要的作用。50]，雄蕊和不育胚珠数量增加吃晚饭突变体(51]。在模式雌雄异株植物中硅宾latifolia,超人同族体SlSUP只在雌花中表达，与雌花发育有关[52]。

bHLH家族成员在营养生长、生殖生长、植物生长、形态发生和胁迫响应等方面发挥着多种作用[53，54]。在芦笋雌雄花蕾中bHLH基因的差异表达中，有5个基因在雄性中显著表达。其中一个基因是自动对盘及成交系统雄性生殖发育所需的基因，在绒毡层中特异表达拟南芥花药，参与绒毡层组织发育[55]。自动对盘及成交系统也可以调节其他几种花药发育相关基因的表达，如bHLH89和bHLH91(55]。我们还发现bHLH89在雄性或超级雄性中优先表达。因此，我们提出的相互作用bHLH89和自动对盘及成交系统与芦笋的花药发育有关。此外，包括光敏色素相互作用因子4(光敏色素相互作用因子4)同源物在内的11个bHLH基因(PIF4)，GLABRA3，bHLH57，bHLH66，bHLH74,bHLH75，参与开花时间调节、激素调节、应激反应和细胞分化，在雌花芽中显著且高表达[56，57，58]。

MYB/MYB相关TF基因家族参与多种重要的植物生物学过程，包括花的发育[59]。MYB/MYB相关转录本分别在雄/超雄和雌花蕾中优先表达12个和5个。在这些转录本中，有两个是同源的瑞士/ SNF它编码了染色质重塑复合体的蛋白质，这对叶和花的发育至关重要。这种TF可以调节多种花器官特征基因的表达模式，如APETALA1，APETALA3,AGL24(60]。最近，Murase等人简要报道了MYB TF基因MSE1这可能与芦笋的性别决定有关。MSE1在雄性芦笋中特异表达。击倒MSE1相同器官中拟南芥导致男性不育[19]。然而，这只是一个初步的研究，该基因的功能验证和潜在机制的研究还需要进一步的研究。

一些其他家族的TF基因，如MADS、WRKY、GATA、NAC、LBD、AP和SAP，在雄性/超雄性或雌性花蕾中也表现出偏表达。因此，这些基因也在芦笋的性别决定或分化中起作用。例如，aptala2 (AP2)样乙烯反应TF基因AINTEGUMENTA（蚂蚁)与男性或超级男性相比，更倾向于在女性中表达。蚂蚁参与女性发育的许多方面，包括胚珠和雌蕊的发育[61，62]。在其他植物中也有一些功能不明确的差异表达tf。因此，这些tf被认为是可能参与雄花和雌花发育的潜在候选者。这些基因的遗传和功能将为进一步阐明这种雌雄异株物种的性别决定和分化提供相关信息。

芦笋中参与性调节的mirna

mirna是关键的基因表达调控因子，在多种细胞过程中发挥重要作用。mirna可能参与动物繁殖和性别分化[63，64]和植物的发育过程，包括花的发育[65，66，67]。最近，Chen等人通过小RNA测序和降解分析发现了一些可能参与芦笋生殖器官发育的mirna及其靶点[21]。本研究对小rna进行测序，鉴定miRNA，并进行miRNA与转录组的相关性分析，以确定miRNA在调节芦笋性别分化中的潜在作用。

多种dem，包括miR156、miR159、miR165、miR172和miR319，据报道与花卉发育有关[68，69]，在雄性和雌性花蕾之间被识别出来。对dem的靶点进行了预测，并鉴定了多个DEG-DEM对，如mir159 - gamyb样蛋白、mir156 -鳞状启动子结合样蛋白、mir171 -稻草人样蛋白和miR5719-MADS-box转录因子AP3。

拟南芥GAMYB比如基因，比如MYB33和MYB65，受miR159调控，对花药发育至关重要[70]。突变体分析表明，花药发育myb33myb65双突变体在减数分裂前期被阻断，并被认为是男性不育的原因。虽然miR159a在拟南芥也引起雄性不育，它们诱导花药表型的形成不同于那些myb33 myb65。这些发现表明，其他成员拟南芥GAMYB样家族被miR159靶向并参与花药发育[70]。许多的GAMYB在本研究中发现了类似转录本(附加文件)9)。与之前的发现一致，我们的结果表明GAMYB在雄性花蕾中，样转录本上调，miR159下调。这些发现提高了mirna调控的可能性GAMYB样基因参与雄性芦笋花药发育，miRNA上调和下调GAMYB类基因可能参与了芦笋雌花芽雄蕊败育过程。一些未知的deg也被民主党盯上了。因此，这些基因可能与性别决定和分化有关，但这种可能性尚未得到证实。

应用RNA-Seq技术系统地研究了雌雄异株芦笋生殖器官中性别偏倚基因的表达。结合之前的研究数据，我们的结果显示雄性/超级雄性和雌性花蕾之间的差异为4876度。许多与植物生殖、植物激素信号和TFs相关的候选基因被表征，并涉及复杂的性别决定和分化网络。我们的转录组分析为芦笋性别表达的遗传调控提供了新的见解。本研究也为进一步开展这一重要的雌雄异株蔬菜作物的基因组研究奠定了基础。

植物材料

本研究使用的雄性和雌性芦笋(品种UC309)在河南师范大学的温室中生长。从一个完整的兄弟后代中抽取了6个雄性和6个雌性芦笋个体，代表两种性别类型的两个重复。花蕾(直径0.5-0.7毫米)，可能代表每个性别类型的关键性别决定阶段[21，71]，立即用液氮冷冻，- 80°C保存备用。

Illumina测序和数据处理

用TRIZOL试剂(Life Technologies, CA, USA)按照生产商的说明书从芦笋花蕾中提取总RNA。在每次重复中，从三株植物的花蕾混合物中分离出等量的RNA。使用NEBNext®Ultra™定向RNA文库准备试剂盒for Illumina®(NEB, USA)按照制造商的推荐使用，每个样品约3 μg RNA构建测序文库。转录组测序在Illumina Hiseq 2500平台上进行，产生125 bp的成对末端reads。

快速格式的原始读取，包括我们的测序数据和先前报道的花蕾数据[6]，通过内部perl脚本进行处理。在这一步中，通过从原始数据中去除包含适配器的读取、长度至少为10ns的读取和低质量读取，获得干净数据(clean reads)。下游分析基于高质量的清洁数据。使用Trinity进行De novo组装，min_kmer_cov默认设置为2，其他参数均默认设置[72]。为了评估转录组的完整性，BUSCO v3 [73]通过使用植物同源数据库，其中包括1440个基因。我们还使用了Transrate v1.0.3 [74]来评估转录组组装质量。

基因注释与功能分类

利用BLAST对以下5个数据库进行基因功能标注，截止e值为10⁻⁵: Nr (NCBI非冗余蛋白序列)，KOG/COG(同源蛋白群)，Swiss-Prot(手动注释和审查的蛋白质序列数据库)，KEGG Ortholog数据库和GO。

花相关基因分析拟南芥数据

一个包含大量花相关基因的数据集拟南芥为研究芦笋花蕾组装单基因中与花相关的基因而产生。这些数据主要来自植物繁殖数据库[75]和花卉相关数据库[27]。剔除冗余基因后，无冗余基因11398个拟南芥该数据集涉及基因位点。的拟南芥从拟南芥信息资源(The Arabidopsis Information Resource, TAIR)下载蛋白质组，并设置为参考，以芦笋组装单基因作为查询进行局部BLASTX分析。

差异表达分析

利用SOAPaligner/soap2将测序reads映射到组装的序列上[76]。通过计算不同样本中每个单基因的RPKM值，将reads计数归一化。使用DESeq R软件包(1.10.1)进行差异表达分析[77]。p-value调整为多次测试的错误发现率(FDR) [78]。在表达上有最小两倍差异的基因₂Ratio|≥1)，FDR≤0.05为deg。我们的数据和Harkess的减数分裂前花蕾数据的交集被确定为雌性和雄性/超雄性偏异基因，用于进一步分析。然后在基于Wallenius非中心超几何分布的GO序列R包中对这些deg进行GO富集分析[79]，可以根据deg的基因长度偏差进行调整。使用KOBAS软件评估KEGG通路中deg的统计富集程度[80]。利用HMMER v3.1b2 (http://hmmer.org/)，以鉴定雌雄花蕾之间差异表达的tf。检测并进一步检查潜在的性别特异性或差异表达的tf。

中存在的验证

我们选择了24个deg进行定量RT-PCR分析，使用来自不同个体的花蕾RNA样本进行转录组测序，以验证我们的Illumina测序数据。利用引物6(附加文件)根据参考单基因序列设计基因特异性引物11)。以18S rDNA基因为内对照进行qPCR分析。cDNA合成和q-PCR使用先前描述的程序进行[81]。用2^−ΔΔCT方法(82]。

小RNA测序和miRNA鉴定

构建了芦笋雌雄花蕾独立的小RNA文库，并对其进行了测序。每个文库中有3株雄性单株或3株相同数量的雌性单株。采用NEBNext Multiplex Small RNA library Prep Set for Illumina (NEB, USA)软件，用总RNA (3 μg)构建每个男性和女性样本的sRNA文库。制备的文库在Illumina Hiseq 2500平台上测序，生成50 bp的单端reads。测序数据包括转录组测序和sRNA测序结果，录入NCBI SRA数据库，登录号为SRR5112683、SRR5112684、SRR5118382和SRR5118383。

已知和新型mirna的鉴定和差异表达分析

使用Illumina的Genome Analyzer Pipeline V1.5对原始读数进行处理和过滤。筛选后的小RNA序列与芦笋mrna、Rfam (http://rfam.xfam.org/)，以及Repbase (http://www.girinst.org/repbase/)，以丢弃mRNA、rRNA、tRNA、snRNA、snoRNA和重复序列。已知mirna和新型mirna是使用先前研究描述的方法鉴定的[83]。

利用DEGseq对雌雄花蕾中的mirna进行差异表达分析[84R包。p-value使用问值(85]。问-value <0.01和|log₂(foldchange)| > 1默认为显著差异表达分析的阈值。

miRNA靶标预测

植物miRNA靶标预测服务器(http://plantgrn.noble.org/psRNATarget/)根据芦笋转录组的组装数据，使用默认参数预测推测的miRNA靶基因。

爆炸:

基本的局部对齐搜索工具

齿轮:

同源的蛋白质群

度:

差异表达基因

民主党:

差异表达的miRNA

小块土地:

乙烯反应转录因子

罗斯福:

错误发现率

走:

基因本体论

嗯:

隐马尔可夫模型

是:

茉莉酸

KEGG:

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

microrna的:

微

Nr:

NCBI非冗余蛋白序列

聚合酶链反应:

聚合酶链反应

RPKM:

每千碱基每百万的读取量

之后:

性别特异性表达基因

TAIR:

拟南芥信息资源

TF:

转录因子

作者感谢Novogene Bioinformatics Technology (Beijing, China)和Dr. Minshan Sun (Genedenovo, Guangzhou, China)对高通量测序和数据分析的协助。

资金

本研究得到国家自然科学基金项目(31470334和31300202)、河南大学科技创新团队项目(17IRTSTHN017)和国家重点研发计划项目(2016YFD0300203-3)的资助。

数据和材料的可用性

原始测序数据包括转录组测序和sRNA测序，录入NCBI SRA数据库，登录号为SRR5112683、SRR5112684、SRR5118382和SRR5118383。种子将根据要求提供。

从属关系

1. 河南师范大学生命科学学院，新乡453007

   李淑芬，张学进，袁金红，邓传亮，高武军

2. 新乡医学院基础医学院，新乡453003

   郭军张

贡献

SL分析数据并撰写稿件，GZ参与数据分析，XZ和SL进行实验，JY制备mRNA进行测序，CD对数据分析提供有益建议，WG设计项目并修改稿件。所有作者都阅读并批准了最终的手稿。

相应的作者

对应到胡军高。

伦理批准并同意参与

不适用。

发表同意书

不适用。

相互竞争的利益

作者宣称他们没有竞争利益。

出版商的注意

施普林格·自然对已出版的地图和机构关系中的管辖权要求保持中立。

开放获取本文遵循知识共享署名4.0国际许可协议(http://creativecommons.org/licenses/by/4.0/)，它允许在任何媒体上不受限制地使用、分发和复制，前提是您要适当地注明原作者和来源，提供到知识共享许可协议的链接，并注明是否进行了修改。创作共用公共领域免责声明(http://creativecommons.org/publicdomain/zero/1.0/)适用于本文中提供的数据，除非另有说明。

转载及权限