沙棘刺和芽的組織解剖與轉(zhuǎn)錄組分析

王建武，楊 濤，嚴(yán)加坤，相微微，馮光惠，劉翠英，尚愛軍*

(1 榆林學(xué)院陜西省陜北礦區(qū)生態(tài)修復(fù)重點實驗室, 陜西榆林 719000；2 陜西省林業(yè)科學(xué)院，西安 710082；3 國家林業(yè)和草原局長柄扁桃工程技術(shù)研究中心，陜西榆林 719000)

自然界中的許多植物都有刺，但不同植物中刺的發(fā)育機制存在顯著差異。英語詞匯里有專門的單詞來區(qū)分不同的刺，即prickle(由皮層發(fā)育而來)、spine(葉子的尖銳部分)和thorn(尖銳的樹枝)[1]，玫瑰刺是植物表皮毛發(fā)育而形成的尖銳突起，通常被誤認(rèn)為是真正的葉刺或枝刺[2]。本研究在討論尖銳的植物器官時統(tǒng)一使用“刺”一詞。刺的明顯功能是防御，刺可以提供機械保護(hù)，防止草食動物啃食；刺可以幫助種子和植物擴散；刺可以增加表皮的厚度，從而減少熱量和水分的散失；一些刺可以防止陽光直射，降低光抑制，并可降低低溫產(chǎn)生的損傷[1-2]。

從生理角度來看，植物刺是有重要功能的器官，刺的發(fā)育和形成除了受自身基因控制外，環(huán)境因素對刺的發(fā)育和形成也有重要影響。動物啃食可以刺激植物產(chǎn)生更多的刺，以防御動物過度啃食[1]；土壤濕度降低會使黑枸杞的莖刺變得堅硬[3]；遮陰可以使鷹爪化的鉤刺數(shù)量增多，有利于攀緣[4]。因此從遺傳的角度來看，植物刺是一個數(shù)量性狀，而不是一或幾個基因控制的質(zhì)量性狀。刺的發(fā)育和形成由多個基因控制，其分子機理應(yīng)該是一個復(fù)雜的調(diào)控網(wǎng)絡(luò)。目前，研究人員把表皮毛作為研究擬南芥和棉花細(xì)胞增殖和生長發(fā)育的模式細(xì)胞，因此，由表皮毛發(fā)育成刺的分子機理相對清晰[2,5-9]。對于葉刺、皮刺、枝刺等多細(xì)胞構(gòu)成的刺，目前已通過遺傳和分子生物學(xué)手段克隆鑒定到一些相關(guān)基因，如玫瑰的RrTT-G1基因[10]、花椒的10個皮刺分化相關(guān)基因[11]、萵苣的LG5467628基因(混池RNA-seq測序后定位到的候選基因)[12]及生菜(做了基因定位分析，尚未具體到基因水平)[13]，但尚未建立一個調(diào)控網(wǎng)絡(luò)，而且已有的研究相互印證和支持性也不強，也進(jìn)一步說明刺的發(fā)育和形成分子機制比較復(fù)雜。不同植物的刺，甚至同種植物不同部位的刺，其發(fā)育和形成分子機理可能都有其特異性。

