花青素合成基因的差异表达是异域花颜色差异的基础杜鹃sanguineum复杂的

背景

的杜鹃sanguineum它是中国云南西北部和西藏东南部高山特有的。在这个复杂的品种显示独特的花朵颜色，甚至在花蕾阶段。然而，花色变化的分子调控机制还没有很好地阐明。在这里，我们通过测量花的反射率剖面和比较转录组分析，对三个共存的品种进行了研究r . sanguineum复杂的，花分别为黄红粉红色、亮深红色和深黑深红色。我们利用RNA-seq和qRT-PCR数据比较了花青素/类黄酮生物合成途径差异表达基因(DEGs)的表达水平。我们基于转录组衍生的单核苷酸多态性(SNPs)数据进行聚类分析，最后分析了DEGs的启动子结构。

结果

在400 ~ 700 nm范围内，三种颜色形态的反射光谱变化明显，在饱和度、亮度、色相和饱和度/色相比(花青素含量的间接测量指标)上存在明显差异。我们鉴定了三个品种共有的15164个正位类群。SNP聚类分析表明，这些品种并非单系。共有40个同源基因编码12种酶参与了花色多态性。这些花青素生物合成相关基因与花青素的合成、修饰和运输特性相关(RsCHS，RsCHI，RsF3H，RsF3′H，RsFLS，RsANS，RsAT，RsOMT，RsGST)，以及参与分解代谢和降解的基因(RsBGLU，负责，RsCAD)．序列和序列的变化独联体这些基因的作用元件可能与花青素的积累有关，从而导致了花青素的分化r . sanguineum复杂。

结论

结果表明，该品种亲缘关系非常密切，花色变异也非常明显r . sanguineum复合体与参与花青素合成和分解合成网络的基因的不同表达水平密切相关。我们的研究提供了一个复杂的场景，包括花的颜色遗传机制之间的关系，伴随基因流动的品种，这可能代表了传粉者介导的早期同域物种形成的一个早期案例。

