瀕危植物秦嶺石蝴蝶花瓣數(shù)量變異機理研究

蔣景龍，孫 旺，李 麗，李 耘，胡鳳成

(1 陜西理工大學(xué) 生物科學(xué)與工程學(xué)院，陜西漢中 723000；2陜西理工大學(xué) 化學(xué)與環(huán)境科學(xué)學(xué)院，陜西漢中 723000； 3 漢中市野生動植物保護(hù)管理站，陜西漢中 723000；4 略陽縣苗圃，陜西漢中 724310)

花是被子植物特有的觀賞和生殖器官，是植物分類學(xué)中的重要依據(jù)，同時也是被子植物進(jìn)化過程中變化最豐富的器官[1]。研究花器官發(fā)育及其分子機制對其植株的發(fā)育及演化過程具有重要意義[2-3]。花器官發(fā)育是由相關(guān)基因經(jīng)過一系列復(fù)雜的調(diào)控機制形成的，Coen等[4]首次于1991年提出了花器官發(fā)育的經(jīng)典ABC模型以及隨著研究的深入，由ABC模型發(fā)展起來的其他模型?；谀Ｊ街参飻M南芥和金魚草各類花器官的研究，在其他物種中大量決定花器官特征屬性的基因相繼被克隆鑒定出來，任何一種基因突變都會導(dǎo)致花器官形態(tài)發(fā)生變化[5-6]。在環(huán)境因素(不同生物和非生物因素)的影響下，花器官形態(tài)常會顯現(xiàn)出特定的變異類型，其分子機制歸于基因結(jié)構(gòu)的改變和表達(dá)模式的多樣化[7-8]。對花器官變異的研究有助于揭示瀕危種群野生植物響應(yīng)環(huán)境變化的遺傳機制，對理解其致瀕機理進(jìn)而采取相應(yīng)的保護(hù)措施具有重要的意義[9]。

秦嶺石蝴蝶(PetrocosmeaqinlingensisW. T. Wang)屬苦苣苔科(Gesneriaceae)石蝴蝶屬(PetrocosmeaOliv.)中華石蝴蝶組(Sect.Petrocosmea)的多年生草本植物，其淡紫色的花冠和翠綠覆絨毛的葉片具有很高的觀賞價值。目前該物種僅在陜西省漢中市勉縣和略陽縣境內(nèi)的陰濕山溝發(fā)現(xiàn)，分布區(qū)域狹窄，數(shù)量稀少，且極易受環(huán)境和人為活動的影響，被列為秦嶺地區(qū)特有的國家二級珍稀瀕危保護(hù)植物[10-12]。本課題組通過無性繁殖技術(shù)人工繁殖秦嶺石蝴蝶10 000余株，在人工繁育的這些秦嶺石蝴蝶群體中發(fā)現(xiàn)了大量的花器官變異，除中國植物志中描述的正常秦嶺石蝴蝶花冠為二唇形，兩側(cè)對稱，上唇2裂、下唇3裂形態(tài)外，共發(fā)現(xiàn)了17種變異的花冠類型[13]。

本研究采用Illumina HiSeq 2500高通量測序技術(shù)對秦嶺石蝴蝶正常花朵2-3型(上唇2裂、下唇3裂，占65.43%)和在所有的變異類型中比例最多的變異類型2-4型(上唇2裂、下唇4裂，占26.15%)的開花早期和晚期進(jìn)行了轉(zhuǎn)錄組測序，為解析秦嶺石蝴蝶的花器變異調(diào)控機制提供一定的參考。

1 材料和方法

1.1 植物材料采集及形態(tài)學(xué)觀察

實驗材料于2019年7月2日采自略陽縣苗圃的人工仿野生栽培大棚中，所取秦嶺石蝴蝶均來自同一苗床，株齡、植株大小、生長狀態(tài)均相對一致。選取處于發(fā)育階段正?；ǘ漕愋?2-3型，N)和典型變異花朵類型(2-4型，V)的花苞(A)和成熟開放的花朵(B)為實驗材料。實驗材料分為4組：正常花苞(NA)、變異花苞(VA)、正?；ǘ?NB)和變異花朵(VB)，每組3個重復(fù)分別編號。保存在干冰中帶回實驗室，-80 ℃保存，選取完整的正?；ê妥儺惢ɑò突ǘ?，拍照保存。