沙棘是多年生灌木，根具有固氮根瘤菌，耐貧瘠和干旱[14]，是優(yōu)良的防風(fēng)固沙植物，具有重要的生態(tài)價值。此外，沙棘果實具有很高的營養(yǎng)和藥用價值，可以加工沙棘汁，并從果實中提取藥用成分[15-17]。沙棘枝上有很多刺，沙棘刺含有維管組織，折斷時斷面粗糙，屬于典型的枝刺(thorn)。從外部形態(tài)來看，沙棘刺具有營養(yǎng)葉，發(fā)育早期的刺和枝芽差異很小，只是刺頂端生長點消失，而且刺的著生位置跟枝芽也沒有明顯差異，總體上刺和芽在莖上有規(guī)律地間隔開。這些特征跟黑枸杞和酸棗的枝刺有較大區(qū)別，黑枸杞和酸棗的刺不具有營養(yǎng)葉，而且長在葉腋中，因此，對沙棘刺的研究可以聚焦于比較分析去掉營養(yǎng)葉的刺和芽之間的差異，在取材上能盡量去掉背景噪音，凸顯差異。從采摘的角度來看，人們希望獲得少刺或無刺的沙棘品種，然而要通過遺傳手段調(diào)控刺的發(fā)育，培育少刺或無刺新品種，需要對沙棘刺的發(fā)育和形成分子機理有深入了解。綜上，對沙棘刺的研究既有理論意義，也有應(yīng)用前景。本研究比較了沙棘刺和芽的形態(tài)和組織差異，通過轉(zhuǎn)錄組測序分析比較了沙棘刺和芽基因表達(dá)譜，挖掘到與木質(zhì)素、細(xì)胞壁、表皮和氣孔相關(guān)的差異表達(dá)基因(differentially expressed genes，DEGs)，并通過qRT-PCR驗證了關(guān)鍵基因的表達(dá)水平。

1 材料和方法

1.1 植物材料

從7年生雄性沙棘(HippophaerhamnoidesL.ssp.sinenisRousi)上采集芽和刺，該沙棘品種多刺。種植地點位于陜西省榆林市榆陽區(qū)陜西省防沙研究院苗圃內(nèi)(海拔1 129 m；N38°19′49″，E109°42′42″)。在2021年4月22日，選擇約2 cm長的芽和刺，在開花前采集，除去營養(yǎng)葉，液氮中速凍，干冰運輸?shù)缴虾Ｃ兰具M(jìn)行轉(zhuǎn)錄組測序。樣品中所有4個生物重復(fù)都來自同一株植株，芽標(biāo)記為bud1、bud2、bud3、bud4，剌標(biāo)記為thorn1、thorn2、thorn3和thorn4。

1.2 組織切片分析

樣品在FAA中固定過夜，將固定的材料分割成0.3 cm大小，梯度乙醇脫水，石蠟包埋，Leica RM2135 石蠟切片機切片，切片厚度 8～12 μm，然后二甲苯脫蠟，梯度乙醇復(fù)水，番紅-固綠染色，中性樹膠封片，Leica DMLB顯微鏡觀察并拍照。切片用開放場或在340～380 nm激發(fā)后，用Leica DMLB熒光顯微鏡觀察，并用攝像機(Leica，DC 300F)記錄[18]。

1.3 RNA-Seq分析

使用TRIzol試劑(美國Invitrogen)從樣品中提取總RNA，并在Illumina HiSeq平臺上進(jìn)行RNA測序。RNA測序數(shù)據(jù)分析按標(biāo)準(zhǔn)化流程進(jìn)行[19-21]。每個轉(zhuǎn)錄本的表達(dá)水平被標(biāo)準(zhǔn)化為每百萬個轉(zhuǎn)錄本(transcripts per million，TPM)值[22]。在本研究中，每個轉(zhuǎn)錄本的差異表達(dá)水平以| log2(fold-change)| ≥ 1和P-value ≤ 0.05為標(biāo)準(zhǔn)，以確定芽和刺之間的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)數(shù)據(jù)庫進(jìn)行比對，對unigenes進(jìn)行注釋[20]?；蚓垲惙治霾捎肕EGA 7.0軟件。

1.4 qRT-PCR驗證

所用樣本與轉(zhuǎn)錄組分析相同。qRT-PCR在FTC-3000儀器(加拿大，F(xiàn)unglyn生物技術(shù)公司)上使用TB-Green方法(TB-Green?Fast qPCR Mix，日本Takara)進(jìn)行。用于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)。內(nèi)參基因采用通過轉(zhuǎn)錄組測序得到的沙棘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，該基因表達(dá)相對穩(wěn)定，適合作為內(nèi)參基因[24]。

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

