青花菜花蕾轉(zhuǎn)錄組測(cè)序分析及蠟粉合成相關(guān)基因挖掘

張振超　陶美奇　潘永飛　戴忠良　姚悅梅

摘要：【目的】對(duì)青花菜花蕾進(jìn)行轉(zhuǎn)錄組測(cè)序分析，并挖掘與蠟粉合成相關(guān)基因，為探明青花菜花球表面蠟粉形成的分子機(jī)制提供理論參考?！痉椒ā糠謩e提取野生型和蠟粉缺失型青花菜花球總RNA，采用Illumina HiSeq^TM2500平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序，獲得高質(zhì)量Clean reads，采用Trinity進(jìn)行序列組裝后獲得青花菜Unigene庫(kù)，將獲得的Unigene序列與Nr、Nt、KEGG、Pfam、KOG/COG、Swiss-Prot和GO數(shù)據(jù)庫(kù)比對(duì)，獲得基因功能注釋信息;使用DESeq2進(jìn)行差異表達(dá)分析。【結(jié)果】共獲得44.68 Gb Clean data，De novo組裝得到41244條Unigenes，N50長(zhǎng)度為1847 bp。從所獲得的Unigenes中篩選出8685個(gè)差異表達(dá)基因（DEGs）（上調(diào)基因5747個(gè)，下調(diào)基因2938個(gè)），共有8038個(gè)基因被注釋到不同數(shù)據(jù)庫(kù)，其中，5220個(gè)基因注釋到Pfam數(shù)據(jù)庫(kù);2066個(gè)基因注釋到COG數(shù)據(jù)庫(kù)，3866個(gè)基因注釋到KOG數(shù)據(jù)庫(kù);2580個(gè)差異表達(dá)基因被注釋到75個(gè)轉(zhuǎn)錄因子家族中，注釋最多的是MYB家族（235個(gè)）;GO數(shù)據(jù)庫(kù)中6095個(gè)差異表達(dá)基因注釋到細(xì)胞組分、分子功能和生物學(xué)過(guò)程三大類(lèi)的52個(gè)功能分類(lèi);KEGG數(shù)據(jù)庫(kù)中，1671個(gè)差異表達(dá)基因富集到138條代謝通路，其中13個(gè)差異表達(dá)基因與脂肪酸合成有關(guān)，7個(gè)差異表達(dá)基因與蠟粉生物合成途徑有關(guān)。【結(jié)論】轉(zhuǎn)錄因子MYB家族在調(diào)控青花菜蠟粉合成中發(fā)揮重要作用。蠟粉合成過(guò)程中相關(guān)酶基因的差異表達(dá)是調(diào)控青花菜蠟粉合成的關(guān)鍵，尤其是野生型和蠟粉缺失突變體中特異性表達(dá)的差異表達(dá)基因，可作為后續(xù)研究青花菜花球表面蠟粉形成分子機(jī)制的對(duì)象。

關(guān)鍵詞： 青花菜;花蕾;蠟粉缺失突變體;轉(zhuǎn)錄組;基因功能注釋

中圖分類(lèi)號(hào)： S635.3? ? ? ? ? ? ? ? ? ? ? ? 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)：2095-1191（2022）03-0607-11

Transcriptome analysis and mining of genes related to wax powder synthesis of broccoli flower buds

ZHANG Zhen-chao， TAO Mei-qi， PAN Yong-fei， DAI Zhong-liang， YAO Yue-mei

（Zhenjiang Institute of Agricultural Sciences in Hilly Area of Jiangsu Province， Jurong， Jiangsu? 212400， China）

Abstract：【Objective】To perform transcriptome sequencing on broccoli buds， and mine their genes related to wax powder synthesis， so as to lay a foundation for the discovery of the molecular mechanism of wax powder formation on the surface of broccoli spherules. 【Method】Total RNA from wild-type and wax-powder-deficient broccoli flower bulbs was extracted， and Illumina HiSeq^TM2500 platform was used for transcriptomes sequencing to obtain high-quality Clean reads， and Trinity software was used for sequence assembly to obtain the broccoli Unigene library. The obtained Unigenes sequence was compared with that in Nr， Nt， KEGG， Pfam， KOG/COG， Swiss-Prot， GO databases to obtain gene function annotation information; DESeq2 software was used for differential expression analysis. 【Result】A total of 44.68 Gb Clean data was generated. And 41244 unigenes with N50 length of 1847 bp， were obtained by De novo assembly method. 8685 differentially expressed genes （DEGs） （5747 up-regulated genes and 2938 down-regulated genes） were obtained in the unigenes， of which 8038 had annotated information. Of these， 5220 genes were annotated to the Pfam database; 2066 genes annotated to COG database and 3866 genes annotated to KOG database annotation analysis showed that 36230 unigenes had homologens in different public protein databases. A total of 2580 DEGs were annotated to 75 transcription factor families， and the most annotated was MYB family （235）. In GO database， 6095 DEGs were annotated to 52 functional groups in three categories：cell components， molecular functions and biological processes. In KEGG database， 1671 DEGs were classified into 138 metabolic pathway branches， and 13 DEGs were annotated to fatty acid synthesis related pathway， and 7 DEGs were annotated to wax powder biosynthesis related pathway. 【Conclusion】The MYB family of transcription factors plays an important role in the regulation of wax powder synthesis in broccoli. The differential expression of related enzyme genes during wax powder synthesis is the key to regulating wax powder synthesis in broccoli. In particular， the DEGs in wild-type and wax powder-deficient mutants can be used for further study on the molecular mechanism of wax powder formation on the surface of broccoli.BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

Key words： broccoli; flower bud; wax powder deletion mutant; transcriptome; gene function annotation

Foundation items： Jiangsu Key Research and Development Project （BE2021376）; Jurong Science and Technology Innovation Fund Project （ZB22105）;Youth Fund Project of Zhenjiang Academy of Agricultural Sciences（QNJJ2021003）

0 引言

