關(guān)鍵詞：文心蘭；唇瓣化；轉(zhuǎn)錄組測(cè)序；差異表達(dá)基因

中圖分類號(hào)：S432.1 文獻(xiàn)標(biāo)志碼：A

蘭科植物因其獨(dú)特的花型倍受人們喜愛，MADS-box 轉(zhuǎn)錄因子是花的起始和發(fā)育過程中的關(guān)鍵調(diào)控成分，在參與花器官的發(fā)育、開花時(shí)間的調(diào)節(jié)等過程中均起著重要作用[1]。ABCDE 模型由MADS-box 基因家族中的五類同源基因組成，這些基因通過相互作用共同決定花器官身份[2]。因此，大多數(shù)關(guān)于蘭花花器官發(fā)育的研究都集中在這些ABCDE 基因上，而關(guān)于調(diào)控蘭花唇瓣發(fā)育的其他轉(zhuǎn)錄因子及下游相關(guān)調(diào)控網(wǎng)絡(luò)仍不清晰，因此，探究其他相關(guān)調(diào)控機(jī)制對(duì)蘭花植物分子育種具有重要意義。

文心蘭是蘭科文心蘭屬植物，因其獨(dú)特的花型和艷麗的顏色是世界上具有較高觀賞價(jià)值的切花品種[3]。在長期的文心蘭生長進(jìn)化過程中，極易產(chǎn)生花型多樣性，這種花型多樣性主要體現(xiàn)在唇瓣的變化上，使其具有更高的觀賞價(jià)值和生物學(xué)價(jià)值[4]。因此揭示文心蘭唇瓣多樣性形成機(jī)制顯得極其重要。“花被密碼模型”是目前被廣泛接受的蘭花花器官發(fā)育模型， 該模型中， 不同MADS-box 基因分別組成唇瓣復(fù)合物和萼片/花瓣復(fù)合物，這2 種蛋白復(fù)合物相互競爭決定蘭花花器官身份形成[5]，而文心蘭花瓣形態(tài)結(jié)構(gòu)多樣化的其他調(diào)控網(wǎng)絡(luò)研究仍鮮有報(bào)道。本研究以檸檬綠文心蘭（Oncidium. flexuosum‘Honey Angel’）和其側(cè)花瓣唇瓣化變異種（突變體）為材料，對(duì)2 種側(cè)花瓣進(jìn)行RNA 轉(zhuǎn)錄組測(cè)序， 通過MADS-box 及其他轉(zhuǎn)錄因子的差異表達(dá)挖掘調(diào)控文心蘭唇瓣化變異關(guān)鍵轉(zhuǎn)錄因子，為建立文心蘭唇瓣化變異調(diào)控關(guān)系網(wǎng)絡(luò)，培育花型獨(dú)特蘭花新品種奠定基礎(chǔ)。

1 材料與方法

1.1 材料

供試材料為種植在華南農(nóng)業(yè)大學(xué)花卉研究中心種質(zhì)資源圃中生長健壯的檸檬綠文心蘭側(cè)花瓣發(fā)生唇瓣化變異種（突變體），對(duì)照為檸檬綠文心蘭（圖1）。對(duì)2種側(cè)花瓣表皮細(xì)胞觀察，發(fā)現(xiàn)在8 mm 花蕾時(shí)期表皮細(xì)胞數(shù)目差異最大，因此在9月和10 月采集8 mm 花蕾并分離其側(cè)花瓣迅速置于液氮中凍存，然后保存于–80 ℃冰箱，用于后續(xù)轉(zhuǎn)錄組測(cè)序及相關(guān)試驗(yàn)。

1.2 方法

1.2.1 總RNA 提取、文庫構(gòu)建與轉(zhuǎn)錄組測(cè)序（1）總RNA 提取。使用天根多糖多酚試劑盒（QIAGEN， Germany）提取檸檬綠文心蘭和唇瓣化變異種的側(cè)花瓣組織樣品的總RNA，并檢測(cè)其總量和完整性。RNA 提取與檢測(cè)均委托北京諾禾致源科技股份有限公司完成。