1.5 統(tǒng)計分析

qRT-PCR試驗采用3個生物學(xué)重復(fù)，數(shù)據(jù)均為“平均值±標(biāo)準(zhǔn)偏差”。統(tǒng)計分析用SPSS 22.0。

2 結(jié)果與分析

2.1 沙棘刺和芽的形態(tài)和組織比較

成熟的沙棘刺長約3 cm，具有營養(yǎng)葉；同一年生長的刺會在秋天時從刺尖開始枯萎(圖1，A)?？偟膩碚f，沙棘芽橫切面的組織邊界明顯且層次分明；而沙棘刺橫切面的組織邊界模糊，染色區(qū)域類型單調(diào)，且刺尖中心的紅色區(qū)域較大(圖1，C～F)。這可能是因為刺尖細(xì)胞開始變得像導(dǎo)管一樣，原生質(zhì)體逐漸解體并死亡，只留下細(xì)胞壁。木質(zhì)素的自發(fā)藍(lán)色熒光顯示刺中的木質(zhì)化程度顯著高于芽，木質(zhì)化程度的增加可增加刺的硬度(圖1，G～J)。

2.2 轉(zhuǎn)錄組數(shù)據(jù)整體分析

轉(zhuǎn)錄組數(shù)據(jù)優(yōu)化組裝結(jié)果(表2)顯示，經(jīng)數(shù)據(jù)組裝共獲得81 560個unigenes(單基因)，其中60.22%的unigenes被定位，unigenes和transcript(轉(zhuǎn)錄本)的TransRate得分分別為0.303 14和0.418 81。統(tǒng)計表明，一個從頭組裝的得分為0.22(優(yōu)化得分為0.35)的TransRate分析結(jié)果將優(yōu)于50%的已發(fā)布在NCBI TSA數(shù)據(jù)庫中的分析結(jié)果[25]，故本研究的轉(zhuǎn)錄組組裝質(zhì)量相對較好。BUSCO得分C:75.0%(S:70.3%；D:4.7%)，也相對較好(表2)[26]。

表2 沙棘轉(zhuǎn)錄組數(shù)據(jù)優(yōu)化組裝結(jié)果評估表Table 2 Results of optimized assembly for sea buckthorn transcriptome

組裝完成后，通過BLAST在6個公共數(shù)據(jù)庫(NR、Swissprot、KEGG、COG、Pfam和GO)中對組裝的unigenes進(jìn)行注釋，其中分別有38 210、31 176、18 728、34 175、30 298和32 173個unigenes在上述數(shù)據(jù)中得到比對注釋。至少在一個數(shù)據(jù)庫得到注釋的unigenes數(shù)目為42 193(占所有unigenes的52%)。GO和KEGG注釋基因在所有單基因中的比例分別為39%和23%。共獲得9 385個DEGs(上調(diào)6 428個、下調(diào)2 957個)。DEGs聚類分析表明，本研究的4個生物學(xué)重復(fù)具有良好的一致性(圖2)。

為了進(jìn)一步了解DEGs的功能，將所有DEGs映射到GO數(shù)據(jù)庫，并進(jìn)行富集分析。GO分析表明DEGs被分為三大功能類別47個項目，其中BP類別、CC類別和MF類別中的項目數(shù)量分別為20、12和15。圖3顯示了前20個GO富集程度最高的項目及其DEGs數(shù)目。排名前三的GO富集項目分別是植物表皮發(fā)育(plant epidermis development)、氣孔復(fù)合體發(fā)育(complex stomatal development)和木質(zhì)素分解代謝(lignin catabolic process)過程。KEGG富集分析得到DEGs參與的138個代謝途徑，圖4顯示了前20個KEGG富集程度最高的項目及其DEGs數(shù)目，其中角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成途徑(cutin, suberine and wax biosynthesis)是排名第一的途徑。

2.3 木質(zhì)素分解代謝過程相關(guān)的DEGs

