沙棘刺和芽的组织解剖与转录组分析

王建武，杨 涛，严加坤，相微微，冯光惠，刘翠英，尚爱军*

(1 榆林学院陕西省陕北矿区生态修复重点实验室, 陕西榆林 719000；2 陕西省林业科学院，西安 710082；3 国家林业和草原局长柄扁桃工程技术研究中心，陕西榆林 719000)

自然界中的许多植物都有刺，但不同植物中刺的发育机制存在显著差异。英语词汇里有专门的单词来区分不同的刺，即prickle(由皮层发育而来)、spine(叶子的尖锐部分)和thorn(尖锐的树枝)[1]，玫瑰刺是植物表皮毛发育而形成的尖锐突起，通常被误认为是真正的叶刺或枝刺[2]。本研究在讨论尖锐的植物器官时统一使用“刺”一词。刺的明显功能是防御，刺可以提供机械保护，防止草食动物啃食；刺可以帮助种子和植物扩散；刺可以增加表皮的厚度，从而减少热量和水分的散失；一些刺可以防止阳光直射，降低光抑制，并可降低低温产生的损伤[1-2]。

从生理角度来看，植物刺是有重要功能的器官，刺的发育和形成除了受自身基因控制外，环境因素对刺的发育和形成也有重要影响。动物啃食可以刺激植物产生更多的刺，以防御动物过度啃食[1]；土壤湿度降低会使黑枸杞的茎刺变得坚硬[3]；遮阴可以使鹰爪化的钩刺数量增多，有利于攀缘[4]。因此从遗传的角度来看，植物刺是一个数量性状，而不是一或几个基因控制的质量性状。刺的发育和形成由多个基因控制，其分子机理应该是一个复杂的调控网络。目前，研究人员把表皮毛作为研究拟南芥和棉花细胞增殖和生长发育的模式细胞，因此，由表皮毛发育成刺的分子机理相对清晰[2,5-9]。对于叶刺、皮刺、枝刺等多细胞构成的刺，目前已通过遗传和分子生物学手段克隆鉴定到一些相关基因，如玫瑰的RrTT-G1基因[10]、花椒的10个皮刺分化相关基因[11]、莴苣的LG5467628基因(混池RNA-seq测序后定位到的候选基因)[12]及生菜(做了基因定位分析，尚未具体到基因水平)[13]，但尚未建立一个调控网络，而且已有的研究相互印证和支持性也不强，也进一步说明刺的发育和形成分子机制比较复杂。不同植物的刺，甚至同种植物不同部位的刺，其发育和形成分子机理可能都有其特异性。

沙棘是多年生灌木，根具有固氮根瘤菌，耐贫瘠和干旱[14]，是优良的防风固沙植物，具有重要的生态价值。此外，沙棘果实具有很高的营养和药用价值，可以加工沙棘汁，并从果实中提取药用成分[15-17]。沙棘枝上有很多刺，沙棘刺含有维管组织，折断时断面粗糙，属于典型的枝刺(thorn)。从外部形态来看，沙棘刺具有营养叶，发育早期的刺和枝芽差异很小，只是刺顶端生长点消失，而且刺的着生位置跟枝芽也没有明显差异，总体上刺和芽在茎上有规律地间隔开。这些特征跟黑枸杞和酸枣的枝刺有较大区别，黑枸杞和酸枣的刺不具有营养叶，而且长在叶腋中，因此，对沙棘刺的研究可以聚焦于比较分析去掉营养叶的刺和芽之间的差异，在取材上能尽量去掉背景噪音，凸显差异。从采摘的角度来看，人们希望获得少刺或无刺的沙棘品种，然而要通过遗传手段调控刺的发育，培育少刺或无刺新品种，需要对沙棘刺的发育和形成分子机理有深入了解。综上，对沙棘刺的研究既有理论意义，也有应用前景。本研究比较了沙棘刺和芽的形态和组织差异，通过转录组测序分析比较了沙棘刺和芽基因表达谱，挖掘到与木质素、细胞壁、表皮和气孔相关的差异表达基因(differentially expressed genes，DEGs)，并通过qRT-PCR验证了关键基因的表达水平。