（2）文庫構(gòu)建、質(zhì)檢與測(cè)序。選取符合質(zhì)量標(biāo)準(zhǔn)的總RNA 通過Oligo（dT）磁珠富集帶有polyA 尾的mRNA。隨后在NEB FragmentationBuffer 中用二價(jià)陽離子將得到的mRNA 隨機(jī)打斷，按照NEBNext? Ultra? RNA Library Prep Kitfor Illumina?試劑盒說明書進(jìn)行cDNA 文庫的構(gòu)建。文庫構(gòu)建完成后，使用Qubit 2.0 Fluorometer熒光分光光度計(jì)進(jìn)行初步定量，隨后使用Agilent2100 bioanalyzer 生物分析儀對(duì)文庫的insert size進(jìn)行檢測(cè)，通過RT-qPCR 對(duì)文庫有效濃度進(jìn)行準(zhǔn)確定量，以保證文庫質(zhì)量。把不同文庫按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling 后，使用Illumina NovaSeq 6000（Illumina，USA）測(cè)序平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序。

1.2.2 測(cè)序數(shù)據(jù)分析 （1）數(shù)據(jù)質(zhì)控與轉(zhuǎn)錄本拼接。為了保證數(shù)據(jù)分析的質(zhì)量和可靠性，將測(cè)序獲得的原始數(shù)據(jù)（raw reads）進(jìn)行過濾，去除帶接頭reads、含N 堿基的reads、低質(zhì)量reads（Q 值lt;20），經(jīng)過原始數(shù)據(jù)過濾、測(cè)序錯(cuò)誤率檢查、GC 含量分布檢查，獲得后續(xù)分析使用的高質(zhì)量clean reads。使用Trinity 程序?qū)lean reads進(jìn)行拼接， 對(duì)拼接獲得的轉(zhuǎn)錄本去除冗余以Corset 層次聚類后得到最長Cluster 序列進(jìn)行后續(xù)分析。

（2）unigene 功能注釋。得到拼接轉(zhuǎn)錄本后，將轉(zhuǎn)錄本通過NCBI Blastx 后比對(duì)到蛋白質(zhì)的7 個(gè)數(shù)據(jù)庫[Nr、Nt、Pfam、KOG、Swissprot、KEGG、GO，（E-value≤0.000001）]，獲得與目標(biāo)序列相似度最高的蛋白質(zhì)，并對(duì)其unigenes 進(jìn)行功能注釋。

1.2.3 差異表達(dá)基因篩選及富集分析 利用RSEM（RNA-seq by Expectation Maximization）軟件[6]對(duì)基因在不同樣本之間進(jìn)行定量分析，并使用DESeq2 軟件對(duì)原始的read count 進(jìn)行標(biāo)準(zhǔn)化處理， 通過統(tǒng)計(jì)學(xué)模型進(jìn)行假設(shè)檢驗(yàn)概率（P-value）計(jì)算后進(jìn)行FDR 矯正獲得padj值。為了得到顯著差異表達(dá)基因，將篩選條件設(shè)為padjlt;0.005 且差異倍數(shù)|log2（fold change）|gt;1。分別使用GOseq[7]和KOBAS 軟件進(jìn)行差異表達(dá)基因的GO 和KEGG 的富集分析，當(dāng)padjlt;0.05 時(shí)，為顯著富集。

1.2.4 差異表達(dá)基因的RT-qPCR 驗(yàn)證 以檸檬綠文心蘭及其唇瓣化變異種為材料，基于轉(zhuǎn)錄組數(shù)據(jù)分析，篩選出與唇瓣變異相關(guān)的差異轉(zhuǎn)錄因子基因進(jìn)行RT-qPCR 驗(yàn)證。提取8 mm 花蕾期檸檬綠和唇瓣化變異種文心蘭的側(cè)花瓣的RNA，利用湖南艾克瑞生物工程有限公司的Evo M-mlV反轉(zhuǎn)錄試劑盒合成cDNA 單鏈。利用Primer Premier5.0 軟件設(shè)計(jì)引物，委托上海生物工程有限公司合成引物。RT-qPCR 反應(yīng)體系和程序按照Hieff?qPCR SYBR Green Master Mix（Low Rox Plus）試劑盒操作說明進(jìn)行，參考文心蘭屬植物相關(guān)研究以Actin 為內(nèi)參基因[8]，通過Applied BiosystemsABI7500 熒光定量PCR 儀檢測(cè)10 個(gè)相關(guān)轉(zhuǎn)錄因子基因的相對(duì)表達(dá)量，每個(gè)樣品設(shè)置3 個(gè)生物學(xué)重復(fù)，反應(yīng)結(jié)束后確認(rèn)擴(kuò)增曲線和熔解曲線，使用2-ΔΔCT 法計(jì)算相關(guān)基因在不同樣品中的表達(dá)情況。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測(cè)序質(zhì)量分析