【研究意義】植物蠟粉，又稱(chēng)蠟質(zhì)，是覆蓋于植物組織外在表皮上脂質(zhì)成分的統(tǒng)稱(chēng)，由脂肪酸、烷烴、醇類(lèi)、醛類(lèi)和酮類(lèi)等親脂性化合物組成，一般呈綠灰色或灰白色霜狀（Domínguez et al.，2011;張曦等，2013）。蠟粉不僅是植物表皮細(xì)胞外的疏水屏障，還是重要的防衛(wèi)機(jī)構(gòu)，在植物生長(zhǎng)發(fā)育過(guò)程中發(fā)揮反射紫外線(xiàn)、抵御病菌入侵、防止食草昆蟲(chóng)侵害等作用（Oliveira et al.，2004;王立山等，2018）。青花菜（Brassica oleracea L. var. italica Plenck）口感清脆、營(yíng)養(yǎng)價(jià)值高且富含抗癌活性成分硫代葡萄糖苷，被譽(yù)為“蔬菜皇冠”。青花菜花球表面因覆蓋一層蠟粉而呈灰綠色或灰白色，其中以蠟粉含量多、顏色灰綠的花球較受市場(chǎng)青睞（張振超等，2021）。因此，開(kāi)展青花菜花球表面蠟粉合成代謝的分子機(jī)制研究及相關(guān)功能基因挖掘，對(duì)提高青花菜優(yōu)異種質(zhì)資源創(chuàng)制和新品種選育具有重要意義?！厩叭搜芯窟M(jìn)展】近年來(lái)，擬南芥、油菜、大白菜、玉米、大麥等作物的蠟粉缺失突變體均有報(bào)道。擬南芥和大麥的蠟粉突變位點(diǎn)被稱(chēng)為eceriferum（cer），油菜和玉米中的突變位點(diǎn)被稱(chēng)為glossy（張曦等，2013;Saet and Mi，2015）。目前從大麥蠟粉突變體中獲得了85個(gè)cer突變位點(diǎn)（Wettstein and Chua，1987）;從擬南芥中篩選出32個(gè)突變位點(diǎn)，其中與蠟粉合成相關(guān)的基因有20多個(gè)，如CER1、CER2、KCS1和WAX2等，這些基因參與蠟粉直接合成、轉(zhuǎn)錄和激素水平調(diào)控等（Pruitt et al.，2000;Jenks et al.，2002;Chen et al.，2003）。對(duì)油菜（Pu et al.，2013）、白菜（Zhang et al.，2013）和結(jié)球甘藍(lán)（Liu et al.，2017，2018）中的蠟粉缺失基因進(jìn)行區(qū)間精細(xì)定位，并將區(qū)間內(nèi)與擬南芥CER1和CER4基因同源的基因確定為候選目的基因。Pyee和Kolattukudy（1995）以青花菜中編碼WAX9蛋白的cDNA序列作為探針從基因文庫(kù)中分離出青花菜表面蠟粉4個(gè)同源基因wax9A、wax9B、wax9C和wax9D，并對(duì)蠟粉遺傳規(guī)律進(jìn)行了系統(tǒng)研究，結(jié)果顯示，所有基因均在葉片和花蕾中表達(dá)，但均不在根中表達(dá)， wax9A、wax9B和wax9D在莖和花中也有表達(dá)。近年來(lái)，隨著高通量測(cè)序技術(shù)快速發(fā)展，該技術(shù)被廣泛應(yīng)用于基因表達(dá)與植物表型關(guān)系的研究。利用該技術(shù)對(duì)同一物種不同樣本進(jìn)行轉(zhuǎn)錄組測(cè)序分析，在獲得大量功能基因的基礎(chǔ)上比較樣本間存在表達(dá)差異的基因，并根據(jù)生物信息學(xué)分析探究差異表達(dá)基因功能分類(lèi)，篩選出與差異表型相關(guān)的基因群，進(jìn)而發(fā)掘新的功能基因（林琿等，2018;吳寧柔等，2018;張振超等，2018）?！颈狙芯壳腥朦c(diǎn)】目前高通量測(cè)序技術(shù)已成功應(yīng)用到植物代謝機(jī)制和基因調(diào)控功能的研究中，如青花菜胚胎發(fā)育初期小孢子經(jīng)高通量測(cè)序篩選得到大量與胚胎發(fā)育相關(guān)的代謝通路和候選基因（張振超等，2018）。但未見(jiàn)利用轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)蠟粉合成基因挖掘及功能分析的研究報(bào)道。【擬解決的關(guān)鍵問(wèn)題】以野生型青花菜花球與蠟粉缺失突變體為試材，采用高通量測(cè)序技術(shù)進(jìn)行轉(zhuǎn)錄組測(cè)序分析，對(duì)差異表達(dá)基因進(jìn)行功能注釋?zhuān)?duì)其中的脂肪酸和蠟粉合成相關(guān)基因進(jìn)行KEGG代謝通路富集分析，以篩選出合成途徑中的關(guān)鍵酶基因，為深入探究青花菜花球蠟粉生物合成代謝的分子機(jī)制、挖掘相關(guān)基因及遺傳育種應(yīng)用打下基礎(chǔ)。

1 材料與方法

1. 1 試驗(yàn)材料

青花菜野生型（WT20）與蠟粉缺失突變體（MT20）（圖1）由江蘇丘陵地區(qū)鎮(zhèn)江農(nóng)業(yè)科學(xué)研究所提供。植物總RNA提取試劑盒購(gòu)自寶生物工程（大連）有限公司。主要儀器設(shè)備：NanoDrop ND-2000微量紫外分光光度計(jì)（NanoDrop technology，美國(guó)）、Qubit2.0核酸蛋白熒光定量?jī)x（美國(guó)賽默飛世爾科技有限公司）和Aglient 2100芯片生物分析儀（美國(guó)安捷科技有限公司）。

1. 2 試驗(yàn)方法

1. 2. 1 總RNA提取 參照植物總RNA提取試劑盒說(shuō)明提取青花菜總RNA。設(shè)置3個(gè)生物學(xué)重復(fù)，每個(gè)重復(fù)3株。采用1%瓊脂糖電泳檢測(cè)RNA完整性，分別用NanoDrop ND-2000微量紫外分光光度計(jì)、Qubit2.0核酸蛋白熒光定量?jī)x和Aglient 2100芯片生物分析儀定量檢測(cè)RNA樣品的濃度、純度和完整性等，以保證使用合格的樣品進(jìn)行轉(zhuǎn)錄組測(cè)序。