1 材料和方法

1.1 植物材料

从7年生雄性沙棘(HippophaerhamnoidesL.ssp.sinenisRousi)上采集芽和刺，该沙棘品种多刺。种植地点位于陕西省榆林市榆阳区陕西省防沙研究院苗圃内(海拔1 129 m；N38°19′49″，E109°42′42″)。在2021年4月22日，选择约2 cm长的芽和刺，在开花前采集，除去营养叶，液氮中速冻，干冰运输到上海美吉公司进行转录组测序。样品中所有4个生物重复都来自同一株植株，芽标记为bud1、bud2、bud3、bud4，剌标记为thorn1、thorn2、thorn3和thorn4。

1.2 组织切片分析

样品在FAA中固定过夜，将固定的材料分割成0.3 cm大小，梯度乙醇脱水，石蜡包埋，Leica RM2135 石蜡切片机切片，切片厚度 8～12 μm，然后二甲苯脱蜡，梯度乙醇复水，番红-固绿染色，中性树胶封片，Leica DMLB显微镜观察并拍照。切片用开放场或在340～380 nm激发后，用Leica DMLB荧光显微镜观察，并用摄像机(Leica，DC 300F)记录[18]。

1.3 RNA-Seq分析

使用TRIzol试剂(美国Invitrogen)从样品中提取总RNA，并在Illumina HiSeq平台上进行RNA测序。RNA测序数据分析按标准化流程进行[19-21]。每个转录本的表达水平被标准化为每百万个转录本(transcripts per million，TPM)值[22]。在本研究中，每个转录本的差异表达水平以| log2(fold-change)| ≥ 1和P-value ≤ 0.05为标准，以确定芽和刺之间的DEGs。使用BLASTx程序，针对NCBI NR(non-redundant protein)、Swiss-Prot、Pfam、KEGG(Kyoto Encyclopedia of Genes and Genomes)、KOG(euKaryotic orthologous groups of proteins)和GO(gene ontology)数据库进行比对，对unigenes进行注释[20]。基因聚类分析采用MEGA 7.0软件。

1.4 qRT-PCR验证

所用样本与转录组分析相同。qRT-PCR在FTC-3000仪器(加拿大，Funglyn生物技术公司)上使用TB-Green方法(TB-Green®Fast qPCR Mix，日本Takara)进行。用于qRT-PCR的引物为unigene DN7364_c0_g2、DN4324_c0_g3、DN10443_c0_g2、DN8484_c0_g2、DN16315_c0_g1、DN8647_c0_g3、DN424_c0_g1、DN8861_c0_g2、DN14231_c0_g1、DN7861_c0_g1、DN11032_c0_g1、DN13014_c0_g1、DN16737_c0_g1和DN13978_c0_g1(表1)。内参基因采用通过转录组测序得到的沙棘ef1α基因片段(unigene DN5364_c0_g6)(表1)[23]。Unigene DN5364_c0_g6被较好注释，为延伸因子ef1-α基因，包含466 bp。bud1、bud2、bud3、bud4、thorn1、thorn2、thorn3和thorn4中DN5364_c0_g6的TPM值分别为68.09、81.53、88.84、66.05、89.7、70.73、84.69和79.32，该基因表达相对稳定，适合作为内参基因[24]。

表1 qRT-PCR引物Table 1 Primers for qRT-PCR

1.5 统计分析

qRT-PCR试验采用3个生物学重复，数据均为“平均值±标准偏差”。统计分析用SPSS 22.0。

2 结果与分析

2.1 沙棘刺和芽的形态和组织比较

成熟的沙棘刺长约3 cm，具有营养叶；同一年生长的刺会在秋天时从刺尖开始枯萎(图1，A)。总的来说，沙棘芽横切面的组织边界明显且层次分明；而沙棘刺横切面的组织边界模糊，染色区域类型单调，且刺尖中心的红色区域较大(图1，C～F)。这可能是因为刺尖细胞开始变得像导管一样，原生质体逐渐解体并死亡，只留下细胞壁。木质素的自发蓝色荧光显示刺中的木质化程度显著高于芽，木质化程度的增加可增加刺的硬度(图1，G～J)。