木質(zhì)素的生物合成和代謝途徑主要分為莽草酸途徑(shikimic acid pathway)、苯丙烷途徑(phenylpropanoid pathway)和木質(zhì)素特異性合成途徑(lignin-specific synthesis pathway)三個階段。迄今為止，莽草酸途徑中參與木質(zhì)素生物合成的調(diào)控機制尚未被確定。然而，在苯丙烷途徑和木質(zhì)素合成特異途徑中發(fā)現(xiàn)了能夠特異性調(diào)節(jié)木質(zhì)素合成的關(guān)鍵酶和相關(guān)轉(zhuǎn)錄因子[27-29]。本研究通過GO富集分析發(fā)現(xiàn)了20個與木質(zhì)素合成有關(guān)的DEGs，其中只有2個屬于苯丙烷途徑，即DN2241_c0_g1(肉桂醇脫氫酶1，Cinnamyl alcohol dehydrogenase 1)和DN8647_c0_g3(咖啡酸3-O-甲基轉(zhuǎn)移酶，Caffeic acid 3-O-methyltransferase)；兩種推定的蛋白(DN7364_c0_g1和DN4324_c1_g1)的具體功能尚不清楚；其他16個DEGs具有漆酶(laccase)活性，這些DEGs屬于木質(zhì)素合成特異途徑(表3)。有些漆酶能降解木質(zhì)素，有些漆酶能合成木質(zhì)素，這與漆酶的來源和合成周期有關(guān)。與芽相比，刺中只有5個具有漆酶活性的DEGs表達(dá)下調(diào)，另外11個具有漆酶活性的DEGs表達(dá)上調(diào)(表3)。其中，DN8484_c0_g2、DN8484_c0_g3、DN8484_c0_g1、DN27398_c0_g1和DN16315_c0_g1在芽中的表達(dá)非常低，幾乎在刺中特異表達(dá)(表3)。接下來，選擇6個基因，通過qRT-PCR驗證其表達(dá)水平。這些基因的相對表達(dá)水平與TPM值一致，與芽相比，DN4324_c0_g3、DN8484_c0_g2和DN16315_c0_g1在刺中的表達(dá)分別上調(diào)了9.5倍、41.1倍和13.7倍； DN7364_c0_g2、DN10443_c0_g2和DN8647_c0_g3的下調(diào)顯著，分別只有芽中表達(dá)量的5%、14%和24%(圖5，A)。

表3 沙棘芽和刺中木質(zhì)素分解代謝相關(guān)DEGsTable 3 DEGs linked with lignin catabolic process of sea buckthorn bud and thorn

2.4 細(xì)胞壁、表皮和氣孔相關(guān)的DEGs

通過KEGG途徑富集分析，角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成是最富集的途徑。因此，我們分析了屬于該途徑的DEGs。共有29個DEGs與角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成有關(guān)；這些DEGs可分為乙醇合成相關(guān)的脂肪酰輔酶A還原酶(Alcohol-forming fatty acyl-CoA reductase)、ω-羥基棕櫚酸酯O-阿魏酰轉(zhuǎn)移酶(Omega-hydroxypalmitate O-feruloyl transferase)、細(xì)胞色素P450(Cytochrome P450)、熱頭蛋白(HOTHEAD protein)和其他類別(表4)。許多基因在芽和刺中的表達(dá)水平具有相反的趨勢，要么在刺中表達(dá)，而在芽中表達(dá)非常低，要么在芽中表達(dá)，而在刺中表達(dá)非常低(表4)。例如，DN53228_c0_g1、DN15725_c0_g4、DN41003_c0_g1、DN14231_c0_g1和DN16701_c0_g1在芽中的表達(dá)水平很低，但在刺中高表達(dá)，而DN424_c0_g1、DN48710_c0_g1和DN27861_c0_g1在刺中的表達(dá)水平很低，但在芽中高表達(dá)；特別是，刺中DN27861_c0_g1的TPM值為0(表4)。接下來，從上述基因中選擇4個基因，通過qRT-PCR驗證其表達(dá)。這些基因的相對表達(dá)水平與TPM值一致，與芽相比，刺中的DN8861_c0_g2和DN14231_c0_g1分別上調(diào)3.6倍和22.1倍；DN424_C0_g1和DN27861_C0_g1的下調(diào)顯著，分別只有芽中表達(dá)量的3%倍和6%倍(圖5，B)。