取檸檬綠文心蘭和唇瓣化變異種8 mm 花蕾時(shí)期側(cè)花瓣進(jìn)行轉(zhuǎn)錄組測(cè)序。測(cè)序共獲得46 436 745個(gè)raw reads，原始數(shù)據(jù)過濾后獲得45 525 415 個(gè)高質(zhì)量clean reads。經(jīng)質(zhì)量評(píng)估，數(shù)據(jù)整體測(cè)序錯(cuò)誤率不超過0.03%，Q20 堿基百分比（測(cè)序錯(cuò)誤率不小于1%）分布在96.61%～98%之間，Q30（測(cè)序錯(cuò)誤率小于0.1%）介于91.63%～94.25%之間，GC 含量介于45.92%～46.67%之間，表明測(cè)序質(zhì)量較好，可用于后期數(shù)據(jù)分析。

2.2 轉(zhuǎn)錄組unigenes 功能注釋

將高質(zhì)量clean reads 進(jìn)行拼接后共獲得52253 個(gè)unigenes，利用NR、Swissport、PFAM、GO、NT、KO、KOG 等7 個(gè)常用蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對(duì)及基因功能注釋，根據(jù)比對(duì)蛋白數(shù)據(jù)庫的序列同源性結(jié)果可知，至少在1 個(gè)數(shù)據(jù)庫注釋成功的unigenes 有36 118，占69.12%，其中與NR數(shù)據(jù)庫匹配的unigenes 數(shù)目最多為31 828 條，占60.91%；其次為Swissport 數(shù)據(jù)庫23 174 條（占44.34%），PFAM 數(shù)據(jù)庫21 369 條（占40.89%），GO 數(shù)據(jù)庫21 367 條（占40.89%），NT 數(shù)據(jù)庫20830 條（占39.86%），KO 數(shù)據(jù)庫10 523 條（占20.13%），KOG 數(shù)據(jù)庫6465 條（占12.37%）；在以上7 個(gè)數(shù)據(jù)庫中都注釋成功的unigenes 數(shù)目為3549 條，占6.79%。

2.3 文心蘭側(cè)花瓣差異表達(dá)基因分布

對(duì)2 種側(cè)花瓣的unigenes 表達(dá)量進(jìn)行比較，以|log2（fold change）|≥1 且padjlt;0.05 作為篩選基準(zhǔn)，在2 種側(cè)花瓣中共檢測(cè)到差異表達(dá)基因（differentially expressed gene， DEG）3120 條，其中與檸檬綠文心蘭相比，在突變體側(cè)花瓣中上調(diào)表達(dá)的DEG 有1526 條，下調(diào)表達(dá)的DEG 有1594條（圖2）。

2.4 差異表達(dá)基因的GO功能分析

將檢測(cè)到的DEG 進(jìn)行GO 功能富集分析，以注釋其生物學(xué)功能。共有1537 條差異表達(dá)基因注釋到細(xì)胞組分（cellular component， CC）、生物學(xué)過程（ biological process， BP ） 和分子功能（molecular function， MF）3 個(gè)功能群組。其中差異表達(dá)基因更多富集在生物學(xué)過程功能中，主要途徑包括生物學(xué)過程功能中的碳水化合物代謝過程（carbohydrate metabolic process）、脂質(zhì)代謝過程（lipid metabolic process）和含核堿基化合物分解代謝過程（ nucleobase-containing compoundcatabolic process）；在細(xì)胞組分功能上，主要富集在核質(zhì)（nucleoplasm）和染色體（chromosome）中；分子功能中差異基因重點(diǎn)富集在氧化還原酶活性（oxidoreductase activity）和裂解酶活性（lyaseactivity）上（圖3）。