2.2 转录组数据整体分析

转录组数据优化组装结果(表2)显示，经数据组装共获得81 560个unigenes(单基因)，其中60.22%的unigenes被定位，unigenes和transcript(转录本)的TransRate得分分别为0.303 14和0.418 81。统计表明，一个从头组装的得分为0.22(优化得分为0.35)的TransRate分析结果将优于50%的已发布在NCBI TSA数据库中的分析结果[25]，故本研究的转录组组装质量相对较好。BUSCO得分C:75.0%(S:70.3%；D:4.7%)，也相对较好(表2)[26]。

表2 沙棘转录组数据优化组装结果评估表Table 2 Results of optimized assembly for sea buckthorn transcriptome

组装完成后，通过BLAST在6个公共数据库(NR、Swissprot、KEGG、COG、Pfam和GO)中对组装的unigenes进行注释，其中分别有38 210、31 176、18 728、34 175、30 298和32 173个unigenes在上述数据中得到比对注释。至少在一个数据库得到注释的unigenes数目为42 193(占所有unigenes的52%)。GO和KEGG注释基因在所有单基因中的比例分别为39%和23%。共获得9 385个DEGs(上调6 428个、下调2 957个)。DEGs聚类分析表明，本研究的4个生物学重复具有良好的一致性(图2)。

为了进一步了解DEGs的功能，将所有DEGs映射到GO数据库，并进行富集分析。GO分析表明DEGs被分为三大功能类别47个项目，其中BP类别、CC类别和MF类别中的项目数量分别为20、12和15。图3显示了前20个GO富集程度最高的项目及其DEGs数目。排名前三的GO富集项目分别是植物表皮发育(plant epidermis development)、气孔复合体发育(complex stomatal development)和木质素分解代谢(lignin catabolic process)过程。KEGG富集分析得到DEGs参与的138个代谢途径，图4显示了前20个KEGG富集程度最高的项目及其DEGs数目，其中角质、木栓质和蜡质生物合成途径(cutin, suberine and wax biosynthesis)是排名第一的途径。

2.3 木质素分解代谢过程相关的DEGs

木质素的生物合成和代谢途径主要分为莽草酸途径(shikimic acid pathway)、苯丙烷途径(phenylpropanoid pathway)和木质素特异性合成途径(lignin-specific synthesis pathway)三个阶段。迄今为止，莽草酸途径中参与木质素生物合成的调控机制尚未被确定。然而，在苯丙烷途径和木质素合成特异途径中发现了能够特异性调节木质素合成的关键酶和相关转录因子[27-29]。本研究通过GO富集分析发现了20个与木质素合成有关的DEGs，其中只有2个属于苯丙烷途径，即DN2241_c0_g1(肉桂醇脱氢酶1，Cinnamyl alcohol dehydrogenase 1)和DN8647_c0_g3(咖啡酸3-O-甲基转移酶，Caffeic acid 3-O-methyltransferase)；两种推定的蛋白(DN7364_c0_g1和DN4324_c1_g1)的具体功能尚不清楚；其他16个DEGs具有漆酶(laccase)活性，这些DEGs属于木质素合成特异途径(表3)。有些漆酶能降解木质素，有些漆酶能合成木质素，这与漆酶的来源和合成周期有关。与芽相比，刺中只有5个具有漆酶活性的DEGs表达下调，另外11个具有漆酶活性的DEGs表达上调(表3)。其中，DN8484_c0_g2、DN8484_c0_g3、DN8484_c0_g1、DN27398_c0_g1和DN16315_c0_g1在芽中的表达非常低，几乎在刺中特异表达(表3)。接下来，选择6个基因，通过qRT-PCR验证其表达水平。这些基因的相对表达水平与TPM值一致，与芽相比，DN4324_c0_g3、DN8484_c0_g2和DN16315_c0_g1在刺中的表达分别上调了9.5倍、41.1倍和13.7倍； DN7364_c0_g2、DN10443_c0_g2和DN8647_c0_g3的下调显著，分别只有芽中表达量的5%、14%和24%(图5，A)。