表4 沙棘芽和刺中角質(zhì)、木栓質(zhì)和蠟質(zhì)合成相關(guān)DEGsTable 4 DEGs linked with cutin, suberine and wax biosynthesis of sea buckthorn bud and thorn

植物表皮發(fā)育和氣孔復(fù)合體發(fā)育是GO富集分析結(jié)果中的前2個項目，因此，本研究進(jìn)一步分析了相關(guān)DEGs。植物表皮發(fā)育與氣孔復(fù)合體發(fā)育密切相關(guān)。本研究鑒定的參與植物表皮發(fā)育和氣孔復(fù)合體發(fā)育的DEGs是相同的；這些DEGs可分為表皮形成因子(Epidermal patterning factor，EPF)、轉(zhuǎn)錄因子(Transcription factor)、多聚半乳糖醛酸酶(Polygalacturonase)和其他類別(表5)。與芽相比，所有DEGs在刺中均下調(diào)表達(dá)，DN55365_c0_g1、DN11032_c0_g1、DN13014_c0_g1、DN29670_c0_g1、DN46942_c0_g1、DN16737_c0_g1、DN13978_c0_g1和DN16781_c0_g1在刺中的表達(dá)水平非常低，但在芽中高表達(dá)；特別是，刺中DN29670_c0_g1的TPM值為0(表5)。接下來，從上述基因中選擇4個基因，并通過qRT-PCR驗證其表達(dá)。這些基因的相對表達(dá)水平與TPM值一致，與芽相比，刺中的DN11032_c0_g1、DN13014_c0_g1、DN16737_c0_g1和DN13978_c0_g1顯著下調(diào)，分別只有芽中表達(dá)量的6%、13%、3%和4%(圖5，C)。

表5 沙棘芽和刺中表皮/氣孔發(fā)育相關(guān)DEGsTable 5 DEGs linked with epidermis/stomatal development of sea buckthorn bud and thorn

3 討 論

前人研究認(rèn)為，無論側(cè)芽最終是發(fā)育成刺還是枝，都是在芽個體發(fā)育的早期階段決定的，刺芽通常比枝芽產(chǎn)生更小、更少的三角形葉原基[30]。刺發(fā)育過程中的顯著變化是頂端分生組織縮小，最終變硬，同時伴隨分生組織活動和葉片起始發(fā)育的停止[30]。

頂端分生組織的分子調(diào)控機制非常復(fù)雜，僅基于轉(zhuǎn)錄組數(shù)據(jù)很難進(jìn)行深入、明確的分析。最近，人們發(fā)現(xiàn)水稻中也存在與擬南芥HTH基因高度同源的HOTHEAD(HTH)基因[31-32]。ONI3在水稻莖尖細(xì)胞中特異表達(dá)，ONI3突變體在莖尖細(xì)胞中表現(xiàn)出發(fā)育缺陷，導(dǎo)致莖尖發(fā)育異常[33]。類似地，Mini1在水稻幼嫩莖尖表達(dá)，但在根中不表達(dá)，Mini1突變體芽的發(fā)育提前終止[34]。因此，推測水稻HTH基因和擬南芥HTH基因具有不同的功能；水稻HTH基因在維持莖尖分生組織(SAM)活性和促進(jìn)水稻莖尖發(fā)育方面起著重要作用[33-34]。在本研究中，挖掘得到5個HTH基因，其中unigene DN27861_c0_g1在刺中不表達(dá)，DN48710_c0_g1在刺中的表達(dá)也非常低，DN34353_c0_g1在刺中表達(dá)增加，而其他2個基因在刺中表達(dá)降低。推測其中一些基因可能參與角質(zhì)合成，類似于擬南芥HTH基因的功能，一些可能參與分生組織活動。