2.5 差異表達(dá)基因的KEGG 富集分析

將差異表達(dá)的基因進(jìn)行KEGG 富集分析，結(jié)果發(fā)現(xiàn)在2 種側(cè)花瓣中有690 個(gè)差異表達(dá)的基因在KEGG 途徑中顯著富集，KEGG 富集的前20個(gè)途徑如圖4 所示，其中代謝途徑（metabolicpathways）、次生代謝物的生物合成（biosynthesisof secondary metabolites ） 和抗生素生物合成（biosynthesis of antibiotics）是富集差異表達(dá)基因個(gè)數(shù)較多的通路。

2.6 文心蘭唇瓣化變異相關(guān)基因的挖掘與分析

MADS-box 基因是一類廣泛存在于植物中且被證明在蘭花花器官發(fā)育中具有極其重要的調(diào)控作用。為了挖掘文心蘭唇瓣化變異相關(guān)調(diào)控基因，首先從差異表達(dá)基因中重點(diǎn)關(guān)注MADS-box轉(zhuǎn)錄因子，在該轉(zhuǎn)錄組中共挖掘5 個(gè)存在顯著差異（ 以|log2（fold change）|gt;1 為篩選依據(jù)） 的MADS-box 基因（MADS1、MADS10、MADS11、MADS7、MADS5）。在植物中MADS-box 基因主要通過調(diào)節(jié)其他基因的轉(zhuǎn)錄來發(fā)揮作用[9]。因此，也對(duì)文心蘭的側(cè)花瓣轉(zhuǎn)錄組測(cè)序結(jié)果進(jìn)行其他轉(zhuǎn)錄因子的數(shù)據(jù)分析，發(fā)現(xiàn)存在差異表達(dá)的其他轉(zhuǎn)錄因子主要包括WRKY、TCP、NAC、MYB等（表1）。

將MADS1、MADS10、MADS11、MADS7、MADS5、TCP1、TCP2、NAC1、MYB1 蛋白與其他物種同源蛋白序列構(gòu)建進(jìn)化樹，結(jié)果表明，MADS-box 蛋白與扇葉文心蘭（Erycina pusilla）的MADS-box 親緣關(guān)系最近，聚為一類，其次與鐵皮石斛（Dendrobium catenatum）親緣關(guān)系較近；NAC 與鐵皮石斛（Dendrobium catenatum）、蝴蝶蘭（Phalaenopsis equestris）親緣關(guān)系較近，聚為一類；TCP 與蝴蝶蘭（Phalaenopsis equestris）、蕙蘭（Cymbidium faberi）親緣關(guān)系最近，聚為一類；MYB 與建蘭（Cymbidium ensifolium）親緣關(guān)系最近，聚為一類，其次與扇葉文心蘭（Erycinapusilla）、鐵皮石斛（Dendrobium catenatum）親緣關(guān)系較近（圖5）。

2.7 轉(zhuǎn)錄組測(cè)序RT-qPCR驗(yàn)證

為了驗(yàn)證文心蘭側(cè)花瓣轉(zhuǎn)錄組測(cè)序數(shù)據(jù)的真實(shí)性，共選擇10 個(gè)差異表達(dá)基因，通過RT-qPCR的方法檢測(cè)其在花蕾期階段的側(cè)花瓣中的表達(dá)水平，與轉(zhuǎn)錄組測(cè)序獲得的基因表達(dá)水平相比，二者變化趨勢(shì)一致（圖6），說明本研究的轉(zhuǎn)錄組測(cè)序結(jié)果具有良好的可靠性。

3 討論

本研究對(duì)檸檬綠文心蘭和其側(cè)花瓣唇瓣化變異種花蕾期側(cè)花瓣進(jìn)行RNA 測(cè)序和轉(zhuǎn)錄組分析，共獲得52 253 個(gè)unigenes。利用NR、Swissprot、PFAM、GO、NT、KO、KOG 七個(gè)數(shù)據(jù)庫分別進(jìn)行基因功能注釋，分別注釋到31 828（60.91%）、23 174（44.34%）、21 369（40.89%）、21367（40.89%）、20 830（39.86%）、10 523（20.13%）、6465（12.37%）條unigenes，說明本轉(zhuǎn)錄組測(cè)序數(shù)據(jù)質(zhì)量較高，為文心蘭唇瓣發(fā)育相關(guān)研究及其分子育種提供理論基礎(chǔ)。