表3 沙棘芽和刺中木质素分解代谢相关DEGsTable 3 DEGs linked with lignin catabolic process of sea buckthorn bud and thorn

2.4 细胞壁、表皮和气孔相关的DEGs

通过KEGG途径富集分析，角质、木栓质和蜡质生物合成是最富集的途径。因此，我们分析了属于该途径的DEGs。共有29个DEGs与角质、木栓质和蜡质生物合成有关；这些DEGs可分为乙醇合成相关的脂肪酰辅酶A还原酶(Alcohol-forming fatty acyl-CoA reductase)、ω-羟基棕榈酸酯O-阿魏酰转移酶(Omega-hydroxypalmitate O-feruloyl transferase)、细胞色素P450(Cytochrome P450)、热头蛋白(HOTHEAD protein)和其他类别(表4)。许多基因在芽和刺中的表达水平具有相反的趋势，要么在刺中表达，而在芽中表达非常低，要么在芽中表达，而在刺中表达非常低(表4)。例如，DN53228_c0_g1、DN15725_c0_g4、DN41003_c0_g1、DN14231_c0_g1和DN16701_c0_g1在芽中的表达水平很低，但在刺中高表达，而DN424_c0_g1、DN48710_c0_g1和DN27861_c0_g1在刺中的表达水平很低，但在芽中高表达；特别是，刺中DN27861_c0_g1的TPM值为0(表4)。接下来，从上述基因中选择4个基因，通过qRT-PCR验证其表达。这些基因的相对表达水平与TPM值一致，与芽相比，刺中的DN8861_c0_g2和DN14231_c0_g1分别上调3.6倍和22.1倍；DN424_C0_g1和DN27861_C0_g1的下调显著，分别只有芽中表达量的3%倍和6%倍(图5，B)。

表4 沙棘芽和刺中角质、木栓质和蜡质合成相关DEGsTable 4 DEGs linked with cutin, suberine and wax biosynthesis of sea buckthorn bud and thorn

植物表皮发育和气孔复合体发育是GO富集分析结果中的前2个项目，因此，本研究进一步分析了相关DEGs。植物表皮发育与气孔复合体发育密切相关。本研究鉴定的参与植物表皮发育和气孔复合体发育的DEGs是相同的；这些DEGs可分为表皮形成因子(Epidermal patterning factor，EPF)、转录因子(Transcription factor)、多聚半乳糖醛酸酶(Polygalacturonase)和其他类别(表5)。与芽相比，所有DEGs在刺中均下调表达，DN55365_c0_g1、DN11032_c0_g1、DN13014_c0_g1、DN29670_c0_g1、DN46942_c0_g1、DN16737_c0_g1、DN13978_c0_g1和DN16781_c0_g1在刺中的表达水平非常低，但在芽中高表达；特别是，刺中DN29670_c0_g1的TPM值为0(表5)。接下来，从上述基因中选择4个基因，并通过qRT-PCR验证其表达。这些基因的相对表达水平与TPM值一致，与芽相比，刺中的DN11032_c0_g1、DN13014_c0_g1、DN16737_c0_g1和DN13978_c0_g1显著下调，分别只有芽中表达量的6%、13%、3%和4%(图5，C)。

表5 沙棘芽和刺中表皮/气孔发育相关DEGsTable 5 DEGs linked with epidermis/stomatal development of sea buckthorn bud and thorn

3 讨 论

前人研究认为，无论侧芽最终是发育成刺还是枝，都是在芽个体发育的早期阶段决定的，刺芽通常比枝芽产生更小、更少的三角形叶原基[30]。刺发育过程中的显著变化是顶端分生组织缩小，最终变硬，同时伴随分生组织活动和叶片起始发育的停止[30]。

