《Frontiers in Plant Science》：Chromatin accessibility directly governs flavonoid biosynthesis and indirectly orchestrates cannabinoid production in Cannabis
编辑推荐：
本文通过整合代谢组学、转录组学和ATAC-seq技术，揭示了大麻花序中染色质开放性通过直接调控黄酮生物合成基因表达，并间接通过影响脂肪酸合成、腺毛发育（茉莉酸信号、花器官特性）等相关基因来协调大麻素合成的表观遗传网络，为通过表观遗传改良提升大麻药用价值提供了新策略。
3.1 两个大麻品种花序的代谢组学分析
大麻富含多种代谢物，不同器官间代谢物的含量和种类存在差异。雌性花序是大麻素生物合成的主要器官。更具体地说，雌性花序中的苞片和萼片覆盖有大量的具柄腺毛。腺毛由分泌盘细胞和分泌腔组成，次级代谢物（包括大麻素和萜类）在此产生和储存。本研究选择了两个形态差异显著的品种。‘H’品种的花序密被腺毛，而‘O’品种的花序仅有少量腺毛。
鉴于两个品种间的显著形态差异，我们推测其花序的代谢物谱存在显著差异。因此，我们利用非靶向代谢组学技术对两个品种的代谢物谱进行了表征。总共鉴定出2153种代谢物，包括921种脂质和类脂分子、250种有机杂环化合物、165种有机酸及其衍生物、147种苯丙素和聚酮化合物以及145种苯环型化合物。代谢组数据中也检测到大麻色素、Δ8-四氢大麻酚等大麻素。大麻花序是一个由不同发育阶段的叶片和花朵组成的异质性集合体，代谢物高度丰富。在H品种中，苞片和叶表面密布着积累广谱次级代谢物的腺毛。同时，花序中含有不同发育阶段的种子；这些种子本身富含脂质和黄烷酮等次级代谢物。这种源自种子的脂质库很可能是我们代谢组数据集中检测到多达921种不同脂质物种的原因。随后，将代谢物注释到KEGG通路。鉴定出甘油磷脂代谢（71）、花生四烯酸代谢（60）、α-亚麻酸代谢（58）和亚麻酸代谢（58）等通路的代谢物。
主成分分析（PCA）显示，根据品种不同，样本明显分离。PC1解释了45.33%的方差，PC2解释了10.36%的方差。这两个主成分共同解释了55.69%的方差，表明这两个主成分能够很好地代表原始数据的变异性。
3.2 两个品种间差异积累的代谢物
接着我们比较了两个品种间的差异积累代谢物（DAMs）。总共有250种代谢物在H品种中高积累，而241种代谢物在O品种中高积累。在Top30代谢物中，黄酮类化合物的含量显示出显著差异。例如，山奈酚和2'-羟基染料木素8-C-葡萄糖苷等黄酮类化合物在H品种中的含量远高于O品种。此外， diketene，一种可以参与黄酮类、生物碱和其他次级代谢物合成的重要合成前体，在H品种中也显著富集。
相比之下，O品种也富含一些黄酮类代谢物，如Tectoridin和Isopongaflavone。此外，膜脂代谢物如PC (O-16:1/18:0)和PC (O-18:1(11Z)/16:0)在O品种中富集。其中，大量的磷脂酰胆碱（PC）化合物，如PC (18:1(11Z)/18:2 (9Z,12Z))、PC (18:0/22:5)和PC (18:0/20:1)，在O品种中富集。这些化合物是细胞膜的重要组成部分。不出所料，代谢组学分析发现大麻色素和Δ8-四氢大麻酚在H品种中的含量更高。KEGG分析显示，DAMs主要富集在亚油酸代谢（ko00591）、甘油磷脂代谢（ko00564）和α-亚麻酸代谢（ko00592）通路。这些结果表明，两个品种花序的代谢物存在显著差异，尤其是在黄酮类和脂肪酸方面。
3.3 利用RNA-seq进行大麻转录组测序
为了进一步探索次级代谢物生物合成的机制，我们对两个品种进行了比较转录组分析。每个生物学重复平均产生4200万条reads。经过质量控制，每个重复平均剩余4100万条reads。超过75%的clean reads成功比对到参考基因组。PCA显示两个品种彼此分离，但重复样本聚集在一起，表明我们的转录组数据是可靠的。
为了验证转录组数据的可靠性，我们选择了七个基因，包括黄烷酮3-羟化酶（LOC115702709）、二氢黄酮醇4-还原酶（LOC115710150）、柚皮素查尔酮合成酶（LOC115724170）、脂肪酸去饱和酶（LOC115705664）、叶绿素A-B结合蛋白（LOC115702322）和铁氧还蛋白–NADP还原酶（LOC115702780）进行qPCR验证。RT-PCR数据显示与RNA-seq相似表达模式，这表明我们的RNA-seq数据是可靠的。
然后我们比较了两个品种间的差异表达基因（DEGs）。O品种中有4926个基因上调，而H品种中有3417个基因上调。随后我们将DEGs注释到GO和KEGG数据库。结果显示，两个品种在不同通路中的DEGs数量存在显著差异，表明两个品种在转录水平上具有不同的基因表达调控模式。在H品种中，DNA模板转录（GO:0006351）、DNA模板转录调控（GO:0006355）、蛋白结合（GO:0005515）和DNA结合转录因子活性（GO:0003700）等通路，以及与黄酮生物合成、甘油磷脂代谢、脂肪酸生物合成和糖异生/糖酵解相关的通路中，上调基因的数量显著高于O品种。相比之下，O品种在蛋白磷酸化（GO:0006468）和蛋白丝氨酸/苏氨酸激酶活性（GO:0004674）等项目，以及卟啉代谢、光合作用、植物昼夜节律和ABC转运蛋白等通路中，上调基因的数量更多。KEGG富集分析显示，DEGs主要富集在光合作用（map00195）、黄酮生物合成（map00941）、MAPK信号通路–植物（map04016）等通路。黄酮生物合成通路中的几个关键基因，例如查尔酮异构酶样蛋白（LOC115706946）、黄酮醇合成酶（LOC115708857）、二氢黄酮醇4-还原酶（LOC115710150）和黄烷酮3-羟化酶（LOC115702709）在H品种中上调。细胞色素P450（CYP）是一个广泛分布于植物中的单加氧酶超家族。这些酶在黄酮类化合物的生物合成中起着至关重要的作用，并显著影响黄酮类化合物的多样性。细胞色素P450（CYP450）是植物中最大的酶超家族之一，参与多种次级代谢物的生物合成，决定了其结构多样性和生物活性。几个P450家族成员（LOC115709845和LOC115709933）在H品种中上调。
基因集富集分析（GSEA）也显示，植物器官形态发生相关条目，如种子发育和花器官特性指定，在H品种中上调。几个与花和种子发育相关的转录因子，如MADS-box（LOC115714657）、B3（LOC115696465）和YABBY（LOC115704942, LOC115704995, LOC115715807）基因在H品种中上调。贮藏蛋白，如大麻球蛋白和细胞周期蛋白，在H品种中也上调。腺毛主要分布在雌花的萼片和花序的幼叶上。这种分布模式表明，花形态发生可能是腺毛和次级代谢物合成的必要前提。在H品种中，大量花发育相关基因表现出高表达水平，表明该品种已进入生殖生长阶段，进而启动了大量次级代谢物（如大麻素）的生物合成。此外，大麻种子富含脂质和蛋白质。在H品种的种子成熟过程中，脂质含量的显著增加是必然结果。大麻素生物合成通路基因，包括OAC（LOC115723438）、TKS（LOC115699293）、PT（LOC115713215）、CBDAS（LOC115697762）在H品种中也表现出高表达。
与染色质组织和染色体结构相关的基因也在H品种中富集。几个染色体结构因子，如染色体凝缩调节因子（RCC1）家族蛋白（LOC115722105和LOC115709118）和染色体结构维持蛋白（LOC115704761）在H品种中高表达。这些染色体相关因子的上调表明细胞可能在染色质结构调控、基因表达调控和细胞周期等生物过程中发生了显著变化。
而光合作用相关基因在O品种中上调。例如，细胞色素b6-f复合体铁硫亚基2（LOC115725694）、光系统II PsbY（LOC115714330）、光系统I反应中心亚基N（LOC115712566, LOC115719040）、psbP样蛋白1（LOC115705010）这些光合作用的核心蛋白在O品种中上调。这些结果表明，光合作用仍然是O品种的主要生物过程。因此，腺毛的数量相对较少，进而导致大麻素等次级代谢物含量降低。由于政府严格的禁毒政策，例如持续清除高密度腺毛品种，只有低密度腺毛品种被保留下来并允许种植。转录组数据表明两个品种之间存在显著差异，导致其代谢物含量存在显著差异。
3.4 染色质开放性分析
代谢组学和转录组学分析揭示了两个品种在代谢物和基因表达方面存在广泛差异，尤其是次级代谢物生物合成通路基因和代谢物。我们假设表观遗传调控导致了关键通路基因的表达差异。一致的是，转录组数据也显示几个染色质相关基因的表达存在显著的品种特异性差异。为了验证这一点，我们利用ATAC-seq系统比较了两个品种花序的染色质开放性。对六个样本（每个品种三个生物学重复）进行了测序，每个样本平均产生1.75-1.8亿条原始reads。经过质量过滤后，约98-99%的高质量reads被保留，证明了测序数据的优异质量和可靠性。clean reads的平均Q30为93.93%。
鉴定出的可及区域主要富集在基因组转录起始位点（TSS）附近，表明ATAC-seq数据质量很高。H品种三个生物学重复之间的相关性超过0.95，而O品种三个生物学重复之间的相关性超过0.79。而不同品种之间的相关性仅为0.56。这些结果表明两个品种间的染色质开放性存在差异，且数据显示出可重复性和可靠性。
随后，我们将peaks注释到基因。结果显示，绝大多数peaks注释在远端基因间区（H品种约50%，O品种约73%）。启动子在调控基因表达中起着至关重要的作用。大量peaks注释在启动子区域（H品种约37%，O品种约20%），并且启动子的近端（<1kb）区域比远端区域富集了更多的peaks。启动子区域的这些peaks可能影响启动子的驱动力。少量peaks注释在UTR和基因体区域（内含子、外显子和下游）。
DARs及相关基因（不同可及性基因，DAGs）详见补充材料。总共有400个基因在O品种中显示出高可及性，而10976个基因在H品种中显示出高可及性。GO富集分析显示，胞质大核糖体亚基（GO:0022625）、核糖体结构成分（GO:0003735）和RNA修饰（GO:0009451）相关基因在两个品种间显示出可及性差异。KEGG注释显示，许多代谢通路基因（淀粉和蔗糖代谢（map00500）和苯丙烷生物合成（map00940））在两个品种间表现出染色质可及性差异，这可能导致通路内基因表达的变化，进而引起代谢物水平的差异。例如，属于黄酮生物合成通路的柚皮素查尔酮合成酶（LOC115724170和LOC115712997）和细胞色素P450（LOC115709845, LOC115709933）在H品种的启动子区显示出高可及性。转录组分析显示，两个品种间的DEGs富集在黄酮生物合成通路。我们的结果表明这些基因的染色质可及性可能导致它们在H品种中高表达。高度参与植物组织中各种碳水化合物化合物代谢的糖苷水解酶在启动子区显示出高可及性。这些DARs可能导致两个品种间淀粉和糖含量的差异。苯丙烷生物合成通路是大麻素合成的上游通路之一。苯丙烷类化合物（如莽草酸、肉桂酸等）的代谢物是合成大麻素的关键前体物质。苯丙氨酸通过苯丙氨酸解氨酶（PAL）转化为肉桂酸。大麻的PAL同源基因（LOC115709862）在H品种的启动子区显示出更高的可及性。这可能是H品种大麻素含量增加的原因之一。
DARs也富集在植物激素信号转导（map04075）、MAPK信号通路–植物（map04016）中，表明两个品种在信号通路转导上存在差异。双组分信号系统元件基因：组氨酸-containing phosphotransferase protein（LOC115705353）和调节因子（LOC115709125）的启动子在H品种中显示出高可及性。脱落酸不敏感5（ABI5, LOC115708807）在H品种的启动子区也比O品种具有更高的可及性。而许多转录因子（碱性螺旋-环-螺旋、GRAS、AP2）的启动子在O品种中显示出高可及性。这些结果表明两个品种在激素信号转导过程中存在表观遗传差异。
3.5 代谢组、转录组和ATAC-seq数据的整合
为了更全面、综合地分析两个品种间的差异，我们进行了整合的转录组-代谢组和转录组-ATAC-seq分析。总共有250种代谢物和4926个基因在H品种中高度富集，而在O品种中，有241种代谢物和3417个基因高度积累。九象限图将差异表达基因和代谢物划分为促进代谢物合成的基因集和抑制代谢物合成的基因集。α-亚麻酸代谢、ABC转运蛋白、光合生物中的碳固定以及黄酮生物合成等通路均存在差异表达基因和代谢物。
对8343个DEGs和11096个DAGs进行整合分析，揭示了6045个基因的染色质可及性水平与其转录丰度显著相关。鉴于ATAC-seq在O品种中检测到的DAGs数量明显较少，我们重点关注了H品种中染色质可及性增加但表达模式出现分歧的基因。具体来说，3361个基因（Hyper-up）在其基因体内表现出增强的染色质可及性，并且在H品种中转录上调。相反，2416个基因（Hyper-down）表现出升高的染色质可及性，但在H品种中转录下调。随后的GO和KEGG富集分析显示，Hyper-up基因显著富集于植物激素信号转导（ko04075）、脂肪酸生物合成（ko00061）和黄酮生物合成（ko00941）等通路。这些基因在GO条目中也过度代表，包括DNA模板转录（GO:0006351）、DNA结合转录因子活性（GO:0003700）和转录顺式调控区域结合（GO:0044212）。相比之下，Hyper-down基因优先富集于光合作用（ko00195）、昼夜节律–植物（ko04712）和α-亚麻酸代谢（ko00592）等通路，以及对几丁质反应（GO:0010200）、对伤害反应（GO:0009611）和DNA结合转录因子活性（GO:0003700）。我们的结果表明，除了与大麻素生物合成相关的基因和与腺毛起始/发育相关的基因外，其他基因如花转变和调控染色质开放的基因也可能影响大麻素的含量。
我们随后对两个大麻品种的黄酮和大麻素生物合成通路进行了系统评估。代谢组学显示，山奈酚、异荭草苷和槲皮素-3-O-葡萄糖苷等黄酮类化合物在H品种中的积累水平显著高于O品种。转录组分析进一步显示，关键的黄酮生物合成基因，包括黄酮3′,5′-甲基转移酶（LOC115712973）、查尔酮异构酶（CHI, LOC115706946）、查尔酮合成酶（LOC115724170）和黄酮醇合成酶（LOC115708857）在H品种中显著上调。ATAC-seq数据显示，这些基因转录起始位点上游2 kb内的染色质可及性在H品种中显著增加。启动子区域染色质可及性的增加使得转录因子、增强子和其他调控元件能够结合启动子 motif，从而激活或抑制下游基因的表达。总之，这些发现表明，H品种关键黄酮生物合成基因启动子的开放染色质状态增强了其转录活性，从而驱动了高效的黄酮积累。
对两个品种大麻素生物合成通路的系统比较显示，大麻色素、Δ8-四氢大麻酚和1-羟基四氢大麻酚等大麻素在H品种中的积累水平高于O品种。转录组分析显示，关键生物合成酶基因——OAC（LOC115723438）、TKS（LOC115699293）和CBDAS（LOC115697762）——在H品种中显著上调。ATAC-seq显示这些基因的基因体或启动子区域存在peaks，然而，peak强度在两个品种间没有显著差异。这表明染色质可及性有助于大麻素生物合成通路基因的表达，但不太可能是观察到的品种间大麻素含量变异的关键驱动因素。启动子区域更大的染色质可及性暴露了转录因子结合位点，促进了它们的快速占据并驱动基因转录。研究发现了三种能够结合THCAS基因启动子并改变其表达的新型转录因子。启动子上的peaks可能是这种转录调控的前提条件。我们的数据显示，位于OAC（LOC115723438）、TKS（LOC115699293）和CBDAS（LOC115697762）启动子上的peaks可能是未来研究大麻素生物合成调控的关键区域。
大麻素生物合成强烈依赖于脂肪酸衍生的前体。定量结果表明，脂肪酸代谢通路基因，如脂肪酸去饱和酶（LOC115705664）、硬脂酰-9-去饱和酶（LOC115698723）和乙酰-CoA载体（LOC115720238）在H品种中其启动子区的染色质可及性增加，导致表达显著升高。这表明染色质可及性通过影响大麻素前体的生物合成基因来改变大麻素的含量。同时，这些脂肪酸也是大麻种子的重要组成部分，表明大麻种子的发育可能也受到表观遗传的调控。大麻素在腺毛中合成和储存；因此，腺毛密度直接决定大麻素含量。腺毛发育在番茄、烟草和拟南芥中得到了广泛研究。SlMXI、NtMYC2a/b、GLABRA1（GL1）和GLABRA2（GL2）等基因被认为是腺毛发育的关键调节因子。核心腺毛调节因子GLABRA 2（LOC115707829）在H品种中显著上调，并表现出升高的启动子可及性。这可能是两个品种腺毛数量差异的关键原因。先前报道表明，外源茉莉酸甲酯（MeJA）会增加腺毛密度。我们的数据显示，关键的MeJA响应基因，如JAZ（LOC115723116）和MYC（LOC115704602）在H品种中也高表达，并且其启动子的染色质开放性显著增强。腺毛主要富集在雌花的苞片和花序的叶片上，将腺毛发育与开花紧密联系起来。在H品种中，开花整合因子FLOWERING LOCUS T（LOC115700781）和花特性转录因子MYB113（LOC115708959）、MADS-box（LOC115716337）和YABBY（LOC115704995）均显著上调，并伴随着其启动子区染色质可及性的增加。
本研究从‘染色质开放-转录-代谢’的角度探讨了工业大麻次级代谢物合成的调控机制。然而，值得注意的是，两个品种在基因组水平上的差异可能是它们之间差异的重要原因。例如，大麻和hemp的基因组中存在大量SNPs，这导致两个物种间基因表达的差异，并且这些差异不仅限于与大麻素生物合成相关的基因如THCAS。研究发现THCAS和CBDAS包含大的富含逆转录转座子的区域，这些区域在drug型和hemp型等位基因之间高度非同源。因此，未来比较O和H品种基因组的差异可能会揭示更多影响大麻素合成的因素。
总之，大麻素生物合成通路基因的染色质可及性并非两个品种大麻素含量差异的原因。相反，H品种中大麻素积累的增加似乎主要由前体供应增强（上调的脂肪酸代谢）和腺毛丰度增加所驱动。相比之下，前体合成、腺毛发育、茉莉酸信号和开花控制相关基因的表达增强与表观遗传激活密切相关，表现为其各自启动子区染色质可及性的升高。