数百年来，植物学家、生态学家和园艺学家都为花卉颜色的显著多样性，尤其是野生植物的多样性而着迷。1，2，3.］．花器官的着色是开花植物的一个显著特征，也是被子植物辐射的一个显著特征[4，5］．花朵颜色多样性被认为是主要与传粉者(如昆虫、鸟类)和动物有关的种子传播的关键适应性特征之一[6，7］．此外，花朵颜色表型是植物分类学家用来分类的一个重要特征。然而，花色在进化上似乎是最不稳定的性状之一，甚至在同一物种的种群中也是如此。7，8］．

花朵的细胞化合物有助于颜色的轮廓和视觉上被人类感知，通常被称为“色素”。一类属于类黄酮的次生代谢物是植物着色色素的主要决定因素，其中花青素负责花、叶、果实、种子和其他组织中从红橙到红、从紫到紫的色素[9，10］．花青素是花卉着色的主要化合物，存在于90%以上的被子植物中[11］．导致花青素积累的类黄酮生物合成途径具有高度的保守性和良好的特征，在许多物种中得到了广泛的研究，其中大多数是模式植物或具有重要农业和园艺价值的植物[12，13，14，15］．很少有研究探讨野生物种花青素形成和积累的分子基础[16，17］．基于这些研究，提出了涉及花青素积累的三个主要相关因素，包括发生在细胞核中的转录调节基因(MYB-bHLH-WD40复合体)、结构基因(CHS，FLS的，DFR,答)作用于生物合成途径，以及转运基因(销售税)将花青素从胞浆转移到液泡[10，18，19］．这些基因的表达也可能受到序列和自然变异的影响独联体-启动子区域的调控元件以及表观遗传修饰(如DNA甲基化)[18，20.］．此外，黄酮醇对花青素的共色素作用可以稳定和增强花的颜色，这种作用被观察到为深色素效应，即花青素含量的强度得到加强[21］．例如,DFR基因和FLS的基因可以通过两个主要分支竞争导致产生不同花青素和黄烷醇的底物[22，23]，从而导致共色素沉着。与生物合成途径相比，植物中花青素分解代谢的知识有限。一些分解代谢基因BGLU而且每已经被证明是导致花青素降解的原因[24，25］．然而，调节花青素合成的分子机制已被证明在不同的植物种类中存在差异，从而导致花青素的结构多样性，因为生物合成途径是通过调节网络受到多个因素的调节[26］．

颜色是可见光谱范围内的一种电磁辐射形式。色素反射的波长决定了花的颜色[27］．颜色可以根据亮度(信号的强度，B)、饱和度(颜色的纯度，S)和色相(颜色的光谱描述符，H)来定义和分类，这些特征通常用于区分颜色[27，28］．亮度是指由花青素含量决定的颜色强度[29，30.]，不同颜色成分组合如B/H、S/H也与花青素含量显著相关[31］．刘等人。[32]提出颜色亮度随花青素含量的增加而降低。研究还表明，饱和度/色相比(S/H)与花青素含量之间存在相关性[31］．利用这些参数，可以快速、无损地测定花青素含量。

常绿杜鹃花(杜鹃)，花青素和黄烷醇是花色素的主要成分，尤其是花青素成分(即花青素、飞燕草苷、麦草苷、天竺葵苷、芍药苷、芍药苷)的组成，其数量决定了花的颜色从浅粉色到紫色不等[11，33］．一些研究已经报道了这一点r . kiusianum紫色的花含有花青素和飞燕草素的衍生物，而红色的花r . kaempferi只含蓝靛素衍生物[34］．勒迈特等人[35，36)研究艾丽卡种，属于同一科的鸢尾科杜鹃，利用qRT-PCR和UPLC-MS，揭示了红色或粉色花与白色或黄色花的花色转移的花青素遗传网络，发现单通路基因表达丢失，或通路基因中某一转录因子表达丢失或功能突变丢失导致整个通路的失效，从而导致显著的花色转移。

在此，我们利用高度的颜色多态性研究了花着色的遗传基础杜鹃sanguineum复杂。复杂的(R．subgen。Hymenanthes)包括花由黄至粉或由深红色至暗红色的植物，它们主要根据花的颜色差异被分为六个品种[37］．该建筑群的成员基本上位于与积雪有关的高海拔地区(> 3000米)[37］．它们是全球生物多样性热点之一的云南西北部和西藏东南部的特有物种[38］．这一地区也被认为是多样化和差异化的中心之一杜鹃(37，39］．本属的花朵颜色多态性传统上被视为一种生态适应性特征，对吸引特定的传粉者至关重要[40，41，42]，亦可能是对环境变化的反应，例如不同海拔、温度和土壤条件下的紫外线辐射[32］．虽然有研究发表了花青素的成分和含量杜鹃花，大部分专门用于用薄层色谱法(TLC)和高效液相色谱法(HPLC)鉴定一些野生和栽培杜鹃花花瓣中的色素成分[11，33］．到目前为止，还没有研究集中在植物的种下颜色多态性的分子机制杜鹃．研究密切相关的实体，如物种复合体，具有相当同质的遗传背景的优势，花的颜色基因不同，同质性的情况是有限的。以前的研究主要集中在单一物种不同发育阶段的颜色变化[14，18]，或涵盖若干相关物种[26，35］．

在本研究中，我们将转录组测序(RNA-seq)和基因组重测序与反射光谱分析相结合，以阐明三种不同颜色的天然品种之间的分子和花青素含量的差异r . sanguineum繁花，黄红粉红至深黑深红色。本研究旨在研究花青素/类黄酮生物合成途径候选基因表达与花种下花色变异的相关性。我们的研究结果可能为花朵颜色变异表达的遗传机制提出假说，并为植物中传粉者介导的早期同域物种形成提供了一个代表性案例r . sanguineum复杂。此外，这是第一个比较非模式物种在自然系统中的转录组谱的研究杜鹃．

反射光谱和颜色形态的差异

采样的个体可以分为三种颜色类别代表r . sanguineumvar。sanguineum（RsS鲜红的花朵。1d),R。var。haemaleum（RsH那里的花是深黑深红色的。1e),R。var。didymoides（标准偏差粉红色的花管和红色的花瓣。1f).三种颜色变体的反射率光谱在波长范围内变化明显，特别是在400 ~ 620 nm之间(图5)。2一个)。RsS红光谱(约620 nm)有一个离散峰RsH在更远的红色范围(> 700nm)有一个峰。的反射标准偏差在蓝到黄光谱范围内增加，并在橙红色范围内达到峰值(~ 600 nm)(图5)。2a)三个品种的花色在亮度、饱和度和色相上有显著差异。标准偏差亮度值最高的是RsS,RsH与最小值，而RsH饱和度最高，是吗RsS与标准偏差(图的最低水平。2b, c).色相值在样品中最高标准偏差,然后RsS而且RsH(无花果。2b).花青素的饱和度和色相的比值(花青素含量的间接测量)RsH是最高的与那些标准偏差而且RsS样品(图。2d).的比值值标准偏差样品的差异比其他两个品种要大得多，这可能与它的双色花的颜色组成的变化有关(图。1)．

转录组从头组装和质量评估

我们从两种组织(花芽和叶芽)中为三个个体分别测序了共18个RNA-seq文库r . sanguineum使用Illumina配对端测序的复合物。经过质量控制后，大约有6.558亿M (M)的清洁读取(约96千兆对，Gbp)，库之间的读取数量非常一致，在32,111,674和43,353,842之间S1)．3个品种的CORSET的contig N50长度相似，碱基对(bp)在903 ~ 1125之间，转录本数在117,976 ~ 171,725之间(表2)S1)．将18个库的reads与品种特异性参考文献进行比对，比对率为84.25 ~ 94.85%，平均比对率分别为92.82、90.30和93.77%R。var。sanguineum，R。var。haemaleum而且R。var。didymoides分别(表S2)．基于1614个保守的BUSCO胚胎植物同源序列，对其转录组完整性进行评估，确定了1540个(95.4%)完整和片段化的BUSCOR。var。sanguineum， 1548 (95.9%)R。var。haemaleum1538人(95.3%)R。var。didymoides(图S1一、表S3)．因此，我们的结果表明，三个转录组组装良好，相对完整。

正交鉴别与功能注释

ORF预测发现，在最终的非冗余转录本中，共有53,207个蛋白质编码转录本R。var。sanguineum46754年,R。var。haemaleum和38548年R。var。didymoides．对于三个品种，OrthoVenn2 [43]在总共138,509个蛋白质编码转录本中鉴定出31,525个簇，其中16,361个是同源基因(包含至少两个物种)，15,164个是单拷贝基因。有54,192个单例未包含在任何集群中(表S4)．基于15,164个一对一的单拷贝正交关系的注释(表S53个品种中，我们在国家生物技术信息中心非冗余(NR)蛋白数据库中找到14441个(95.23%)匹配项，UniProt/Swiss-Prot中找到10648个(70.22%)匹配项，COG/KOG中找到13089个(86.32%)匹配项，eggNOG中找到14109个(93.04%)匹配项，基因本体(GO)中找到7445个(49.1%)匹配项，京都基因与基因组百科(KEGG)中找到4302个(28.37%)匹配项(图)S1b,图S2)．

SNP检测和聚类分析

以15164个同源单拷贝基因为参考序列集，构建了该9个个体的cDNA文库r . sanguineum复合物共产生50853个snp。基于这个SNP数据集，snphyo树和PCA分析都没有将样本聚类到特定品种的聚类中(图S3)．只有RsH样本落在SNPhylo树的一个分支中，但也有一个样本落在标准偏差(废话= 100%)。标准偏差而且RsS样品混合在两个枝中。

整体差异表达基因和功能富集

我们进行了三个成对的转录组比较，以确定在花的颜色变体中的DEGsr . sanguineum通过计算基因的FPKM值。在三向比较中共检测到2148个差异表达基因(DEGs)。特别是对比深黑深红色的变形RsH带着黄色的粉红色的变形标准偏差鉴定出701个DEGs(413个上调，288个下调)RsH）;的比较RsH对比深黑深红色的变形RsS共发现1378个DEGs，其中851个表达上调，527个表达下调RsH）;的比较RsS与标准偏差给出1034个DEGs(513个上调，521个下调)RsS)(图。3.a, b).差异表达基因最高RsS而且RsH．所有deg模式均在补充图中给出S4＆S5．基于PCA和热图结果的基因表达谱(图S4；数字S6a)，每个品种聚集在一起(显示种特异性模式)，这表明RNA-seq文库是可靠的。

所有的DEGs被划分为三个主要的基因本体类别:生物过程(BP)、细胞成分(CC)和分子功能(MF)，其中GO术语大部分属于生物过程，其次是分子功能和细胞成分。三个两两比较结果显示，共有40个显著GO术语，包括27个BPs、8个CCs和5个mf(图5)S6b, c).来自不同比较的基因清楚地表明了相同的分子和细胞事件，如植物抗性(GO:0006952，防御反应)和代谢过程(GO:0009737，脱落酸反应;GO:0009753，茉莉酸反应)。KEGG通路富集分析结果显示，DEGs与多种代谢和生物合成通路相关:分别有20、17和29个通路富集，分别对应701、1378和1034个DEGsRsHvs。标准偏差，RsHvs。RsS,RsSvs。标准偏差,分别。特别是编码苯丙烷类生物合成酶的基因(ko00940)均显著富集。此外，与类黄酮生物合成(ko00941)和细胞色素P450 (ko00980)代谢异种生物相关的途径也被丰富RsHvs。标准偏差而且RsHvs。RsS．图中显示了一些代表性的最显著富集的KEGG通路。3.c.为了更深入地了解所有deg的表达模式，使用TMM归一化表达值生成了一个热图(图1)S6A)，这里我们关注的是与色素相关的基因。

与颜色多态性相关的候选基因r . sanguineum复杂的

基于DEGs的KEGG途径富集的结果，并考虑去除极低表达的基因(FPKM < 1)，我们从40个同源基因中确定了13个候选基因，这些基因可能与花的色素沉积相关，有助于该复合物中的花青素积累和/或共色素沉积(表1)S6)．这些都是花青素相关基因，参与花青素生物合成、花青素修饰和花青素转运。其中包括花青素合酶(答)、酰基转移酶(在), beta-glucosidase (BGLU)、4-香豆酸辅酶a连接酶(4 cl)、肉桂醇脱氢酶(计算机辅助设计)、查尔酮异构酶(气)、查尔酮合酶(CHS)、黄酮3-羟化酶(F3H)、类黄酮3 ' -羟化酶(F3′H)、黄酮醇合酶(FLS的)、谷胱甘肽s -转移酶(销售税), O-methyltransferase (OMT)和过氧化物酶(每)．在这其中,Rs4CL，RsCHS，RsF3H，RsF3'H而且RsOMT只有一本，其他的(RsANS，RsAT，RsBGLU，RsCAD，RsCHI，RsFLS，RsGST,负责)代表多基因家族有2个(例如:RsAT，RsCAD，RsCHI)至8或9份(例如:RsBGLU，RsGST，负责)．对这些基因序列进行两两比较，鉴定出23个DEGs, 25个DEGs和12个DEGsRsHvs。标准偏差，RsHvs。RsS而且RsSvs。标准偏差分别(无花果。3.d;数据S7，S8，S9)．获得了13个候选基因的代表性拷贝的qRT-PCR验证候选基因的可靠性表达谱，除Rs4CL未能放大(表S7)．

如图所示。4，结果表明，绝大多数花青素合成、修饰和转移酶的基因在中表达较高RsH,在标准偏差而且RsS，分配给水解酶的基因下调RsH．两个类群共有6种酶共有9个基因(RsHvs。标准偏差而且RsHvs。RsS),包括RsBGLU1/3，RsPER1/3/4，RsCHI1，RsFLS1，RsANS2而且RsAT2(无花果。3.d;数字S7，S8)．然而，只有少数的deg，包括Rs4CL，RsBGLU7而且RsPER1表现出更高的表达RsS在标准偏差．此外,在标准偏差，RsBGLU4/5，RsPER5/7/8，RsCAD1/2，RsAT1，RsOMT的表达高于RsS(无花果。3.d;无花果。4；数字S9)．假设花青素生物合成途径与花的颜色有关的初步示意图r . sanguineum如图所示。4．

花青素通路相关基因启动子区序列多样性

每个基因的启动子序列(翻译起始位点上游2 kb)在序列相似性(98.6 - 99.7%)和GC含量上都表现出高度的保守性，尽管在某些位置上存在差异(表2)S8；数字S10)．的详细分析顺式以MYB转录结合位点为靶点的这些启动子的作用元件显示，除了RsANS而且RsF3'H在美国，其他11个候选基因都至少在一个基因上有变异顺式图案(图S10)．这些启动子在其他类型中也有变异顺式主题(表S8)．因此，这些数据提供了启动子区域(序列的变异和独联体-作用元件)和各自花青素基因的差异表达。

类黄酮/花青素代谢被认为是促进花色素形成的最重要途径之一，由一种多酶复合物催化[10］．花卉颜色的遗传基础最近在物种水平上得到了研究艾丽卡，与…密切相关的属杜鹃在杜鹃花科(35，36］．这些研究表明，在颜色变化的基础上，单通路基因的表达丢失，或转录因子丢失或功能突变丢失导致的通路崩溃[35］．虽然这些研究中的颜色变化涉及物种之间的变异，但在这项研究中，我们努力阐明更密切相关的进化实体之间的颜色变化的遗传学，在不同的物种之间r . sanguineum复杂。这些变化可能更加微妙，因为进化距离比物种之间的距离更小，可能更不受同形性的影响，反映了更准确的进化路径变化。之前的研究的一个缺点是缺乏转录组数据，无法检测到基因同源谬误[35，36］．事实上，在目前的研究中，我们发现许多花青素基因是小基因家族的一部分(图。3.D)，因此更精确地解开花青素途径。

通过比较转录组和确证qRT-PCR实验，我们发现在具有相似反射谱差异特征的类群中，差异表达基因的模式相似。RsHvs。标准偏差而且RsHvs。RsS．因此，当比较两个反射率光谱更相似的类群时，发现的deg数量更少，RsSvs。标准偏差．因此，密切相关的类群之间表型差异(即花青素含量差异)的大小可能是各自花组织中颜色基因表达水平变化的结果。在这里，我们确定了两个类群中共有的6种酶的9个基因，RsHvs。标准偏差而且RsHvs。RsS(图S7，S8)．在我们的反射光谱结果中，亮度参数(B)(通常指强度)是最低的，饱和度(S)和色相(H)的比值是最高的RsH(无花果。2C, d)，表明花青素含量高[29，31］．同样，之前的研究表明，花的颜色强度在很大程度上是由花青素的数量决定的。44，45]，即花青素含量越高，花色越深。

高水平的表达会使颜色更加鲜艳。46我们发现了不同的基因，可能促成了深黑深红色的花R。var．haemalum（RsH)．我们的研究结果表明RsCHI1，RsFLS1，RsANS2而且RsAT2在类群比较中高度表达RsH在该品种中显著表达，对应其深黑深红色的花(图。3.d;无花果。4)．这一发现与蓝莓果实发育的情况类似，在蓝莓果实发育过程中发现了通路基因表达与花青素产生之间的相关性[46］．在类黄酮/花青素途径中，气高度表达的RsH，被证明可以将查尔酮转化为黄酮。的答基因，也只在RsH，在催化无色白花青素合成为彩色花青素的过程中起着关键作用。例如，张等人。[47)报道,气在生菜的红叶中表达量显著高于在浅色叶中表达量(摘要以)．赵等人。[48的低表达答的芍药导致花青素积累减少，导致颜色变黄。在我们的研究中,RsANS可能是花青素积累较弱的原因标准偏差，而高积累在RsH在后者达到相对较高的花瓣色彩强度。一般来说，花青素的结构由色素积累而成答催化是不稳定的，通过花青素酰基转移酶(在)形成稳定的花青素[10］．据报道，当在被抑制，花瓣变得苍白[49］．这在我们的结果中得到了证明RsAT基因在RsH亮度低(花色强度高)，其他两个花冠较浅的品种亮度较低。有趣的是,FLS的基因在RsH．FLS的竞争相同的衬底DFR，从而产生不同的类黄酮和花青素[10，49]，表明共色素作用[50］．由于共色素沉淀而引起的颜色变化增加了花的颜色强度[51］．在我们的研究结果中，辅色素沉着，通常与黄酮醇有关RsFLS，可能导致强烈的颜色花青素，使颜色转变为深黑深红色。这与研究的红紫色花瓣是一致的r . simsii花(52］．

在比较RsSvs。标准偏差对比转录组分析表明，与花色生物合成相关的6个基因中只有12个片段存在差异表达(图1)S9)．这很可能是相对较小的光谱反射差(即花色差)之间的差异RsS而且标准偏差，这就是与花青素生物合成相关的DEGs数量较低的原因。高表达的六个基因中的三个标准偏差,即RsBGLU，RsCAD,负责它们属于分解代谢酶，因此可能是颜色较浅的原因。我们的发现与Yang等人之前的研究一致。[53他对两个牡丹品种的花色变化进行了研究，结果表明，负责花色从红色变为橙色和黄色的花青素含量急剧下降。同样的，朱等人。[54的差异表达NnOMTs可能与花瓣的颜色差异有关吧莲属椰子具有黄色和白色花的栽培品种。但这些可能与Le Maitre等人的发现不同。[35，36]艾丽卡开黄色花的种。对于这些，他们发现要么基因表达正常，要么基因表达缺失F3′HMYB识别元件的缺失。其他色素，如类胡萝卜素也会影响花的颜色[10］．然而，这些在植物中很少被检测到杜鹃(11，33］．其他基因，比如RsCHS，RsF3H，RsF3′H而且RsGST也被发现在RsH与两者相比标准偏差或RsS．这些是单拷贝基因，很可能解释了这个品种的深色花朵。

我们的结果进一步表明RsBGLU1/3，RsPER1/3/4与代谢酶相关的基因显著升高标准偏差和更高的RsS相比RsH．在类黄酮/花青素途径(图。4)，这两个基因被划分为与苯丙类早期生物合成途径相关的“早期生物合成基因”(egs)。它们分别参与了肉桂酸和p -香豆素酰辅酶a水解为香豆素和木质素[13］．β-glucosidases (BGLU)及过氧化物酶(每)已经被证明是许多植物花青素降解的原因[24，25］．一般来说，花青素的积累是由生物合成和降解的平衡所决定的。在我们的研究中，品种标准偏差而且RsS可能由于早期底物的分解代谢导致花青素积累减少[24]和黄酮/花青素通路中成熟花青素的降解，从而导致花的颜色强度低。很可能这些品种的花的颜色与通路中的通量转移(花青素含量)有关，从而导致了颜色强度的差异。

然而，上述所有候选基因都被归类为结构基因。MBW (MYB-bHLH-WD40)蛋白复合体基因表达差异不显著。一种可能的解释可能指向这里的采样方法，只使用一个阶段(花蕾后期)。杨等人。[19]表明，MBW转录因子在转录水平上共同调节花青素的积累，特别是在花着色的最初阶段，而可能不存在于花蕾后期。此外，MBW复合体可以与花青素基因的启动子结合，启动子序列的MYB结合位点在花青素合成中起着至关重要的作用[55］．在本研究中，三个品种各基因的启动子序列总体上是保守的，但在基因序列上有明显的变异独联体-基序如MYB结合位点在11个基因中，除RsANS而且RsF3′H(表S8；数字S10)．我们的结果也可能在某种程度上支持今晚的假设[56在那里，表型差异，特别是对于密切相关的物种，更有可能由独联体-调控区域比基因的蛋白质编码区域更明显，但它们被认为相对不受适应度的负多效性影响，而转录因子基因的突变变化是最不可能的，因为其潜在的影响更广泛[35］．此外，一些研究表明甲基化水平在启动子区域MdGST(18),OgCHS(57与基因表达水平和花青素多样性显著相关。然而，各基因启动子区域的GC含量在品种间具有相似的水平(图S10)，这表明DNA甲基化不太可能是颜色多样性的原因r . sanguineum复杂。因此，我们认为花青素通路基因的调控和转录可能不是独立的过程，虽然它们分别表达，但它们共同促进了花青素的积累[55，58］．尽管如此，我们在黄酮类/花青素通路方面有了初步的发现r . sanguineum复杂的，但需要更多的实验来进一步研究每个遗传成分在这种情况下的贡献。

有人提出，花朵颜色的进化转变通常归因于传粉者介导的选择，这可能有助于生态位的占领和随后的生殖隔离和多样化/物种形成[5］．例如，蜜蜂有三种光感受器，在光谱的UV、蓝色和绿色范围内达到峰值，分别对应于344、438和560 nm [59]，而鸟类是四色的，并有对600至620纳米红光敏感的受体[60］．

在一些授粉研究中杜鹃研究发现，鸟类、蜜蜂、蝴蝶和鞘形蛾是主要的传粉者杜鹃物种(40，41，42］．宋等人[61]提出太阳鸟和大黄蜂是潜在的传粉者R．成为而且R．edgeworthii分别有红花和白花，而Epps等[62]观察到蝴蝶在开橙黄色花的r . calendulaceum．这说明传粉者表现出对花朵颜色的偏好，颜色相近的花朵往往会吸引特定的传粉者，即使是在相同的种群之间或密切相关的物种之间[6，59］．因此，传粉者组合的时空波动可能导致花朵颜色的变化，反之亦然[5］．然而，在一些密切相关的同域类群中，不完全生殖障碍有助于杂交，并在基因流动中发挥关键作用[63］．低繁殖隔离的杂交可以频繁地传递颜色基因，特别是在有限的基因库中，如这里所示的品种r . sanguineum因为它们不形成单系基团(图S3)．我们的研究结果与一些报告一致，这些报告显示，相对简单的基因变化导致了花朵颜色的重要变化。64］．我们已经在这里的密切相关的变种中推断了这些。这种情况可能会引起传粉者聚集的转变，繁殖隔离随着时间的推移而加强。因此，我们的研究可能证明了一种由花朵颜色分化诱导的早期同域物种形成模式，正如通过行为隔离观察到的鸟类[65]，以及植物物种的地质-土壤分化[66]，或地理差异[67］．

本研究采用Illumina转录组测序(RNA-seq)和基因组重测序技术，对田间采集的3个品种的花色基因表达变化进行分析r . sanguineum复杂。本研究初步揭示了黄花菜花色分化的遗传机制r . sanguineum复杂。我们的研究结果表明，花色的变异与表达水平的差异有关，并在一定程度上影响花色的表达顺式作用调节花青素生物合成基因和花青素降解基因，而不是功能突变丢失。深黑的深红色花开了R．var。haemalum表现出几乎所有花青素基因的高表达，而明亮的深红色开花R．var。sanguineum在花青素合成的最初步骤中，是否有高水平的基因表达和黄色的粉红色花朵R．var。didymoides大多数花青素基因有低或中表达水平。后者花青素降解基因表达量较高，可能增加了其花中花青素含量较低的原因。这一发现与之前的一项研究不同艾丽卡花青素基因和MYB识别元件的移码突变被发现是红黄颜色转移的原因。可以提出一些假设来解释复合体中的花色变化:这些品种关系太近，花青素基因或转录因子的编码区域没有获得性突变，它们的花色变化是由于可诱导的表达水平的差异独联体-作用调控，以及合成代谢和分解代谢花青素途径中基因之间的平衡。该复合体目前在品种间存在复杂的杂交和基因流动关系，但如果传粉者偏好强化进一步发展，可能发展为传粉者驱动的早期同域物种形成。

样品收集

三个品种的样品r . sanguineum复杂的(图。1),即r . sanguineumvar。sanguineum（RsS，开着鲜红的花)，R。var。haemaleum（RsH，开着深黑深红色的花)，R。var。didymoides（标准偏差于2018年6月在云南西北部的高黎贡山(N 27°47′11.40″，E 98°27′35.28″)采集。花蕾后期的花组织。1b)和叶芽(图。1C)从每个品种的三个个体中取样，立即用液氮冷冻。来自相同个体的叶片组织被取样，同时用硅胶干燥，花也是如此，用于花冠反射光谱测量。这9个个体的采样范围很小(20米)²)在类似的气候条件和环境因素(例如土壤、温度、降水和阳光辐射)下。每个人的代金券被收集并存放在中国科学院昆明植物研究所植物标本室。

花反射光谱测量

植物色素对可见光谱中的光的吸收产生了独特的光谱反射特征。为了获得三个品种不同花色的反射率光谱，采用分光光度法定量花青素花冠色素[11，68］．在田间测量了5个新鲜的、健康的、完全开放的花(每个品种15个，共45个样品)。两个独立的光谱读数从花的相同位置，花冠管的一半。在200-800 nm范围内使用USB2000+光谱仪测量漫反射光谱，采用氘/钨卤素光源(Ocean Optics，达尼丁，美国佛罗里达州)，集成时间为3 s，箱座为12。在每次样品测量前扫描真黑和真白对照参考文献。颜色数据处理采用光学2009 SpectraSuite(海洋光学)软件。光谱被截断到300-700 nm，并使用孔雀座包中R(69，70］．用饱和度与色相的比值(S/H)来测定花的花青素含量[31］．

RNA和DNA提取和Illumina测序

根据制造商的协议，使用Spectrum TM植物总RNA Kit (STRN250, Sigma)从花芽和叶芽组织中分别提取和纯化总RNA。对每个品种的花芽和叶芽材料进行三次生物重复(来自三株植物)。用DNeasy Plant kit (QIAGEN)从叶片中提取基因组DNA。RNA和DNA的质量和数量用NanoDrop 2000分光光度计(Thermo Fisher Scientific, Waltham, MA, USA)进行评估。根据制造商的协议，使用NEBNext Ultra™RNA/DNA文库准备试剂盒(NEB, MA, USA)生成转录组和基因组文库，并在Illumina HiSeq X Ten测序平台(San Diego, CA, USA)上进行测序，分别为每个RNA和DNA文库生成约6gb和7gb (~ 10x)的对端reads (2 × 150bp)。文库制备和Illumina测序在诺健生物信息技术有限公司(北京)完成。

数据处理，从头组装和映射

首先对原始数据进行过滤，去除带有适配器序列的reads、含有poly-Ns的reads、含有歧义核苷酸的reads和低质量的reads，然后使用SOAPnuke计算phred评分(Q20, Q30) [71］．使用FastQC评估剩余读取的质量，包括gc含量、序列长度分布和清洁数据的序列重复水平[72］．所有后续的分析都是基于干净的数据。

为了获得参考水平的转录组组装，从每个品种(叶芽、花芽，在三个重复中)的清洁reads进行组合，并使用Trinity v. 2.6.5重新组装[73］．除了最小kmer覆盖率(设置为2)和最小contig长度(设置为200)外，所有Trinity参数都设置为默认值。程序集统计信息是使用Trinity包中的TrinityStats.pl脚本获得的。我们还使用了HISAT2 v. 2.1.0 [74]来评估组装质量，通过将读取映射回组装的转录本来计算总体对齐率。

程序集筛选和完整性评估

为减少冗余和复杂性，我们对所有组装的转录本进行过滤，方法如下:首先，我们使用CD-HIT-EST v. 4.7.0 [75]，将字长设置为10，序列标识阈值设置为0.95，以去除重复。然后，紧身胸衣v. 1.07 [76]用于聚类转录本序列，过滤多余的转录本，每个基因提取一个具有代表性的转录本。在这种情况下，我们只保留了每个基因最长的转录本。对最终筛选的转录本进行下游分析。为了确定每个组合的转录组完整性，Benchmarking Universal Single-Copy Orthologs工具(BUSCO, v. 4.0.6) [77]用于获取胚胎植物数据库中单拷贝同源的百分比，也用于评估转录本组装的完整性。

正交线预测和函数注释

使用TransDecoder v. 5.5.0从每个过滤后的汇编转录本中预测开放阅读框(orf) [73］．这个管道主要包括两个步骤。在第一步中，预测每个转录本最长的ORF，截断最小长度为100个氨基酸。然后，对预测的orf进行扫描，使用BLASTP v. 2.5.0搜索，以截断e值1e-10对策划的蛋白质数据库查找同源谱r .成为与…密切相关的物种r . sanguineum，从存入GigaDB的全基因组测序项目下载[78］．所有命中最好的编码肽被保留用于最终的氨基酸序列预测。当一个文本中有多个预测时，我们为每个文本选择得分最高的ORF。我们使用CD-HIT v. 4.7.0 [75]进一步减少最终预测氨基酸的冗余，序列一致性阈值设置为0.95。用OrthoVenn2对3个品种的蛋白质序列进行同源簇(正交群)鉴定[43，使用e-value阈值1e-10。用三个品种注释ORF数据集之间的一对一单拷贝正交关系进行后续分析。的同源蛋白序列R。var。didymoides使用NCBI非冗余(NR)和UniProtKB/Swiss-Prot with BLASTP v. 2.5.0作为蛋白质数据库的代理搜索，设置e-值截止值为1e-10。我们还使用DIAMOND[执行了额外的函数注释]79命中eggNOG数据库[80]，总结了不同蛋白质数据库中可用的功能信息，包括GO、COGs/KOGs和KEGG。最好的结果被用作最终的注释。

阅读映射，SNP调用和集群

对于每个品种，使用HISAT2 v. 2.1.0将三个复制植物的清洁reads(花芽和叶芽的组合)分别与参考转录组(同源簇)对齐[74参数设置为默认值。SAMtools v. 1.9 (https://github.com/samtools/)和Picard tools v. 2.21.8 (http://broadinstitute.github.io/picard/)用于对重复读取进行排序、标记和删除，并对9个样本中的每个样本的bam对齐结果进行重新排序。我们使用GATK v. 4.1.5 [81]执行SNP呼叫和Plink v. 1.9 (http://pngu.mgh.harvard.edu/purcell/plink/)使用设置的参数筛选snp——基因族群0.1——加0.01——indep-pairwise10、0.2。最后，将符合条件的snp合并到单个VCF文件中，作为SNPhylo v. 20180901的输入[82]，以最大似然重建样本之间的系统发育关系，使用1000个自举复制进行分支支持。为了进一步研究遗传变异的分布，本文采用主成分分析(PCA)方法对黄曲霉的SNP变异进行了分析R(69］．

转录丰度和差异表达分析

基因表达水平(丰度估计)通过使用RSEM v. 1.3.1将每个生物复制的花芽样本的所有配对端reads分别映射回每个品种的参考转录组的一对一单拷贝正交图来计算[83]，采用Bowtie2 v. 2.3.5进行对齐[84］．因为每个品种集合都是独立从头生成的，没有参考，所以Trinity ID头被随机分配给每个品种。为了确保每个量化文件被分配相同的ID头，从而可以集成，我们替换了Trinity生成的抄本ID，用其各自的单一副本对应名称(对应的集群)。在获得每个生物复制花芽样本的基因表达丰度后，我们基于每千碱基的外显子每百万片段映射reads (FPKM)的片段生成了基因表达矩阵。用Trinity脚本生成的读取计数矩阵用于差异表达基因(DEGs)分析。用DESeq2包进行差异分析[85在R以两两比较的方法鉴定差异基因。deg被认为是那些调整了错误发现率(FDR)的p取值≤0.05,log2(倍数变化)的绝对值≥1。为了比较三个品种的基因表达值，我们使用了m值归一化的修剪平均值(TMM)，如在R包磨边机(86］．所有下游分析均基于归一化表达式数据矩阵(TMM归一化)实现。从各品种间的比较中筛选出deg进行进一步的功能富集分析。在对所有同源基因进行功能注释的基础上，通过内部Perl脚本提取所有注释项，并导入到一个R包AnnotationForge生成一个OrgDB生物体注释对象，该对象包含上述数据库中基因ID和其他标识符之间的对应映射。差异表达基因(DEGs)的GO(基因本体)和KEGG(京都基因与基因组百科全书)富集分析由clusterProfiler包实现[87在R．

表达谱的统计分析

所有的表达分析分别对每个品种的三个生物重复进行。构建了9列15164行的基因表达矩阵。每一列代表一个样本，每一行对应一个同源基因的表达。用数据矩阵计算Pearson相关系数(r)。对所有样品进行了对称热图和主成分分析RPheatmap v. 1.0.12 (https://CRAN.R-project.org/package=pheatmap)和PCAtools v. 1.1.0 (https://github.com/kevinblighe/PCAtools),分别。GO术语和代谢途径p≤0.05的细胞被DEGs显著富集。

实时定量PCR (qRT-PCR)验证

为了验证RNA-seq结果与花青素生物合成途径的可靠性，选取了13个基因进行qRT-PCR分析。使用BioRT Master HiSensi cDNA First Stand Synthesis Kit (Bioer, Hangzhou, China)，稀释20倍作为模板合成第一链cDNA文库。在QuantStudio™7 Flex实时PCR系统(Applied Biosystems, CA, US)上，使用BioEasy主混合SYBR Green (Bioer，杭州，中国)对3个生物重复(植物)进行QRT-PCR。qRT-PCR扩增条件为:95°C变性1 min, 95°C变性15 s, 60°C退火延伸60 60s，循环40次。12个基因的引物列于表中S7．GLYCERALDEHYDE-3-PHOSPHATE脱氢酶（GAPDH)r . simsii(基因库acc。不。FN552706)作为规范化的内部控制[50]，相对表达量用2^−ΔΔCT方法(88］．

启动子序列多样性的鉴定与表征

进一步了解表达水平是否受独联体-调控元件和启动子的DNA甲基化，根据已发表的cDNA序列进行了检索杜鹃基因组(78]，以序列相似度最高的最佳命中基因作为相应参考。然后我们将每个基因组测序品种(三个个体组合)的干净读数映射到每个杜鹃引用序列，并通过sam2consensus.py脚本(https://github.com/edgardomortiz/sam2consensus)．生成每个基因的每个品种的一致性序列(fasta文件)。翻译起始位点上游的每个2kb启动子序列从一致序列和中缩短独联体-作用调控元件由PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/)数据库。TBtools v. 1.068 [89，然后用来可视化的分布独联体表演元素。利用在线平台mVISTA将三个品种中每个基因的启动子序列成对比对到每个参考序列[90]，并在50 bp的滑动窗口内绘制出基对识别图，范围为75 ~ 100%。滑动窗口为50 bp，计算各品种各基因启动子的GC含量，并用ggplot2包进行可视化R(69］．

原始数据可从NCBI序列读取档案(SRA)平台访问，登录号为PRJNA720238 (http://www.ncbi.nlm.nih.gov/bioproject/720238)．

4 cl:

4-coumarate-CoA连接酶

ANR:

花青素还原酶;

答:

花青素合成酶

:

乙酰转移酶

BGLU:

Beta-glucosidase

C4H:

肉桂acid4-hydroxylase

计算机辅助设计:

肉桂醇脱氢酶

气:

查耳酮异构酶

CHS:

查耳酮合酶

DFR:

Dihydroflavonol-4-reductase

F3'H:

黄烷酮3 ' -羟化酶

F3H:

黄烷酮3-hydroxylase

读者:

黄酮醇合成酶

销售税:

谷胱甘肽s-transferase

MBW:

MYB-bHLH-WD40复杂

OMT:

O-methyltransferase

朋友:

苯丙氨酸ammonia-lyase

/:

过氧化物酶

GAPDH:

甘油醛3磷酸脱氢酶

RNA-seq:

RNA序列

度:

差异表达基因

爆炸:

基本的本地对齐搜索工具

车身:

对标通用单拷贝矫形仪

NR:

NCBI冗余

UniProtKB:

通用的蛋白质知识库

齿轮:

同源基团的簇

蛋酒:

基因的进化谱系:非监督的同源组

GATK:

基因组分析工具箱

走:

基因本体论

KOG:

真核蛋白质的同源基团

KEGG:

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

罗斯福:

错误发现率

FPKM:

每kb每百万片段

主成分分析:

主成分分析

SNP:

单核苷酸多态性

TMM:

m值的修剪平均值。

我们感谢博士们。钟美才、蒋晓东等同事(中国科学院昆明植物研究所)和Kanae Nishii博士(神奈川大学和爱丁堡皇家植物园，RBGE)对数据分析和早期版本的手稿提出了建设性的建议。实验室工作和数据分析在中国科学院昆明植物研究所西南野生物种种质库分子生物学实验室和iFlora高性能计算中心进行。RBGE由苏格兰政府的农村和环境科学和分析服务部支持。

本研究得到国家自然科学基金项目(91631101,31670213)、中国科学院战略重点研究项目(XDB31000000)、中国科学院大型科学设施项目(2017-LSFGBOWS-02)、中国科学院国际合作计划项目(151853KYSB20190027)和云南省科技人才培养计划项目(2017HA014)资助。

从属关系

1. 中国科学院昆明植物研究所中国科学院东亚植物多样性与生物地理重点实验室，云南昆明650201

   叶林江，罗亚黄，邹佳云，郑伟，刘杰，胡金勇，高连明

2. 中国科学院昆明植物研究所野生物种种质资源库，云南昆明650201

   叶林江，朱安丹，李德洙

3. 云南大学生命科学学院，云南昆明650091

   叶林江，王跃华

4. 中国科学院大学，北京10049

   叶林江，邹佳云，郑伟，李德柱

5. 英国爱丁堡皇家植物园，EH3 5LR

   迈克尔·Mӧ噢

6. 中国科学院昆明植物研究所云南丽江森林生态系统国家观测研究站，云南丽江674100

   Lian-Ming高

贡献

LMG和DZL构思和设计了这项研究。LJY, YHL, JYZ和WZ进行了实地工作和样品采集。LJY和JYZ进行了实验。LJY分析了这些数据。YHW, JL, ADZ, JYH和LMG引导数据分析。LJY, MM和LMG在其他作者的输入下撰写了手稿。所有作者阅读并批准了最终稿件。

相应的作者

对应到De-Zhu李或Lian-Ming高．

伦理批准和同意参与

本研究所用植物均采自中国国家级自然保护区高黎贡山。对植物的实地研究符合国家指导方针和立法，不涉及任何濒危或受保护物种。

同意出版

不适用。

相互竞争的利益

作者声明他们没有竞争利益。

出版商的注意

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

开放获取本文遵循创作共用署名4.0国际许可协议(Creative Commons Attribution 4.0 International License)，该协议允许在任何媒体或格式中使用、分享、改编、分发和复制，只要您给予原作者和来源适当的署名，提供创作共用许可协议的链接，并说明是否有更改。本文中的图片或其他第三方材料包含在文章的创作共用许可中，除非在材料的信用额度中另有说明。如果材料不包含在文章的创作共用许可中，并且您的预期用途不被法律法规允许或超出了允许的用途，您将需要直接从版权所有者那里获得许可。欲查看此许可证的副本，请访问http://creativecommons.org/licenses/by/4.0/．创作共用公共领域奉献放弃书(http://creativecommons.org/publicdomain/zero/1.0/)适用于本文提供的数据，除非在数据的信用额度中另有说明。

再版和权限