1. 2. 2 轉(zhuǎn)錄組測(cè)序及生物信息學(xué)分析 轉(zhuǎn)錄組測(cè)序參照張振超等（2018）的方法。利用Illumina HiSeqTM 2500測(cè)序平臺(tái)進(jìn)行高通量測(cè)序，形成雙端測(cè)序數(shù)據(jù)，測(cè)序讀長(zhǎng)為PE150。將獲得的原始數(shù)據(jù)進(jìn)行過(guò)濾，去除其中的接頭序列及低質(zhì)量reads，獲得高質(zhì)量Clean reads，采用Trinity進(jìn)行序列組裝后獲得青花菜Unigenes庫(kù)。利用Blast Version 2.2.26和HMMER 3.1進(jìn)行測(cè)序飽和度、基因覆蓋度和冗余序列分析，將獲得的Unigenes序列在Nr（NCBI non-redundant protein sequences）、Nt（NCBI non-redundant nucleotide sequence）、KEGG（Kyoto encyclopedia of genes and genomes）、Pfam（Protein family）、KOG/COG（Clusters of orthologous groups of proteins）、Swiss-Prot（A manually annotated and revie-wed protein sequence database）、KEGG、GO（Gene ontology）數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，獲得基因功能注釋信息。BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

1. 2. 3 差異表達(dá)基因篩選及分析 采用FPKM值（Fragments Per Kilobase of transcript per Million mapped reads）表示基因表達(dá)水平，使用DESeq2進(jìn)行差異表達(dá)分析。將|log2FC（Fold change，差異倍數(shù)）|≥1且FDR（False discovery rate，錯(cuò)誤發(fā)現(xiàn)率）<0.05作為差異表達(dá)基因（Differentially expressed genes，DEGs）篩選標(biāo)準(zhǔn)。采用Benjamini-Hochberg法對(duì)差異顯著性P值進(jìn)行校正，并將FDR作為差異表達(dá)基因篩選的關(guān)鍵指標(biāo)。參照GO數(shù)據(jù)庫(kù)，將差異表達(dá)基因按照參與的生物學(xué)過(guò)程、細(xì)胞組成和分子功能進(jìn)行分類(lèi);參照KEGG數(shù)據(jù)庫(kù)分析差異表達(dá)基因參與的代謝途經(jīng)，并篩選出與蠟粉合成代謝相關(guān)的基因;對(duì)差異表達(dá)基因的轉(zhuǎn)錄因子進(jìn)行預(yù)測(cè)分析。

2 結(jié)果與分析

2. 1 轉(zhuǎn)錄組測(cè)序數(shù)據(jù)組裝結(jié)果

本研究獲得的原始數(shù)據(jù)經(jīng)過(guò)測(cè)序質(zhì)量控制和數(shù)據(jù)分析后，共得到44.68 Gb Clean data，GC含量在45%～47%，Q30堿基百分比超過(guò)92.92%，表明測(cè)序質(zhì)量較好。經(jīng)De novo拼接組裝，共得到164914條Transcripts和41244條Unigenes，N50長(zhǎng)度分別為1954和1847 bp（表1），表明數(shù)據(jù)組裝完整性較高，可進(jìn)行后續(xù)生物信息學(xué)分析。對(duì)組裝的基因片段長(zhǎng)度進(jìn)行分析，結(jié)果發(fā)現(xiàn)Transcripts和Unigenes中1000～2000 bp所占比例最高，分別為37.68%和32.33%。

2. 2 Unigene功能注釋結(jié)果

使用BLAST和HMMER將獲得的Unigenes序列分別與Nr、Swiss-Prot、GO、COG、KOG、KEGG、Pfam數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果發(fā)現(xiàn)有36230條Unigenes（≥300 bp）獲得注釋信息，占比87.8%。在獲得注釋的Unigenes中，在Nr數(shù)據(jù)庫(kù)注釋的Unigenes最多，為35971條，其次是GO數(shù)據(jù)庫(kù)，為27966條，再次是SWISS-PROT數(shù)據(jù)庫(kù)，為26540條，而在KEGG數(shù)據(jù)庫(kù)注釋的Unigenes最少，有8566條（表2）。在獲得注釋的Unigenes中，長(zhǎng)度為300～1000 bp的Unigenes數(shù)為14999條，長(zhǎng)度≥1000 bp的Unigenes有21231條（表2）。

2. 3 差異表達(dá)基因篩選及功能分析結(jié)果

2. 3. 1 差異表達(dá)基因的篩選 采用DESeq2進(jìn)行樣品組間差異表達(dá)分析，以|log2FC|≥1且FDR<0.05作為篩選標(biāo)準(zhǔn)獲得差異表達(dá)基因，并進(jìn)行層次聚類(lèi)分析。將表達(dá)行為相同或相似的基因進(jìn)行歸類(lèi)，獲得差異基因聚類(lèi)圖譜（圖2和圖3）。在WT20和MT20樣本的41244條Unigenes中，差異表達(dá)基因共8685個(gè)，其中上調(diào)基因5747個(gè)，下調(diào)基因2938個(gè)，這些基因是后續(xù)研究的重點(diǎn)。

2. 3. 2 差異表達(dá)基因的功能注釋 將篩選到的8685個(gè)差異表達(dá)基因與數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，共有8038個(gè)基因被注釋?zhuān)渲?095個(gè)基因注釋到GO數(shù)據(jù)庫(kù)，1671個(gè)基因注釋到KEGG數(shù)據(jù)庫(kù)，2066個(gè)基因注釋到COG數(shù)據(jù)庫(kù)，3866個(gè)基因注釋到KOG數(shù)據(jù)庫(kù)，5220個(gè)基因注釋到Pfam數(shù)據(jù)庫(kù)，6059個(gè)基因注釋到Swiss-Prot數(shù)據(jù)庫(kù)，8020個(gè)注釋到Nr數(shù)據(jù)庫(kù)（表3）。