1.2 方 法

1.2.1 秦嶺石蝴蝶花的總RNA提取、文庫構(gòu)建及轉(zhuǎn)錄組測序秦嶺石蝴蝶花總RNA提取方法參照Trizol Reagent方法進(jìn)行，用RNA專用瓊脂糖電泳檢測RNA的濃度和純度。通過Oligo(dT)磁珠富集總RNA中帶有polyA結(jié)構(gòu)的mRNA，采用離子打斷的方式，將mRNA打斷到300 bp。以片段mRNA為模板，用6堿基隨機引物和逆轉(zhuǎn)錄酶合成cDNA第一鏈，并以第一鏈cDNA為模板合成第二鏈cDNA，文庫構(gòu)建完成后，采用PCR擴增進(jìn)行文庫片段富集，之后根據(jù)片段大小進(jìn)行文庫選擇，文庫大小在450 bp。通過Agilent 2100 Bioanalyzer對文庫進(jìn)行質(zhì)檢，再對文庫總濃度及文庫有效濃度進(jìn)行檢測，形成單鏈文庫。建庫后，由上海派森諾公司完成測序，測序平臺為Illumina HiSeq。

1.2.2 數(shù)據(jù)處理、轉(zhuǎn)錄本拼接及功能注釋首先對原始下機的FASTQ文件數(shù)據(jù)進(jìn)行過濾，將帶接頭、長度小于50 bp、序列平均質(zhì)量在Q20以下的Reads進(jìn)行去除，使用Trinity軟件對高質(zhì)量序列進(jìn)行拼接得到轉(zhuǎn)錄本序列(Transcript)，挑選每個基因下最長的轉(zhuǎn)錄本作為該基因的代表序列Unigene。對Unigene與六大數(shù)據(jù)庫進(jìn)行基因功能注釋，基因功能注釋所用到的數(shù)據(jù)庫包括NR(NCBI non-redundant protein sequences)、GO(Gene Ontology)、KEGG(Kyoto Encyclopedia of Genes and Genome)、eggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)、Swiss-Prot和Pfam。

1.2.3 表達(dá)分析使用R語言的DESeq軟件包，根據(jù)表達(dá)量對各樣品進(jìn)行PCA主成分分析。采用DESeq對基因表達(dá)進(jìn)行差異分析，篩選差異表達(dá)基因條件為：表達(dá)差異倍數(shù)|log2FoldChange|>1，顯著性P<0.05。對不同分組之間的差異表達(dá)基因進(jìn)行統(tǒng)計，根據(jù)差異分析的結(jié)果，統(tǒng)計各比較組之間的共有和特有差異基因數(shù)量。