檢測(cè)到差異表達(dá)基因3120 條，與檸檬綠文心蘭相比，在突變體側(cè)花瓣中上調(diào)基因有1526 條，下調(diào)基因有1594 條。其中1537 條差異表達(dá)基因經(jīng)GO 功能富集分析注釋到生物過程、細(xì)胞組分和分子功能3 個(gè)功能群組，經(jīng)KEGG 功能富集主要分布在代謝途徑、次生代謝物生物合成途徑等。表明文心蘭唇瓣化變異是一個(gè)復(fù)雜的過程，涉及多條代謝通路并受多種轉(zhuǎn)錄因子基因共同協(xié)調(diào)影響。

現(xiàn)有的研究揭示了MADS-box 基因尤其是B類MADS-box 基因在蘭花花器官中的重要作用[10]。同時(shí)，同源性蘭花Tepal（HOT）模型可以解釋蘭花MADS-box 基因如何調(diào)控花器官的發(fā)育并推測(cè)B 類基因中的不同基因與其他類別的MADS-box基因結(jié)合，以調(diào)節(jié)萼片、側(cè)花瓣和唇瓣身份的復(fù)雜性[11]。截至目前，B 類MADS-box 基因通過基因復(fù)制事件產(chǎn)生大量同源基因，有學(xué)者從南茜文心蘭中分離到4 個(gè)B 類MADS-box 基因（OMADS3、OMADS5、OMADS8 和OMADS9），這些基因相互作用形成同源或異源二聚體共同參與調(diào)控花被器官的形成，當(dāng)OMADS5 缺失而OMADS9 存在時(shí)會(huì)導(dǎo)致大唇瓣以及萼片/側(cè)花瓣唇瓣化突變體的產(chǎn)生[12]。除了B 類MADS-box 基因，“花被密碼”模型認(rèn)為B 類AP3 基因與E 類AGL6 基因相互組成萼片/側(cè)花瓣復(fù)合物和唇瓣復(fù)合物，這2 種蛋白復(fù)合物相互競爭調(diào)控著蘭花復(fù)雜花被器官的形成[5]。在蝴蝶蘭唇狀花瓣突變體中，B 類MADS-box 基因（PeMADS4）在3 個(gè)唇瓣?duì)罨ò曛芯哂忻黠@的表達(dá)差異[13]，說明B 類MADS-box 基因在蝴蝶蘭唇瓣發(fā)育中具有重要調(diào)控作用。此外，建蘭花器官類型非常多樣，并通過形態(tài)學(xué)觀察和轉(zhuǎn)錄組分析發(fā)現(xiàn)其花型突變體與MADS-box 基因的異常表達(dá)有關(guān)。在建蘭唇狀花瓣突變體以及唇狀萼片突變體中，CeAP3-3、CeAP3-4 和CeAGL6-2 基因的高表達(dá)導(dǎo)致了這種唇瓣突變表型[14]。上述研究結(jié)果進(jìn)一步證實(shí)了“花被密碼”模型及該模型在蘭科植物的廣泛適應(yīng)性。本研究通過對(duì)檸檬綠文心蘭和其唇瓣化變異種的側(cè)花瓣進(jìn)行轉(zhuǎn)錄組測(cè)序，從差異基因中篩選到5 個(gè)差異表達(dá)MADS-box 基因MADS1、MADS10、MADS11、MADS7、MADS5，其中MADS5為B 類基因且在突變體側(cè)花瓣中下調(diào)表達(dá)，這與在南茜文心蘭中該基因負(fù)向調(diào)控唇瓣器官的形成結(jié)果一致，推測(cè)其可能參與組成文心蘭萼片/側(cè)花瓣復(fù)合物促進(jìn)檸檬綠文心蘭正常花瓣的形成。MADS1、MADS10、MADS11、MADS7 為E 類基因且在突變體側(cè)花瓣中上調(diào)表達(dá)（除MADS7），大量研究已經(jīng)證實(shí)MADS1（AGL6-2）參與唇瓣復(fù)合物組成并促進(jìn)文心蘭大唇瓣、蝴蝶蘭唇瓣?duì)罨ò晖蛔凅w及建蘭花瓣/萼片唇瓣化突變體的產(chǎn)生[13-15]，其他E 類基因是否同樣參與唇瓣復(fù)合物的組成，或與其他調(diào)控元件一起共同參與文心蘭唇瓣器官的發(fā)育還需要通過大量實(shí)驗(yàn)進(jìn)一步證實(shí)。