2. 3. 3 差異表達(dá)基因轉(zhuǎn)錄因子注釋 在8038個(gè)被注釋的差異表達(dá)基因中，有2580個(gè)被注釋為轉(zhuǎn)錄因子（表3），分布在75個(gè)轉(zhuǎn)錄因子家族中（圖4），其中，注釋為MYB家族的基因最多，為235個(gè);富集基因數(shù)為100～200個(gè)的轉(zhuǎn)錄因子家族有7個(gè)：AP2/ERF（162個(gè)）、B3（107個(gè)）、C2C2（131個(gè)）、C2H2（118個(gè)）、NAC（137個(gè)）、WRKY（105個(gè)）和bHLH（119個(gè)）;富集基因數(shù)為50～99個(gè)的轉(zhuǎn)錄因子家族有6個(gè)，富集基因數(shù)為0～49個(gè)的轉(zhuǎn)錄因子家族有61個(gè);注釋到其他轉(zhuǎn)錄因子家族的基因有120個(gè)。

2. 3. 4 差異表達(dá)基因的GO功能注釋分析結(jié)果

6095個(gè)差異表達(dá)基因注釋到細(xì)胞組分（Cell component）、分子功能（Molecular function）和生物學(xué)過(guò)程（Biological process）三大類(lèi)別的52個(gè)功能分類(lèi)。其中，共有4710個(gè)基因富集到細(xì)胞組分中，其中富集程度最高的是細(xì)胞（Cell），為4576個(gè)基因，其次為細(xì)胞部件（Cell part），為4562個(gè)基因;共有4497個(gè)基因富集到分子功能，富集程度最高的是結(jié)合（Binding），為3183個(gè)基因;共有4771個(gè)基因被富集到生物學(xué)過(guò)程，其中富集程度最高的是細(xì)胞過(guò)程（Cellular process），為3252個(gè)基因，其次是代謝過(guò)程（Metabolic process），為3237個(gè)基因。

2. 3. 5 差異表達(dá)基因的KEGG代謝通路富集分析結(jié)果 KEGG是系統(tǒng)分析基因功能和基因組信息的數(shù)據(jù)庫(kù)。對(duì)差異表達(dá)基因的代謝通路進(jìn)行分析有助于進(jìn)一步解讀基因的功能。將獲得的8685個(gè)顯著差異表達(dá)基因與KEGG數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果（圖6）發(fā)現(xiàn)，有1671個(gè)顯著差異表達(dá)基因富集到代謝（Meta-bolism）、遺傳信息處理（Genetic information proces-sing）、環(huán)境信息處理（Environmental information processing）、細(xì)胞過(guò)程（Cellular processes）、人類(lèi)疾?。℉uman diseases）和有機(jī)系統(tǒng)（Organismal systems）六大類(lèi)的21個(gè)功能分類(lèi)，共138條代謝通路。其中顯著富集通路有9條（FDR≤0.05），由于P越小，富集程度越高，因此最顯著富集通路是硫甙生物合成（Glucosinolate biosynthesis），P為2.3948e-07;富集基因最多的是次生代謝物的生物合成（Biosynthesis of secondary metabolites），為245個(gè)基因，占19.25%，其次是植物信號(hào)轉(zhuǎn)導(dǎo)（Plant hormone signal transduction），為110個(gè)基因，占8.64%（表4）。BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

2. 3. 6 脂肪酸和蠟粉合成相關(guān)基因表達(dá)分析結(jié)果

擬南芥蠟質(zhì)的生物合成一般發(fā)生于表皮細(xì)胞內(nèi)，具體過(guò)程為含有16～18個(gè)C的?；d體蛋白（C16-C18 acyl-acyl carrier proteins，ACP）前體在質(zhì)體中被合成后便被酰基載體蛋白硫酯酶（Fatty acyl-ACP thioesterases， FATA and FATB）水解生成含有16～18個(gè)C的脂肪酸，后經(jīng)長(zhǎng)鏈酰基輔酶A合成酶（Long chain acylcoenzyme A synthetases， LACS）催化合成胞質(zhì)?；o酶A后運(yùn)輸?shù)絻?nèi)質(zhì)網(wǎng)中（Bonaventure et al.，2003）。根據(jù)上述蠟粉生物合成過(guò)程，本研究重點(diǎn)分析脂肪酸代謝（Fatty acid metabolism）、脂肪酸降解（Fatty acid degradation）、脂肪酸生物合成（Fatty acid biosynthesis）、不飽和脂肪酸生物合成（Biosynthesis of unsaturated fatty acids）、脂肪酸延伸（Fatty acid elongation）及角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成（Cutin，suberine and wax biosynthesis）通路。根據(jù)KEGG數(shù)據(jù)庫(kù)的功能釋義和已知研究報(bào)道對(duì)青花菜蠟粉合成相關(guān)差異表達(dá)基因進(jìn)行分析，結(jié)果如表5所示。有6個(gè)基因（TRINITY_DN20763_c0_g5、TRINITY_ DN18264_c0_g2、TRINITY_DN17853_c0_g1、TRINITY_DN20008_c0_g1、TRINITY_DN20271_c2_g1和TRINITY_DN 22032_c1_g1）注釋到K01897，對(duì)應(yīng)酶名稱(chēng)為ACSL（?；o酶A合成酶長(zhǎng)鏈家族成員），其在數(shù)據(jù)庫(kù)中分別注釋為L(zhǎng)ACS1（長(zhǎng)鏈?；o酶A合成酶）、LACS3、LACS7、LACS8和AAE16（?；罨?6）。TRINITY_DN20589_c1_g7注釋到K10258，對(duì)應(yīng)酶為CER10（ECERIFERUM 10）。TRINITY_ DN19765_c1_g2注釋到K10781，其在數(shù)據(jù)庫(kù)中注釋為FATB（脂肪酰基ACP硫酯酶B）。TRINITY_ DN19127_c1_g5注釋到K10782，對(duì)應(yīng)的酶為FATA（脂肪酰基ACP硫酯酶A），其在數(shù)據(jù)庫(kù)中注釋為AtFaTA（擬南芥硫脂酶）。4個(gè)差異表達(dá)基因（TRINITY_DN22206_c1_g4、TRINITY_DN23613_c0_g1、TRINITY_DN18550_c2_g3和TRINITY_DN18550_ c2_g6）注釋為KCS（酮酰輔酶A合酶），其在數(shù)據(jù)庫(kù)中注釋分別為KCS7、KCS9、KCS16和KCS21。