顶端分生组织的分子调控机制非常复杂，仅基于转录组数据很难进行深入、明确的分析。最近，人们发现水稻中也存在与拟南芥HTH基因高度同源的HOTHEAD(HTH)基因[31-32]。ONI3在水稻茎尖细胞中特异表达，ONI3突变体在茎尖细胞中表现出发育缺陷，导致茎尖发育异常[33]。类似地，Mini1在水稻幼嫩茎尖表达，但在根中不表达，Mini1突变体芽的发育提前终止[34]。因此，推测水稻HTH基因和拟南芥HTH基因具有不同的功能；水稻HTH基因在维持茎尖分生组织(SAM)活性和促进水稻茎尖发育方面起着重要作用[33-34]。在本研究中，挖掘得到5个HTH基因，其中unigene DN27861_c0_g1在刺中不表达，DN48710_c0_g1在刺中的表达也非常低，DN34353_c0_g1在刺中表达增加，而其他2个基因在刺中表达降低。推测其中一些基因可能参与角质合成，类似于拟南芥HTH基因的功能，一些可能参与分生组织活动。

细胞壁木质化可显著增加刺的硬度。因此，沙棘刺的木质化程度显著高于芽，这与沙棘刺的硬化一致。木质素分解代谢过程是GO富集分析结果中的第3项，表明刺和芽中的木质素代谢存在显著差异。大多数与木质素分解代谢过程有关的DEGs属于漆酶；漆酶负责将氨酚醇聚合成木质素聚合物[35-36]。单一木质素的组成在植物种类、组织和细胞类型之间变化很大，甚至在单个细胞内的独立细胞壁层之间差异也很大，该组成主要影响最终木质素聚合物的化学结构和物理化学性质[35]。漆酶的组织定位以及它们的生化特性决定了木质素聚合的组织模式[35]。沙棘漆酶生化特性的深入研究，可将其作为经基因工程手段调控沙棘刺硬度的候选基因。

在刺的发育过程中，细胞的功能类型逐渐趋于单一化，也就是说，它们变硬并起到机械支撑的作用。因此，与芽相比，刺细胞的细胞壁成分发生了显著变化。本研究结果支持这一观点。首先，组织切片分析表明，刺细胞的染色结果相对单一，层次性不如芽切片明显，这表明刺细胞的功能在形态学上得到了简化。其次，KEGG富集分析表明，角质、木栓质和蜡质生物合成是富集程度最高的途径，乙醇形成相关的脂肪酰辅酶A还原酶、ω-羟基棕榈酸酯O-阿魏酰转移酶、细胞色素P450、HOTHEAD蛋白家族在角质、木栓质和蜡质生物合成途径中有多种异构体，这些蛋白质异构体属于进化中的旁系同源物，应具有明显的功能分化。因此，进一步研究这些同源蛋白的生物化学功能，分析它们参与的代谢途径，对阐明沙棘刺和芽之间角质、木栓质和蜡质生物合成的差异具有重要意义。第三，沙棘刺中表皮和气孔的发育可能逐渐停止，本研究发现参与表皮和气孔发育的所有DEGs在刺中均下调表达。bHLH转录因子SPEECHLESS(SPCH)在气孔发育中扮演关键角色[37-39]，SPCH整合了内源信号，并作为激素和环境信号的受体，在发育过程中调节气孔密度和发育模式[37-39]。SPCH的表达水平受分泌型的表皮形成因子(EPF)调节，SPCH与EPF在功能上是相关的[40-47]，SPCH和EPF在沙棘中的表达谱相似，这与它们的功能相关性一致。

综上所述，本研究发现沙棘刺和芽具有明显的形态和组织差异，二者中的差异表达基因(DEGs)主要与木质素合成以及细胞壁、表皮和气孔发育有关，参与沙棘刺顶端分生组织活动消失和刺硬化过程相关的关键基因得到初步明确。本研究结果为进一步探讨沙棘刺发育的分子机制奠定了基础。