1.2.4 差異基因富集分析使用topGO進(jìn)行GO富集分析，分析時利用GO term注釋的差異基因?qū)γ總€term的基因列表和基因數(shù)目進(jìn)行計算，然后通過超幾何分布方法計算P(顯著富集的標(biāo)準(zhǔn)為P<0.05)，找出與整個基因組背景相比，差異基因顯著富集的GO term，從而確定差異基因行使的主要生物學(xué)功能。Pathway注釋主要采用KEGG的KAAS(KEGG Automatic Annotation Server)在線自動化注釋系統(tǒng)完成(http://www.genome.jp/tools/kaas/)。

1.2.5 差異表達(dá)基因的實時熒光定量PCR分析為了驗證基因差異表達(dá)的測序結(jié)果，選取差異表達(dá)的基因進(jìn)行熒光定量PCR驗證。用SPSS19.0軟件進(jìn)行生物統(tǒng)計學(xué)分析，通過雙側(cè)t檢驗方法進(jìn)行分析，P<0.05為差異有統(tǒng)計學(xué)意義。

2 結(jié)果與分析

2.1 兩個時期正常(2-3型)和典型變異(2-4型)花型分析

在解剖秦嶺石蝴蝶正常小花時發(fā)現(xiàn)，花器官形成早期上唇和下唇出現(xiàn)明顯的深裂，形成了早期的上唇2裂，下唇3裂(圖1，A)，成熟花朵花冠二唇形，上唇2裂，下唇3裂，命名為2-3型(圖1，B)；變異小花早期出現(xiàn)了上唇2裂，下唇4裂的深裂(圖1，C)，最終形成了變異類型中的上唇2裂，下唇4裂，簡稱2-4型(圖1，D)。

2.2 秦嶺石蝴蝶的花瓣數(shù)量變異的轉(zhuǎn)錄組分析

2.2.1 樣本基礎(chǔ)數(shù)據(jù)分析對12個樣本上機測序得到的原始數(shù)據(jù)分別進(jìn)行篩查和過濾，去除帶接頭、低質(zhì)量的序列，獲得高質(zhì)量序列34.55～47.27 Mb個，并統(tǒng)計其占原始序列的比例，占比為88.44%～93.42%，高質(zhì)量序列堿基數(shù)為5.18～7.09 Gb，GC含量為42.78%～44.56%(表1)。對高質(zhì)量序列進(jìn)行拼接共獲得Unigene序列113 889個，其總長度為99 295 924 bp，平均長度872 bp，N50長度為1 592 bp，GC含量為38.57%(表2)。

A.正常2-3型花苞(NA)；B.正常2-3型花朵(NB)；C.變異2-4型花苞(VA)；D.變異2-4型花朵(VB); a.上唇瓣；b.下唇瓣；c.花藥；d.子房；NA、NB、VA、VB下同圖1 秦嶺石蝴蝶正常2-3型和變異2-4型花苞和花朵A. Normal 2-3 type flower bud(NA); B. Normal 2-3 type flower(NB); C. Variation 2-4 type flower bud(VA); D. Variation 2-4 type flower(VB); a. Upper lip; b. Lower lip; c. Anther; d. Ovary； NA, NB, VA and VB are same as belowFig.1 Flower buds and flowers of normal type 2-3 and variation type 2-4 of Petrocosmea qinlingensis

2.2.2 秦嶺石蝴蝶花器官測序比對本研究將秦嶺石蝴蝶的Unigene與六大功能數(shù)據(jù)庫進(jìn)行比對，并根據(jù)基因的相似性進(jìn)行基因功能注釋，最終分別有52 677(NR: 46.25%)、46 750(eggNOG: 41.05%)、35 069(Swissprot: 30.79%)、21 088(GO: 18.52%)、19 663(KEGG: 17.27%)和18 215(Pfam: 15.99%)個Unigene獲得功能注釋，其中在所有數(shù)據(jù)庫中都被注釋到的Unigene數(shù)目有5 367個。通過與NR庫的同源序列相比，E值較低(1×10-30≤E<1×10-5)的Unigenes占Unigenes總數(shù)的42.81%，E值極低(1×10-100≤E<1×10-30)的Unigenes占Unigenes總數(shù)的35.19%，表明注釋的Unigenes序列與NR數(shù)據(jù)庫的同源序列具有高度相似性。

與NR庫進(jìn)行比對注釋，可以獲取本物種基因序列與近緣物種基因序列的相似性信息。在NR數(shù)據(jù)庫中，秦嶺石蝴蝶與旋蒴苣苔(Dorcocerashygrometricum) Unigene序列匹配度最高，占注釋Unigene總數(shù)的54.29%，其次是芝麻(Sesamumindicum)，占比10.38%，最低的是甜橙(Citrussinensis)，占比0.89%，還有22.21% Unigene與其他物種序列無匹配(圖2)。

2.2.3 主成分分析與基因差異表達(dá)分析對12個樣本進(jìn)行PCA主成分分析(圖3)，NB與VB距離較近，且組內(nèi)分布較為集中，而NA與VA雖然也可分為一組，但組內(nèi)分布較為分散?？傮w而言，花苞與花朵這兩個不同發(fā)育時期基因表達(dá)差異明顯，但4組材料基本分開，可以進(jìn)行后續(xù)分析。采用DESeq對基因表達(dá)進(jìn)行差異分析，統(tǒng)計各比較組之間的差異基因數(shù)量。韋恩圖(圖4)展示了各比較組之間的共有和特有差異基因數(shù)量。其中4個比較組中所共同擁有的差異基因數(shù)為62個，NA vs NB、VA vs VB、NA vs VA與NB vs VB獨有的差異基因數(shù)量分別為3 079、2 914、676與1 382個。NA vs NB與VA vs VB差異表達(dá)基因分別為10 665(5 248個上調(diào)表達(dá)基因和5 417個下調(diào)表達(dá)基因)和10 608個(5 075個上調(diào)表達(dá)基因和5 533個下調(diào)表達(dá)基因)，明顯高于NA vs VA與NB vs VB的2 917個(1 829個上調(diào)表達(dá)基因和1 088個下調(diào)表達(dá)基因)與3 500個(1 721個上調(diào)表達(dá)基因和1 779個下調(diào)表達(dá)基因)(表3)。

表1 樣本測序數(shù)據(jù)統(tǒng)計

表2 轉(zhuǎn)錄組測序的數(shù)據(jù)序列統(tǒng)計

圖2 NR數(shù)據(jù)庫E值分布圖和物種分布圖Fig.2 E-value distribution and species distribution of NR Database

圖3 12個樣本的PCA分析Fig.3 PCA analysis of 12 samples

2.3 秦嶺石蝴蝶花瓣數(shù)量變異的差異代謝通路分析

對差異表達(dá)基因進(jìn)行GO富集分析，按照分子功能、生物過程和細(xì)胞組分進(jìn)行GO分類，對每個比較組GO分類中p-value最小即富集最顯著的前20個GO條目，取-lg(p-value)進(jìn)行展示(表4)。其中NA vs VA差異基因在每個GO分類中富集最顯著的GO條目分別為細(xì)胞膜、跨膜轉(zhuǎn)運和反向轉(zhuǎn)運蛋白活性，NB vs VB差異基因在每個GO分類中富集最顯著的GO條目分別為酶抑制劑活性、細(xì)胞膜和催化活性的負(fù)調(diào)控，NA vs NB和VA vs VB差異基因在每個GO分類中富集最顯著的GO條目均為細(xì)胞膜、細(xì)胞壁組織或生物合成和催化活性。

圖4 組間差異表達(dá)基因韋恩圖Fig.4 Venn diagram of differentially expressed genes between groups

對差異表達(dá)基因進(jìn)行KEGG通路富集分析，將排名前20的KEGG pathway進(jìn)行展示(圖5)。通過富集因子(Rich factor)、FDR值和富集到此pathway上的基因個數(shù)來衡量富集的程度。Rich factor越大，表示富集的程度越大。FDR一般取值范圍為0～1，越接近于零，表示富集越顯著。NA vs VA差異基因主要富集在植物激素信號轉(zhuǎn)導(dǎo)和脂肪酸延長等途徑中；NB vs VB差異基因主要富集在戊糖和葡萄糖醛酸的相互轉(zhuǎn)化、苯丙烷類生物合成、硫代謝、類黃酮生物合成、玉米素的生物合成等途徑中；NA vs NB差異基因富集最顯著的KEGG pathway主要有植物激素信號轉(zhuǎn)導(dǎo)、植物-病原菌互作、氧化磷酸化、苯丙烷類生物合成、戊糖和葡萄糖醛酸的相互轉(zhuǎn)化及淀粉和蔗糖代謝等；VA vs VB差異基因富集最顯著的KEGG pathway主要有植物激素信號轉(zhuǎn)導(dǎo)、苯丙烷類生物合成、戊糖和葡萄糖醛酸酯的相互轉(zhuǎn)化及氧化磷酸化等。

表3 表達(dá)差異分析結(jié)果統(tǒng)計表

表4 差異表達(dá)基因GO富集分析

續(xù)表4 Continued Table 4

2.4 花瓣數(shù)量變異差異基因的表達(dá)分析

由于秦嶺石蝴蝶目前無參考基因組，因此注釋得到的許多基因功能未知或未預(yù)測，明確功能的基因很少。以校正后的P值為參數(shù)，將4個比較組(NA vs NB、VA vs VB、NA vs VA與NB vs VB)中的表達(dá)差異基因進(jìn)行交叉比較，結(jié)合功能注釋與富集結(jié)果，選出部分顯著性強的表達(dá)差異基因。進(jìn)一步對4個比較組進(jìn)行交叉比對分析，獲得6個可能與秦嶺石蝴蝶花器官發(fā)育相關(guān)且功能已知的基因PqMIF2、PqMYB340、PqMYB305、PqGATA12、PqCCD4和PqZBED。

將篩選出的6個基因利用qRT-PCR對轉(zhuǎn)錄組測序結(jié)果進(jìn)行驗證(圖6)，轉(zhuǎn)錄組結(jié)果中上調(diào)表達(dá)的基因其驗證結(jié)果也為上調(diào)表達(dá)，轉(zhuǎn)錄組結(jié)果中下調(diào)表達(dá)基因其驗證結(jié)果也為下調(diào)表達(dá)，表明基于轉(zhuǎn)錄組測序數(shù)據(jù)的基因差異表達(dá)分析結(jié)果是可信的。6個基因qRT-PCR分析得到的相對表達(dá)量中PqMIF2、PqMYB340、PqMYB305與PqCCD4在開放的花朵中表達(dá)量極為豐富，但在花苞組織中表達(dá)量極低；而PqGATA12則明顯在花苞中表達(dá)量高于開放花朵中。此外，PqZBED在正?；ㄆ鞴俨牧现械谋磉_(dá)量高于變異花器官。這與轉(zhuǎn)錄組表達(dá)譜分析趨勢一致，說明基因表達(dá)譜的分析結(jié)果可靠。

圖6 篩選出的6個基因在4個樣品中的表達(dá)Fig.6 Expression of six genes screened in four samples

3 討 論

花器官數(shù)目的變化可通過原基分裂、原基融合或次生原基的形成而產(chǎn)生，花原基細(xì)胞數(shù)目的增減可能導(dǎo)致花器官數(shù)目增減，而引起花器官變異卻是基因和環(huán)境條件共同作用的結(jié)果[14]。在獲得許可的情況下于2019年8月中旬對略陽縣境內(nèi)的野生秦嶺石蝴蝶進(jìn)行了野外觀察和記錄。遺憾的是在所有觀察到的花朵中未發(fā)現(xiàn)變異現(xiàn)象，所有花朵與中國植物志中的記錄相一致，而將采集到的野生秦嶺石蝴蝶植株帶回人工控制的環(huán)境和營養(yǎng)土進(jìn)行栽培，發(fā)現(xiàn)部分植株在移栽后開花時開始表現(xiàn)出花瓣數(shù)量增多的現(xiàn)象。這表明，豐富的營養(yǎng)物質(zhì)可能是導(dǎo)致秦嶺石蝴蝶花器變異的主要原因。為了探明秦嶺石蝴蝶花瓣數(shù)量變異原因，本研究采用Illumina HiSeq 2500高通量測序技術(shù)對秦嶺石蝴蝶兩種花型發(fā)育的早期和晚期進(jìn)行轉(zhuǎn)錄組測序，挖掘參與其花發(fā)育相關(guān)的差異基因并探討花器變異的可能機制。

高通量測序結(jié)果顯示，秦嶺石蝴蝶基因注釋與旋蒴苣苔植物的基因數(shù)據(jù)庫匹配度最高，達(dá)到了54.29%。秦嶺石蝴蝶屬于石蝴蝶屬植物，旋蒴苣苔屬于旋蒴苣苔屬，二者同屬于苦苣苔科，表明秦嶺石蝴蝶的大部分基因測序結(jié)果可靠，可以用于后續(xù)基因功能分析。PCA和組間差異表達(dá)基因分析結(jié)果顯示，秦嶺石蝴蝶花器官發(fā)育早期和晚期的差異表達(dá)基因數(shù)量遠(yuǎn)遠(yuǎn)大于同一時期的正?；ㄆ鞴俸妥儺惢ㄆ鞴俚牟町惐磉_(dá)基因數(shù)量。這表明，相對于整個花器官發(fā)育的過程中差異表達(dá)的基因，影響花器官向著某一方向變異的基因數(shù)目并不是很大。

GO富集分析結(jié)果顯示，早期的花器官變異可能與蛋白質(zhì)的跨膜轉(zhuǎn)運、定位和細(xì)胞分裂素的合成等生物過程有關(guān)，例如轉(zhuǎn)運蛋白、半乳糖-棉子半乳糖轉(zhuǎn)移酶、Rho鳥嘌呤核苷酸交換因子等具有很高的活性，而后期花器官的變異可能與果膠、糖類、蛋白質(zhì)的水解有關(guān)，例如內(nèi)肽酶、肽酶抑制劑活性、核苷酸酶等有關(guān)，同時在這一過程中水楊酸的調(diào)控也起到了重要作用。無論是正常還是變異的花器官，他們在不同發(fā)育的時期都與細(xì)胞壁的形成和組裝、糖代謝和氧化還原生物過程有關(guān)。而KEGG富集分析也進(jìn)一步證實了這些結(jié)果。在花器官發(fā)育的早期變異過程中，植物激素信號相關(guān)的基因差異表達(dá)非常顯著，而在后期的花器官變異過程中則表現(xiàn)出戊糖和葡萄糖類和苯丙烷類的代謝調(diào)節(jié)相關(guān)基因發(fā)生了差異表達(dá)。

值得關(guān)注的是植物激素信號轉(zhuǎn)導(dǎo)通路差異基因富集最顯著，推測這些差異基因可能參與多種激素代謝調(diào)控植物器官發(fā)育。深入研究關(guān)鍵基因所參與這些代謝的反應(yīng)，對揭示秦嶺石蝴蝶花器官變異的分子機制具有重要意義。MIF在花期結(jié)束過程中抑制細(xì)胞增殖，在花的發(fā)育過程中控制心皮的數(shù)量，MIF的表達(dá)不足可阻止配子體組織的形成[15]。此外，研究發(fā)現(xiàn)，擬南芥的Mini Zinc Finger (AtMIF2)蛋白和其在番茄(Solanumlycopersicum)中的同源蛋白SlIMA作為結(jié)合蛋白參與了花分生組織的終止，在控制心皮數(shù)量的發(fā)育階段起到重要作用[16]。在本研究中MIF在秦嶺石蝴蝶的花器官變異中表達(dá)比較低，但在花器官發(fā)育的調(diào)控過程中卻表達(dá)很高，表明在引起花器官變異過程中此基因并未發(fā)揮重要的作用，但在整個花器官發(fā)育過程中發(fā)揮著重要作用。同樣，轉(zhuǎn)錄因子MYB的變化趨勢與MIF基因一致。GATA轉(zhuǎn)錄調(diào)節(jié)因子廣泛存在于真核生物中。在調(diào)節(jié)植物生長、開花時間、花的發(fā)育以及光周期等過程中發(fā)揮著重要作用[17]。本研究中發(fā)現(xiàn)GATA基因無論在花器官早期還是后期的變異過程中均有非常高的表達(dá)量，可能是引起秦嶺石蝴蝶花器官變異的一個重要基因，在后期的基因功能驗證中值得關(guān)注。大量研究表明，類胡蘿卜素裂解雙加氧酶4(CCD4)在調(diào)節(jié)植物花色和果實顏色中具有重要作用[18]。在秦嶺石蝴蝶轉(zhuǎn)錄組分析結(jié)果中，CCD4在開放后的花朵中過量表達(dá)，而在花苞中表達(dá)量較低，這與兩者的顏色差異相一致。

對于秦嶺石蝴蝶而言，這種人工條件下的花器官變異是有利的，增加的花器官數(shù)目可以增加其觀賞性，在今后觀賞價值的開發(fā)方面具有積極意義。本研究結(jié)果表明，環(huán)境改變可能在花器官發(fā)育的早期引起了某些花器官發(fā)育調(diào)控基因的改變，從而導(dǎo)致花器官發(fā)育和表達(dá)調(diào)控異常，但其分子機制仍需進(jìn)一步研究、分析和驗證。