7個(gè)差異表達(dá)基因注釋到角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成通路，其中TRINITY_DN16904_c0_g1注釋到K15404，注釋為無(wú)蠟粉基因CER1（ECERIFERUM 1）;3個(gè)基因（TRINITY_DN21697_c0_g1、TRINITY_ DN21913_c1_g10和TRINITY_DN10321_c0_g1）注釋到K13356，均注釋為FAR（脂肪酰輔酶A還原酶），但在數(shù)據(jù)庫(kù)中分別注釋為FAR2、CER4和FAR7;TRINITY_DN21660_c0_g4注釋到K15398，注釋為CYP86A4S（細(xì)胞色素P450 86A8），在數(shù)據(jù)庫(kù)中注釋為CYP86A8;TRINITY_DN21090_c3_g4和TRINITY_DN21090_c3_g3注釋為CYP86B1（細(xì)胞色素P450 86B1）（表5）。

根據(jù)WT20和MT20的FPKM值（表5）可知，有5個(gè)差異表達(dá)基因TRINITY_DN17853_c0_g1、TRINITY_ DN20008_c0_g1、TRINITY_DN19765_c1_g2、TRINITY_DN19127_c1_g5和TRINITY_DN16904_c0_g1表達(dá)下調(diào)，其余均上調(diào)表達(dá)。在WT20中，有3個(gè)差異表達(dá)基因TRINITY_DN18550_c2_g6、TRINITY_ DN10321_c0_g1和TRINITY_DN21090_c3_g4的FPKM值為0，在MT20中，TRINITY_ DN16904_c0_g1的FPKM值為0。TRINITY_DN21697_c0_g1的log2 （MT20/WT20）絕對(duì)值最高為10.68，其次是TRINITY_DN16904_c0_g1為10.45，最小的是TRINITY_ DN17853_c0_g1為1.02。

3 討論

本研究將野生型青花菜和蠟粉缺失實(shí)變體的轉(zhuǎn)錄組測(cè)序結(jié)果與Nr、Nt、KEGG、Pfam、KOG/COG、Swiss-Prot和GO數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果發(fā)現(xiàn)有36230條（占比87.8%）Unigenes得到注釋?zhuān)传@得注釋的基因有5014條，可能是由于測(cè)序片段過(guò)短、獲得的特殊基因在數(shù)據(jù)庫(kù)中的注釋信息缺乏等原因，與在青花菜（張振超等，2018）、花椰菜（林琿等，2018）、油菜（吳寧柔等，2018）等植物轉(zhuǎn)錄組結(jié)果相似。對(duì)青花菜樣本的表達(dá)譜進(jìn)行比對(duì)，獲得差異表達(dá)基因8685個(gè)，共有2580個(gè)差異表達(dá)基因被注釋到75個(gè)轉(zhuǎn)錄因子家族中，其中注釋為MYB家族的基因最多，為235個(gè)，占9.10%，表明MYB家族轉(zhuǎn)錄因子參與了調(diào)控蠟粉合成。Dubos等（2010）的研究結(jié)果也證實(shí)轉(zhuǎn)錄因子如MYB可通過(guò)獨(dú)立或協(xié)同的方式調(diào)控某些基因的表達(dá)從而實(shí)現(xiàn)蠟粉的生物合成。