細(xì)胞壁木質(zhì)化可顯著增加刺的硬度。因此，沙棘刺的木質(zhì)化程度顯著高于芽，這與沙棘刺的硬化一致。木質(zhì)素分解代謝過程是GO富集分析結(jié)果中的第3項，表明刺和芽中的木質(zhì)素代謝存在顯著差異。大多數(shù)與木質(zhì)素分解代謝過程有關(guān)的DEGs屬于漆酶；漆酶負(fù)責(zé)將氨酚醇聚合成木質(zhì)素聚合物[35-36]。單一木質(zhì)素的組成在植物種類、組織和細(xì)胞類型之間變化很大，甚至在單個細(xì)胞內(nèi)的獨立細(xì)胞壁層之間差異也很大，該組成主要影響最終木質(zhì)素聚合物的化學(xué)結(jié)構(gòu)和物理化學(xué)性質(zhì)[35]。漆酶的組織定位以及它們的生化特性決定了木質(zhì)素聚合的組織模式[35]。沙棘漆酶生化特性的深入研究，可將其作為經(jīng)基因工程手段調(diào)控沙棘刺硬度的候選基因。

在刺的發(fā)育過程中，細(xì)胞的功能類型逐漸趨于單一化，也就是說，它們變硬并起到機械支撐的作用。因此，與芽相比，刺細(xì)胞的細(xì)胞壁成分發(fā)生了顯著變化。本研究結(jié)果支持這一觀點。首先，組織切片分析表明，刺細(xì)胞的染色結(jié)果相對單一，層次性不如芽切片明顯，這表明刺細(xì)胞的功能在形態(tài)學(xué)上得到了簡化。其次，KEGG富集分析表明，角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成是富集程度最高的途徑，乙醇形成相關(guān)的脂肪酰輔酶A還原酶、ω-羥基棕櫚酸酯O-阿魏酰轉(zhuǎn)移酶、細(xì)胞色素P450、HOTHEAD蛋白家族在角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成途徑中有多種異構(gòu)體，這些蛋白質(zhì)異構(gòu)體屬于進(jìn)化中的旁系同源物，應(yīng)具有明顯的功能分化。因此，進(jìn)一步研究這些同源蛋白的生物化學(xué)功能，分析它們參與的代謝途徑，對闡明沙棘刺和芽之間角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成的差異具有重要意義。第三，沙棘刺中表皮和氣孔的發(fā)育可能逐漸停止，本研究發(fā)現(xiàn)參與表皮和氣孔發(fā)育的所有DEGs在刺中均下調(diào)表達(dá)。bHLH轉(zhuǎn)錄因子SPEECHLESS(SPCH)在氣孔發(fā)育中扮演關(guān)鍵角色[37-39]，SPCH整合了內(nèi)源信號，并作為激素和環(huán)境信號的受體，在發(fā)育過程中調(diào)節(jié)氣孔密度和發(fā)育模式[37-39]。SPCH的表達(dá)水平受分泌型的表皮形成因子(EPF)調(diào)節(jié)，SPCH與EPF在功能上是相關(guān)的[40-47]，SPCH和EPF在沙棘中的表達(dá)譜相似，這與它們的功能相關(guān)性一致。

綜上所述，本研究發(fā)現(xiàn)沙棘刺和芽具有明顯的形態(tài)和組織差異，二者中的差異表達(dá)基因(DEGs)主要與木質(zhì)素合成以及細(xì)胞壁、表皮和氣孔發(fā)育有關(guān)，參與沙棘刺頂端分生組織活動消失和刺硬化過程相關(guān)的關(guān)鍵基因得到初步明確。本研究結(jié)果為進(jìn)一步探討沙棘刺發(fā)育的分子機制奠定了基礎(chǔ)。