盡管已有研究表明，MADS-box 轉(zhuǎn)錄因子在蘭花花瓣及唇瓣多樣性的形成上扮演著重要的角色，但該基因家族是如何調(diào)控蘭花唇瓣多樣性形成機(jī)理尚不清楚?；ㄆ鞴偕矸莼虻南掠伟谢蚴腔ㄆ鞴侔l(fā)生的關(guān)鍵調(diào)節(jié)因子，直接或間接作用于花器官組織[16]。通過高通量測(cè)序數(shù)據(jù)及整合多種生物學(xué)數(shù)據(jù)可以推斷出相關(guān)的基因表達(dá)調(diào)控網(wǎng)絡(luò)。因此，本研究也對(duì)轉(zhuǎn)錄組中差異表達(dá)的其他轉(zhuǎn)錄因子進(jìn)行分析，發(fā)現(xiàn)NAC、TCP、MYB、WRKY 等轉(zhuǎn)錄因子在文心蘭唇瓣發(fā)育過程中存在差異表達(dá)。有研究發(fā)現(xiàn)，在玫瑰花瓣中NAC轉(zhuǎn)錄因子控制花朵花瓣中乙烯調(diào)節(jié)的細(xì)胞擴(kuò)張，為NAC 轉(zhuǎn)錄因子的功能提供新見解[17]。在小蘭嶼蝴蝶蘭中，挖掘到PeNAC67 是調(diào)控小蘭嶼蝴蝶蘭唇瓣發(fā)育的關(guān)鍵轉(zhuǎn)錄因子，且該基因與唇瓣發(fā)育相關(guān)的MADS3 基因存在互作，揭示了蘭花唇瓣發(fā)育的分子機(jī)制[18]。將蝴蝶蘭TCP 轉(zhuǎn)錄因子PhCYCs 基因在夏堇中進(jìn)行異源表達(dá)，發(fā)現(xiàn)蝴蝶蘭PhCYCs 基因在夏堇中能正常轉(zhuǎn)錄，且所有轉(zhuǎn)基因株系分株增多且花型改變，初步說明蝴蝶蘭TCP 轉(zhuǎn)錄因子可能既參與花型發(fā)育又調(diào)控分枝過程[19]。植物MYB 轉(zhuǎn)錄因子除了在調(diào)控植物花青素生物合成上有重要貢獻(xiàn)外，也有研究發(fā)現(xiàn)該基因家族在水稻中會(huì)調(diào)控花器官和小穗的分化[20]，而該轉(zhuǎn)錄因子在花卉植物花器官形成的相關(guān)調(diào)控作用尚未見報(bào)道。本研究通過轉(zhuǎn)錄組分析鑒定出的這些可能調(diào)控文心蘭唇瓣發(fā)育的其他轉(zhuǎn)錄因子基因是否會(huì)與MADS-box 轉(zhuǎn)錄因子相互作用共同調(diào)控文心蘭唇瓣器官的發(fā)育仍需要通過大量實(shí)驗(yàn)進(jìn)一步深入研究，以確定其具體功能和調(diào)控方式。

因此，本研究對(duì)唇瓣化變異文心蘭側(cè)花瓣進(jìn)行轉(zhuǎn)錄組分析，獲得了大量差異表達(dá)轉(zhuǎn)錄因子，為深入闡明這些轉(zhuǎn)錄因子的功能，完善文心蘭唇瓣發(fā)育相關(guān)調(diào)控網(wǎng)絡(luò)及分子機(jī)理奠定基礎(chǔ)。