植物表皮蠟粉主要是由特長(zhǎng)鏈飽和脂肪酸的衍生物組成，還包括萜類(lèi)化合物和次生代謝物等（Bourdenx et al.，2011），其合成過(guò)程需要多種酶的共同參與，如?；d體蛋白硫脂酶（FATA和FATB）、長(zhǎng)鏈?；o酶A合成酶（LACS）、β-酮脂酰CoA合成酶（KCS）、β-酮脂酰CoA還原酶（KCR）、β-羥?；o酶A水解酶（HCD）和烯?；?CoA還原酶（ECR）、脂肪酸延伸酶（FAE）等（Kunst and Samuels，2003;Bourdenx et al.，2011）。本研究根據(jù)青花菜野生型（WT20）和蠟粉缺失突變體（MT20）的轉(zhuǎn)錄組信息與KEGG數(shù)據(jù)庫(kù)比對(duì)結(jié)果，重點(diǎn)分析脂肪酸代謝、脂肪酸降解、脂肪酸生物合成、不飽和脂肪酸生物合成、脂肪酸延伸及角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成6個(gè)通路，并參照GO、Nr、Nt、Pfam、KOG/COG、Swiss-Prot、KEGG的注釋信息進(jìn)行分析，結(jié)果發(fā)現(xiàn)有6個(gè)差異表達(dá)基因注釋為ACSL，其中4條基因在數(shù)據(jù)庫(kù)中注釋為L(zhǎng)ACS1、LACS3、LACS7和LACS8，2條基因注釋為AAE16。LACS1、LACS3、LACS7和LACS8是長(zhǎng)鏈脂肪酸合成關(guān)鍵酶，但行使的功能各不相同。LACS1編碼長(zhǎng)鏈酰基輔酶 A 合成酶，可催化合成具有活性的脂酰輔酶A，并參與蠟粉代謝（Lü et al.，2009）;LACS3參與脂質(zhì)的跨膜運(yùn)輸和細(xì)胞內(nèi)運(yùn)輸（Pulsifer et al.，2012）;LACS7編碼氧化物酶體LACS蛋白，參與脂肪酸的β-氧化，可催化大多數(shù)脂肪酸（Fulda et al.，2002）;LACS8激活長(zhǎng)鏈脂肪酸以合成細(xì)胞脂質(zhì)和通過(guò)β-氧化降解（田亞英等，2018）。AAE16為?；罨福赡軈⑴c脂肪酸向?；d體蛋白的活化。此外，本研究中有1個(gè)差異表達(dá)基因（TRINITY_DN20589_c1_g7）在數(shù)據(jù)庫(kù)中注釋為CER10。CER10為烯酰輔酶 A 還原酶，參與角質(zhì)層蠟、貯藏脂和鞘脂代謝所需的所有超長(zhǎng)鏈脂肪酸（VLCFA）延伸反應(yīng)，該基因突變體表現(xiàn)出異常的器官形態(tài)和莖光澤度（Rashotte et al.，2001）。本研究中有2個(gè)差異表達(dá)基因（TRINITY_DN19765_c1_g2和TRINITY_DN19127_c1_g5）分別注釋FATB和FATA。FATB和FATA可催化16～18個(gè)C的ACP形成16～18個(gè)C的脂肪酸，從而參與植物蠟質(zhì)生物合成過(guò)程，F(xiàn)ATB功能缺失的擬南芥葉片和莖部的蠟質(zhì)含量會(huì)降低約20%和50% （Bonaventure et al.，2003）。本研究中4個(gè)差異表達(dá)基因注釋為KCS，其在數(shù)據(jù)庫(kù)中分別注釋為KCS7、KCS9、KCS16和KCS21。前人研究發(fā)現(xiàn)，KCS催化丙二酰CoA和長(zhǎng)鏈?；鵆oA的聚合反應(yīng)，具有嚴(yán)格的底物特異性，目前已在擬南芥基因組中注釋了21個(gè)KCS基因（Joubès et al.，2008），現(xiàn)已有7個(gè)KCS基因的功能及其編碼蛋白的底物特異性被確認(rèn)，其中已確認(rèn)與蠟質(zhì)的生物合成相關(guān)基因包括KCS1、KCS2、KCS6、KCS9和KCS20（Franke et al.，2010;Lee et al.，2010;Kim et al.，2013），其他KCS的功能未見(jiàn)報(bào)道。本研究發(fā)現(xiàn)，7個(gè)差異表達(dá)基因注釋到角質(zhì)、木栓質(zhì)和蠟質(zhì)生物合成通路，其中與蠟質(zhì)合成相關(guān)的酶包括CER1、FAR、CYP86A4S和CYP86B1。其中，CER1基因編碼醛脫羰酶，在烷烴生物合成通路中催化醛脫羰形成烷烴，其突變體中烷烴的含量顯著減少，而過(guò)表達(dá) CER1 基因時(shí)，擬南芥中烷烴的含量會(huì)增加，器官呈現(xiàn)蠟粉合成減少（Bernard et al.，2012;Oshima et al.，2013）。目前在擬南芥中已檢測(cè)出8個(gè)FAR編碼基因，其中只有CER4基因參與了蠟質(zhì)的合成（Doan et al.，2009），Liu等（2018）也研究確定甘藍(lán)蠟粉缺失突變體LD10GL是由擬南芥CER4的同源基因Bol013612突變?cè)斐傻?。CYP86A4S和CYP86B1是細(xì)胞色素P450單加氧酶CYP第86家族中A和B亞族成員，均在木栓質(zhì)脂肪族成分的ω-羥基化中發(fā)揮重要作用（Vincent et al.，2009;Bruckhoff et al.，2016）。由上述推斷蠟粉合成過(guò)程中相關(guān)酶基因的差異表達(dá)是調(diào)控青花菜蠟粉合成的關(guān)鍵。BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

本研究對(duì)WT20和MT20植株中脂肪酸及蠟粉合成相關(guān)基因的表達(dá)情況進(jìn)行分析，結(jié)果發(fā)現(xiàn)20個(gè)差異表達(dá)基因的FPKM值在WT20和MT20中存在明顯差異，其中3個(gè)差異表達(dá)基因（TRINITY_DN18550_ c2_g6、TRINITY_DN10321_c0_g1TRINITY_DN210 90_c3_g4）的FPKM值在WT20中為0，1個(gè)差異表達(dá)基因（TRINITY_DN16904_c0_g1）的FPKM值在MT20值為0，差異極顯著，存在顯著和極顯著差異的差異表達(dá)基因在野生型和蠟粉缺失突變體中表達(dá)具有特異性，可作為下一步研究的重點(diǎn)。

4 結(jié)論

轉(zhuǎn)錄因子MYB家族在調(diào)控青花菜蠟粉合成中發(fā)揮重要作用。蠟粉合成過(guò)程中相關(guān)酶基因的差異表達(dá)是調(diào)控青花菜蠟粉合成的關(guān)鍵，尤其是野生型和蠟粉缺失突變體中特異性表達(dá)的差異表達(dá)基因，可作為后續(xù)研究青花菜花球表面蠟粉形成分子機(jī)制的對(duì)象。

參考文獻(xiàn)：

林琿，薛珠政，李永平，李大忠，劉建汀，朱海生，溫慶放. 2018. 青?；ㄒ撕桶坠；ㄒD(zhuǎn)錄組分析[J]. 核農(nóng)學(xué)報(bào)，32（9）：1708-1720. [Lin H，Xue Z Z，Li Y P，Li D Z，Liu J T，Zhu H S，Wen Q F. 2018. Transcriptome analysis of the green stalk and the white stalk of cauliflower[J]. Journal of Nuclear Agricultural Sciences，32（9）：1708-1720.] doi：10.11869/j.issn.100-8551.2018.09.1708.

田亞英，傅念，吳清. 2018. 長(zhǎng)鏈脂酰輔酶A合成酶及其調(diào)控因素的研究及進(jìn)展[J]. 醫(yī)學(xué)研究雜志，47（10）：13-15. [Tian Y Y，F(xiàn)u N，Wu Q. 2018. Research and progress of long-chain acyl-CoA synthase and its regulatory factors[J]. Journal of Medical Research，47（10）：13-15.] doi：10.11969/j.issn.1673-548X.2018.10.005.

王立山，丁兵，李玉花，張旸. 2018. 植物表皮蠟質(zhì)合成轉(zhuǎn)運(yùn)調(diào)控相關(guān)基因與干旱響應(yīng)的研究進(jìn)展[J]. 園藝學(xué)報(bào)，45（9）：1831-1843. [Wang L S，Ding B，Li Y H，Zhang Y. 2018. Reaserch progress of plant cuticular wax biosynthesis，export and regulation related genes responsed to drought[J]. Acta Horticulturae Sinica，45（9）：1831-1843.] doi：10.16420/j.issn.0513-353x.2017-0548.

吳寧柔，韋云婷，官梅. 2018. 甘藍(lán)型油菜不同花色轉(zhuǎn)錄組分析[J]. 分子植物育種，16（22）：7234-7240. [Wu N R，Wei Y T，Guan M. 2018. Transcriptome analysis of Brassica napus with different flower colors[J]. Molecular Plant Breeding，16（22）：7234-7240.] doi：10.13271/j.mpb.016. 007234.

張曦，王秋實(shí)，鄒春蕾，劉志勇，王一衡，馮輝. 2013. 大白菜花莖蠟粉基因的遺傳分析與初步定位[J]. 分子植物育種，11 （6）：804-808. [Zhang X，Wang Q S，Zou C L，Liu Z Y，Wang Y H，F(xiàn)eng H. 2013. Genetic analysis and preliminary mapping of wax gene on stem in Chinese cabbage[J]. Molecular Plant Breeding，11（6）：804-808.] doi：10. 3969/mpb.011.000804.

張振超，潘永飛，戴忠良，秦文斌，山溪. 2021. 青花菜越冬設(shè)施高產(chǎn)栽培技術(shù)[J]. 長(zhǎng)江蔬菜，（11）：24-26. [Zhang Z C，Pan Y F，Dai Z L，Qin W B，Shan X. 2021. High-yield cultivation techniques of broccoli overwintering in facility[J]. Journal of Changjiang Vegetables，（11）：24-26.]

張振超，姚悅梅，毛忠良，孫國(guó)勝，秦文斌，戴忠良. 2018. 基于高通量測(cè)序的青花菜早期發(fā)育小孢子轉(zhuǎn)錄組分析與基因功能注釋[J]. 核農(nóng)學(xué)報(bào)， 32（5）：848-855. [Zhang Z C，Yao Y M，Mao Z L，Sun G S，Qin W B，Dai Z L. 2018. Transcriptome analysis and gene function annotation of early developmental broccoli microspores based on high-throughput sequencing technology[J]. Journal of Nuclear Agricultural Sciences，32（5）：848-855.] doi： 10.11869/j.issn.100-8551.2018.05.0848.BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

Bernard A，Domergue F，Pascal S，Jetter R，Renne C，F(xiàn)aure J D，Haslam R P，Napier J A，Lessire R，Joubes J. 2012. Reconstitution of plant alkane biosynthesis in yeast de-monstrates that Arabidopsis ECERIFERUM1 and ECERIFERUM3 are core components of a very-long-chain alkane synthesis complex[J]. The Plant Cell，24（7）：3106-3118. doi：10.1105/tpc.112.099796.

Bonaventure B，Salas J J，Pollard M R，Ohlrogge J B. 2003. Disruption of the FATB gene in Arabidopsis demonstrates an essential role of saturated fatty acids in plant growth[J]. The Plant Cell，15（4）：1020-1033. doi：10.1105/tpc. 008946.

Bourdenx B，Bernard A，Domergue F，Pascal S，Léger A，Roby D，Pervent M，Vile D，Haslam R，Napier J A，Lessire R，Joubès J. 2011. Overexpression of Arabidopsis ECERIFERUM1 promotes wax very-long-chain alkane biosynthesis and influences plant response to biotic and abiotic stresses[J]. Plant Physiology，156（1）：29-45. doi：10.1104/ pp.111.172320.

Bruckhoff V，Haroth S，F(xiàn)eussner K，Konig S，Brodhun F，F(xiàn)eussner I. 2016. Functional characterization of CYP94-genes and identification of a novel jasmonate catabolite in flowers[J]. PLoS One，11（7）：e0159875. doi：10.1371/journal.pone.0159875.

Wettstein D V，Chua N H. 1987. Plant molecular biology[M]. New York ：Plenum Press，305-314.

Chen X B，Goodwin M，Boroff V L，Liu X L，Jenks M A. 2003. Cloning and characterization of the WAX2 gene of Arabidopsis involved in cuticle membrane and wax production[J]. The Plant Cell，15（5）：1170-1185. doi：10. 1105/tpc.010926.

Doan T T，Carlsson A S，Hamberg M，Bülow L，Stymne S，Olsson P. 2009. Functional expression of five Arabidopsis fatty acyl-CoA reductase genes in Escherichia coli[J]. Journal of Plant Physiology，166（8）：787-796. doi：10.1016/ j.jplph.2008.10.003.

Domínguez E，Heredia-Guerrero J A，Heredia A. 2011. The biophysical design of plant cuticles：An overview[J]. New Phytologist，189（4）：938-949. doi：10.1111/j.1469-8137. 2010.03553.

Dubos C，Stracke R，Grotewold E，Weisshaar B，Martin C，Lepiniec L. 2010. MYB transcription factors in Arabidopsis[J]. Trends in Plant Science，15（10）：573-581. doi：10. 1016/j.tplants.2010.06.005.

Franke R，Hfer R，Briesen I，Emsermann M，Schreiber L. 2010. The DAISY gene from Arabidopsis encodes a fatty acid elongase condensing enzyme involved in the biosynthesis of aliphatic suberin in roots and the chalaza-micropyle region of seeds[J]. The Plant Journal，57（1）：80-95. doi：10.1111/j.1365-313X.2008.03674.x.

Fulda M，Shockey J，Werber M，Wolter F P，Heinz E. 2002. Two long-chain acyl-CoA synthetases from Arabidopsis thaliana involved in peroxisomal fatty acid beta-oxidation[J]. The Plant Journal，32（1）：93-103. doi：10.1046/j.1365- 313x.2002.01405.xBF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

Jenks M A，Eigenbrode S D，Lemieux B. 2002. Cuticular waxes of Arabidopsis[J]. Arabidopsis Book，1：e0016. doi：10. 1199/tab.0016.

Joubès J，Raffaele S，Bourdenx B，Garcia C，Laroche-Traineau J，Moreau P，Domergue F，Lessire R. 2008. The VLCFA elongase gene family in Arabidopsis thaliana：Phylogene-tic analysis， 3D modelling and expression profiling[J]. Plant Molecular Biology，67（5）：547. doi：10.1007/s11103- 008-9339-z.

Kim J，Jung J H，Lee S B，Go Y S，Kim H J，Cahoon R，Markham J E，Cahoon E B，Suh M C. 2013. Arabidopsis 3-ketoacyl-coenzyme a synthase9 is involved in the synthesis of tetracosanoic acids as precursors of cuticular waxes，suberins， sphingolipids，and phospholipids[J]. Plant Physiology，162（2）：567-580. doi：10.1104/pp.112.210450.

Kunst L，Samuels A L. 2003. Biosynthesis and secretion of plant cuticular wax[J]. Progress in Lipid Research，42（1）：51-80. doi：10.1016/s0163-7827（02）00045-0.

Lee S B，Jung S J，Go Y S，Kim H U，Kim J K，Cho H J，Park O K，Suh M C. 2010. Two Arabidopsis 3-ketoacyl CoA synthase genes， KCS20 and KCS2/DAISY， are functionally redundant in cuticular wax and root suberin biosynthesis，but differentially controlled by osmotic stress[J]. The Plant Journal，60（3）：462-475. doi：10.1111/j. 1365-313X.2009.03973.x.

Liu D M，Tang J，Liu Z Z，Dong X，Zhuang M，Zhang Y Y，Lü H H，Sun P T，Liu Y M，Li Z S，Ye Z B，F(xiàn)ang Z Y，Yang L M. 2017. Cgl2 plays an essential role in cuticular wax biosynthesis in cabbage（Brassica oleracea L. var. capitata）[J]. BMC Plant Biology，17（1）：223. doi：10. 1186/s12870-017-1162-8.

Liu D M，Dong X，Liu Z Z，Tang J，Zhuang M，Zhang Y Y，Lv H H，Liu Y M，Li Z S，F(xiàn)ang Z Y，Yang L M. 2018. Fine mapping and candidate gene identification for wax biosynthesis locus， BoWax1 in Brassica oleracea L. var. capitata[J]. Frontiers in Plant Science，9：309. doi：10.3389/fpls.2018.00309.

Lü S，Song T，Kosma D K，Parsons E P，Rowland O，Jenks M A. 2009. Arabidopsis CER8 encodes LONG-CHAIN ACYL COA SYNTHETASE 1（LACS1） and has overlapping functions with LACS2 in plant wax and cutin synthesis[J]. The Plant Journal，59（4）：553-564. doi：10.1111/ j.1365-313X.2009.03892.

Oliveira A F，Meirelles S T，Salatino A. 2004. Epicuticular waxes from caatinga and cerrado species and their efficiency against water loss[J]. Anais da Academia Brasileira de Ciências，75（4）：431-439. doi：10.1590/s0001-3765200 3000400003.

Oshima Y，Shikata M，Koyama T，Ohtsubo N，Ohme-Takagi M M. 2013. MIXTA-Like transcription factors and WAX INDUCER1/SHINE1 coordinately regulate cuticle deve-lopment in Arabidopsis and Torenia fournieri[J]. The Plant Cell，25（5）：1609-1624. doi：10.1105/tpc.113.110783.BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0

Pruitt R E，Vielle-Calzada J P，Ploense S E，Grossniklaus U，Lolle S J. 2000. FIDDLEHEAD， a gene required to suppress epidermal cell interactions in Arabidopsis， encodes a putative lipid biosynthetic enzyme[J]. Proceedings of the National Academy of Sciences，97（3）：1311-1316. doi：10.1073/pnas.97.3.1311.

Pu Y Y，Gao J，Guo Y L，Liu T T，Zhu L X，Xu P，Yi B，Wen J，Tu J X，Ma C Z，F(xiàn)u T D，Zou J T，Shen J X. 2013. A novel dominant glossy mutation causes suppression of wax biosynthesis pathway and deficiency of cuticular wax in Brassica napus[J]. BMC Plant Biology，13：215-228. doi：10.1186/1471-2229-13-215.

Pulsifer I P，Kluge S，Rowland O. 2012. Arabidopsis LONG-CHAIN ACYL-COA SYNTHETASE 1（LACS1）， LACS2，and LACS3 facilitate fatty acid uptake in yeast[J]. Plant Physiology and Biochemistry，51：31-39. doi：10.1016/j.plaphy.2011.10.003.

Pyee J，Kolattukudy P E. 1995. The gene for the major cuticular wax-associated protein and three homologous genes from broccoli（Brassica oleracea） and their expression patterns[J]. The Plant Journal，7（1）：49-59. doi：10.1046/j.1365-313x.1995.07010049.

Rashotte A M，Jenks M A，F(xiàn)eldmann K A. 2001. Cuticular waxes on eceriferum mutants of Arabidopsis thaliana[J]. Phytochemistry，57（1）：115-123. doi：10.1016/s0031-9422（00）00513-6.

Saet B L，Mi C S. 2015. Advances in the understanding of cuticular waxes in Arabidopsis thaliana and crop species[J]. Plant Cell Reports，34（4）：557-572. doi：10.1007/s00299-015-1772-2.

Vincent C，Patrik D，Irène B，Denise M，Hubert S，Lukas S，Rochus F，F(xiàn)ranck P. 2009. CYP86B1 is required for very long chain ω-hydroxyacidand a，ω-dicarboxylic acid synthesis in root and seed suberin polyester[J]. Plant Physio-logy，150：1831-1843. doi：10.1104/pp.109.141408.

Zhang X，Liu Z，Wang P，Wang Q，Yang S，F(xiàn)eng H. 2013. Fine mapping of BrWax1， a gene controlling cuticular wax biosynthesis in Chinese cabbage（Brassica rapa L. ssp. pekinensis）[J]. Molecular Breeding，32（4）：867-874. doi：10.1007/s11032-013-9914-0.

（責(zé)任編輯 陳 燕）BF9E10E1-A8E8-4AA4-96C2-0B323DA83DF0