轉(zhuǎn)錄組和表達(dá)譜分析揭示PEBP基因家族成員參與調(diào)節(jié)滴水珠的珠芽發(fā)育

摘 要：為進(jìn)一步揭示半夏屬植物珠芽發(fā)育的分子機(jī)理，本文利用Illumina Hiseq平臺對滴水珠進(jìn)行轉(zhuǎn)錄組和表達(dá)譜測序。轉(zhuǎn)錄組測序獲得6.68 Gb的數(shù)據(jù)，組裝去冗余后得到66，907個(gè)Unigene，平均長度為893 bp。將Unigene比對到Nr、KOG、GO和KEGG數(shù)據(jù)庫中，分別獲得了34，947、27，886、9，149和27，006個(gè)注釋信息。葉柄、珠芽和塊莖表達(dá)譜測序分別獲得24.08、24.02和24.06 Mb的數(shù)據(jù)庫，其中6，961個(gè)Unigene的表達(dá)量在葉柄、珠芽與塊莖之間同時(shí)存在差異（6，576和4，021個(gè)差異Unigene分別比對到GO和KEGG數(shù)據(jù)庫）。轉(zhuǎn)錄組測序獲得13個(gè)Unigene為PEBP基因家族成員相關(guān)序列，其中4個(gè)為差異表達(dá)。結(jié)果表明轉(zhuǎn)錄組和表達(dá)譜測序可用于篩選調(diào)控珠芽發(fā)育的候選基因，PEBP基因家族成員在珠芽發(fā)育調(diào)控中的作用值得進(jìn)一步研究。

關(guān)鍵詞：滴水珠；珠芽發(fā)育；轉(zhuǎn)錄組；表達(dá)譜；PEBP基因家族

中圖分類號：Q943

文獻(xiàn)標(biāo)識碼： A

滴水珠（Pinellia cordata N. E. Brown）是主要分布于我國安徽、福建、廣東、廣西、貴州、湖北、湖南等地區(qū)的天南星科植物[1]，具有解毒消腫、散痰止痛等功效[2]。滴水珠的自然繁殖主要通過在葉柄下部和上部（葉片基部）形成珠芽進(jìn)行。滴水珠及同屬植物半夏的珠芽發(fā)育形態(tài)、解剖學(xué)已經(jīng)報(bào)道[3-4]，其發(fā)育特征除啟動(dòng)位置和后期成熟過程外與龍舌蘭、臺閩苣苔、山藥和淡黃花百合等類似[5-7]。對于植物珠芽發(fā)育的分子調(diào)控機(jī)制研究而言，Wang 等[8]在臺閩苣苔上發(fā)現(xiàn)GFLO基因與珠芽形成有關(guān)，珠芽形成過程中該基因的表達(dá)量呈現(xiàn)出下調(diào)趨勢；Abraham-Juarez等[9]研究表明在龍舌蘭中KNOX I基因表達(dá)量隨著珠芽成熟而逐漸增加；Sandoval等[10]發(fā)現(xiàn)在龍舌蘭珠芽形成過程中MADS1， 2， 4， 6， 7基因均呈現(xiàn)出表達(dá)量減少的相同模式。目前在珠芽發(fā)育的分子基礎(chǔ)上研究報(bào)道少且不同植物珠芽發(fā)育的共同特性尚不明確。此外，磷脂酰乙醇胺結(jié)合蛋白基因（PEBP）家族包含F(xiàn)T、MFT和TFL三個(gè)亞家族，其不同家族成員可參與調(diào)節(jié)植物開花[11]、馬鈴薯塊莖形成[12]和洋蔥鱗莖形成[13]等有性和無性生殖發(fā)育過程。但是，目前尚未見報(bào)道PEBP基因家族成員參與調(diào)節(jié)珠芽發(fā)育。

轉(zhuǎn)錄組和表達(dá)譜測序可用于揭示植物器官發(fā)育過程中的調(diào)節(jié)基因，如程立寶等[14]發(fā)現(xiàn)了蓮藕根狀莖發(fā)生過程中的86個(gè)相關(guān)基因。轉(zhuǎn)錄組和表達(dá)譜測序也用于植物珠芽發(fā)育研究：彭斌[15]在山藥的轉(zhuǎn)錄組測序中篩選出23個(gè)可能與珠芽發(fā)生有關(guān)的候選基因；姜福星等[16]在萬年青的轉(zhuǎn)錄組測序中發(fā)現(xiàn)EMB基因、TATA結(jié)合蛋白基因和FKBP蛋白等基因在珠芽的形成過程中表達(dá)量發(fā)生了明顯變化；葉德等[17]在半夏珠芽和葉柄的轉(zhuǎn)錄組測序中篩選出了19個(gè)與激素有關(guān)的差異基因，其中玉米素、茉莉酸、脫落酸4種激素合成代謝相關(guān)的基因在珠芽中的表達(dá)量高于莖中，赤霉素相反。但是，不同植物珠芽發(fā)育的轉(zhuǎn)錄組特征異同的揭示需要更多的轉(zhuǎn)錄組數(shù)據(jù)。

本研究以滴水珠為材料進(jìn)行轉(zhuǎn)錄組和表達(dá)譜測序，分析珠芽發(fā)育過程中的表達(dá)譜宏觀差異，并進(jìn)一步分析PEBP家族成員的表達(dá)差異，為進(jìn)一步研究半夏屬植物的珠芽發(fā)育機(jī)制提供基礎(chǔ)。

1 材料與方法

1.1 實(shí)驗(yàn)材料

滴水珠采集自浙江溫州，盆栽于實(shí)驗(yàn)室備用。在植物的珠芽發(fā)育初期采集植株，純水沖洗后在冰上進(jìn)行取樣：葉柄（pc1）、珠芽（≤3 mm， pc2）和塊莖（pc3）。取得鮮樣用液氮速凍后，放入-80 ℃冰箱保存。

1.2 總RNA的提取和文庫的構(gòu)建

CTAB法提取總RNA，DNase I處理去除DNA，然后用磁珠法富集mRNA。富集后mRNA等量分為兩份：一份將三個(gè)組織鮮樣的mRNA等量混合后進(jìn)行轉(zhuǎn)錄組測序的相關(guān)分析；另一份將各組織mRNA進(jìn)行單獨(dú)表達(dá)譜分析。向獲得的mRNA中加入打斷試劑使得mRNA片段化，以mRNA片段為模板進(jìn)行反轉(zhuǎn)錄生成cDNA。之后對cDNA的末端進(jìn)行修復(fù)、并在3'末端加上“A”堿基并連接接頭。之后對連接產(chǎn)物進(jìn)行擴(kuò)增并將PCR產(chǎn)物環(huán)化，最終文庫構(gòu)建成功。構(gòu)建好的文庫利用Agilent 2100 Bioanalyzer進(jìn)行質(zhì)檢，合格后使用Illumina Hi Seq平臺進(jìn)行測序。

1.3 測序數(shù)據(jù)的加工、組裝和注釋

等量混合mRNA建庫后測序得到的raw reads進(jìn)行過濾處理獲得Clean reads，再利用Trinity軟件[18]對Clean reads進(jìn)行de novo組裝獲得轉(zhuǎn)錄組序列。之后使用Tgicl軟件[19]對轉(zhuǎn)錄本進(jìn)行處理后獲得Unigene。最后，將所有組裝的Unigene利用BLAST[20]比對到non redundant protein database （Nr）、the Kyoto Encyclopedia of Genes and Genomes （KEGG）、clusters of eukaryotic Orthologous Groups（KOG）蛋白質(zhì)數(shù)據(jù)庫（E<10-5）；利用Blast2GO[21]獲得GO注釋；利用InterProScan5[22]將Unigene注釋到InterPro數(shù)據(jù)庫。

1.4 差異基因的GO注釋和KEGG注釋

單獨(dú)建庫mRNA（1.2）表達(dá)譜測序后以de nevo組裝、注釋的轉(zhuǎn)錄組序列為參考并利用FPKM法（Fragments Per kb per Million fragments）計(jì)算各個(gè)樣本Unigene表達(dá)量[23]；利用FDR法[24]來進(jìn)行差異表達(dá)的真實(shí)性（FDR≤0.001且倍數(shù)差異在2倍以上的基因被定義為差異基因）。利用Blast2GO[21]軟件將差異基因進(jìn)行GO功能注釋，從而獲得差異表達(dá)基因的GO注釋和分類。利用BLAST[20]軟件將差異基因比對到KEGG數(shù)據(jù)庫中，進(jìn)一步獲得差異表達(dá)基因的Pathway的注釋信息。

2 結(jié)果

2.1 轉(zhuǎn)錄組序列的de novo組裝

去除質(zhì)量較低、接頭污染以及N含量過高的序列后，轉(zhuǎn)錄組測序共獲得44.54 Mb高質(zhì)量的Clean reads，Q20和Q30的百分比分別為98.86%和96.53%。Clean reads組裝后獲得180，862個(gè)Contigs，Contigs的序列長度范圍是≥200 bp，序列的平均長度和N50的長度分別為547 bp和940 bp，GC%為48.73%。Contigs的片段大小分布：200-500 nt： 73.96%； 600-1000 nt： 12.36%； 1100-1500 nt： 5.61%； 1600-2000 nt： 3.46%； 2100-2500 nt： 1.96%； 2600-3000 nt： 1.04%； ≥3000： 1.60%。轉(zhuǎn)錄本進(jìn)行聚類去冗余得到了66，907個(gè)Unigene，Unigene的平均長度和N50長度分別為893 bp和1，540 bp。Unigene的長度分布：300-500 nt： 50%； 600-1000 nt： 19.82%； 1100-1500 nt： 11.55%； 1600-2000 nt： 7.76%； 2100-2500 nt： 4.55%； 2600-3000 nt： 2.51%；≥3000： 3.80%。

2.2 轉(zhuǎn)錄組序列的KOG功能注釋

利用BLAST將Unigene比對到Nr、Nt、Swiss-prot、KEGG、KOG、Interpro和GO 七個(gè)公共數(shù)據(jù)庫中（E<10-5）。結(jié)果顯示，有34，947（52.23%）個(gè)Unigene比對到Nr數(shù)據(jù)庫中，22，133（33.08%）個(gè)Unigene比對到Nt數(shù)據(jù)庫，23，915（35.74%）個(gè)Unigene比對到Swiss-Prot數(shù)據(jù)庫，比對到KEGG、KOG、Interpro 和 GO數(shù)據(jù)庫的Unigene數(shù)分別為：27，006（40.36%）、 27，886（41.68%）、 25，962（40.30%）和 9，149（13.67%）。

比對到KOG數(shù)據(jù)庫中的共有27，886個(gè)Unigene根據(jù)功能劃分屬于25個(gè)功能類型（圖1）。其中，“一般功能預(yù)測”的Unigene數(shù)最多，有8，669個(gè)Unigene，其次是“信號轉(zhuǎn)導(dǎo)機(jī)制”，有5，022個(gè)Unigene，最少的是“細(xì)胞運(yùn)動(dòng)”，僅有43個(gè)。

2.3 轉(zhuǎn)錄組序列的GO功能注釋和表達(dá)譜差異

比對到GO數(shù)據(jù)庫的Unigene參與生物過程（Biological process）、細(xì)胞組成（Cellular component）和分子功能（Molecular function）3大類別55個(gè)小類。生物過程中：“代謝過程”的Unigene數(shù)最多（4，628個(gè)）， “細(xì)胞進(jìn)程”次之（4，445個(gè)），最少的是“細(xì)胞凋亡”和“移動(dòng)”（分別僅有1個(gè)）。細(xì)胞組分中，“細(xì)胞”所含的Unigene數(shù)最多，有3，823個(gè)，其次是“細(xì)胞部分”，有3，787個(gè)Unigene，“細(xì)胞核”是Unigene數(shù)最少的類別，僅有11個(gè)。分子功能中，“催化活化”的Unigene數(shù)最多（4，389個(gè)）。其次是“結(jié)合”（3，985個(gè)），“蛋白標(biāo)簽”僅有1個(gè)Unigene，是分子功能中Unigene數(shù)最少的類別。

表達(dá)譜測序一共測了3個(gè)樣品，葉柄、珠芽和塊莖各獲得24.08、24.02和24.06 Mb的數(shù)據(jù)，過濾后的reads數(shù)目占總raw reads的比例分別為99.78%、99.53%和99.70%。差異篩選獲得6，961個(gè)表達(dá)量在葉柄與珠芽、珠芽與塊莖之間同時(shí)存在差異的Unigene，其中有6，576個(gè)比對到了GO數(shù)據(jù)庫：細(xì)胞組分2，680個(gè)、生物學(xué)過程2，556個(gè)、分子功能1，340個(gè)。而轉(zhuǎn)錄組序列GO數(shù)據(jù)庫比對時(shí)分別有14，150、18，825和9，724個(gè)Unigene比對到生物學(xué)過程、細(xì)胞組分和分子功能三大類。轉(zhuǎn)錄組與表達(dá)譜差異基因GO比對結(jié)果都表現(xiàn)出“從細(xì)胞組分到生物學(xué)過程再到分子功能Unigene數(shù)量逐漸減少”的趨勢。轉(zhuǎn)錄組注釋中“代謝過程”（4，628個(gè)）所含的Unigene數(shù)是最多的，其次是“細(xì)胞進(jìn)程”（4，445個(gè)），“細(xì)胞”（3，823個(gè)）和“細(xì)胞部分”（3，787個(gè)）。表達(dá)譜差異基因功能富集結(jié)果發(fā)現(xiàn)，所含Unigene數(shù)最多的是“催化活性”（666個(gè)），其次是“代謝過程”（646個(gè)），“細(xì)胞進(jìn)程”（611個(gè)）和“細(xì)胞”（528個(gè)）。GO功能注釋的單獨(dú)類別中，轉(zhuǎn)錄組注釋和表達(dá)譜差異基因注釋的Unigene數(shù)存在不同（圖2）。

2.4 轉(zhuǎn)錄組序列的KEGG功能注釋和表達(dá)譜差異

從轉(zhuǎn)錄組的27，006個(gè)Unigene中發(fā)現(xiàn)了137條KEGG代謝通路。其中參與“全局路徑”的Unigene最多（6，426個(gè)），其次是“翻譯”途徑（2，580個(gè)）、“糖代謝”（2，287個(gè)）、“折疊、分類和降解”（1，705個(gè)）和“信號轉(zhuǎn)導(dǎo)”（1，423個(gè)），Unigene數(shù)最少的是“抗藥性：抗菌劑”代謝途徑（僅8個(gè) Unigene）。表達(dá)譜差異篩選獲得的6，961個(gè)在葉柄、珠芽、塊莖之間同時(shí)存在差異的Unigene中有4，021個(gè)差異基因比對到了KEGG數(shù)據(jù)庫中。其中參與“全局路徑”的最多（1，044個(gè)），其次是“碳水化合物代謝”（369個(gè)）、“翻譯”（353個(gè)）、“其他次生代謝”（216個(gè)）和“信號轉(zhuǎn)導(dǎo)”（206個(gè)），最少的是“抗藥性：抗菌劑”僅有1個(gè)Unigene。比對到KEGG數(shù)據(jù)庫各代謝途徑的轉(zhuǎn)錄組注釋和表達(dá)譜篩選Unigene數(shù)量和比例存在差異（圖3）。

2.5 PEBP家族差異基因篩選

轉(zhuǎn)錄組序列中含有13個(gè)Unigene屬于PEBP基因家族成員，其中有10個(gè)Unigene屬于FT亞家族，TFL亞家族有1個(gè)Unigene，MFT亞家族有2個(gè)Unigene。表達(dá)譜差異基因篩選結(jié)果表明4個(gè)PEBP基因家族的Unigene在不同組織樣品中具有表達(dá)差異：其中屬于FT亞家族的有2個(gè)，屬于的MFT亞家族有2個(gè)，而沒有屬于TFL亞家族的成員（表1）。

3 討論

滴水珠轉(zhuǎn)錄組測序獲得了180，862個(gè)Contig和66，907個(gè)Unigene。其中有27，886（41.68%）個(gè)Unigene成功地比對到KOG數(shù)據(jù)庫，GO數(shù)據(jù)庫和KEGG數(shù)據(jù)庫分別有9，149（13.67%）和27，006（40.36%）個(gè)Unigene比對上。葉德等[17]在半夏轉(zhuǎn)錄組測序時(shí)獲得90，175個(gè)Unigene，KOG、GO和KEGG比對的Unigene數(shù)分別為46，237（51.27%）、31，622（35.07%）和22，961（25.46%）。與半夏轉(zhuǎn)錄組相比，滴水珠比對到KOG和GO數(shù)據(jù)庫的Unigene數(shù)比半夏的多，但比對到KEGG數(shù)據(jù)庫結(jié)果則相反，推測可能是因?yàn)槲锓N特異性和測序量大小不同造成的。

滴水珠表達(dá)譜測序獲得6，961個(gè)表達(dá)量在葉柄、珠芽、塊莖之間同時(shí)存在差異表達(dá)的Unigene，數(shù)量低于轉(zhuǎn)錄測序差異篩選的半夏（15，484）和山藥（7，647）[15，17]。差異Unigene數(shù)的不同可能由物種差異及測序量大小、比較的取樣組織和差異篩選方法不同導(dǎo)致。滴水珠、半夏和山藥的測序量分別是66，907、90，175和70，480個(gè)Unigene；三種植物比較的取樣組織分別是葉柄vs珠芽vs塊莖、葉柄vs珠芽和鐵棍珠芽vs鐵棍葉腋vs花籽葉腋及下表皮；而表達(dá)量定量方法分別是表達(dá)譜、轉(zhuǎn)錄組和轉(zhuǎn)錄組。

PEBP基因家族（包含有FT、TFL和MFT三個(gè)亞家族）廣泛存在于動(dòng)物、植物和微生物中[25]。PEBP基因家族參與植物的許多發(fā)育過程，如植物開花[11]、控制氣孔的開放[26]、花序分生組織的穩(wěn)定[27]等。雖然滴水珠、半夏的珠芽在結(jié)構(gòu)上類似于塊莖[3]，淡黃花百合的珠芽在結(jié)構(gòu)上類似于鱗莖[7]，但是珠芽和塊莖/鱗莖發(fā)育的共同分子基礎(chǔ)目前尚無報(bào)道。滴水珠表達(dá)譜差異Unigene中具有4個(gè)PEBP家族同源序列（FT： 2； MFT： 2； TFL：0），結(jié)果顯示珠芽和塊莖/鱗莖發(fā)育具有共同分子基礎(chǔ)，PEBP基因家族在其中的作用值得進(jìn)一步研究。

參考文獻(xiàn)：

[1]Li H.， G.H. Zhu， P.C. Boyce， et al. Flora of China[C]. Vol 23， Beijing：Science press， 2010：3-79.

[2]宋立人， 洪恂. 丁緒亮，等.現(xiàn)代中藥學(xué)大辭典[M]. 北京： 人民衛(wèi)生出版社， 2001： 2354.

[3]羅睿， 杜禹珊， 孫瑩瑩， 等. 半夏珠芽發(fā)育過程的形態(tài)學(xué)和解剖學(xué)研究[J]. 西北植物學(xué)報(bào)， 2014， 34（9）： 1776-1781.

[4]朱燕燕， 羅睿， 陳海麗， 等. 滴水珠珠芽發(fā)育過程研究[J]. 廣西植物， 2018（2）：225-232.

[5]謝德镕， 曹祥志， 張培源. 薯蕷零余子發(fā)育與形態(tài)學(xué)本質(zhì)探討[J]. 西北農(nóng)業(yè)大學(xué)學(xué)報(bào)， 1993， 21（1）： 77-80.

[6]余炳生， 張儀. 薯蕷地下塊莖發(fā)生的解剖學(xué)研究[J]. 植物學(xué)報(bào)， 1994， 11： 52.

[7]李騰， 李少群， 羅睿. 淡黃花百合珠芽發(fā)育過程的形態(tài)學(xué)與解剖學(xué)研究[J]. 西北植物學(xué)報(bào)， 2012（1）： 85-89.

[8]Wang， C.N.，M. Mller， Q.C.B. Cronk. Altered expression of GFLO， the Gesneriaceae homologue of Floricaula/Leafy， is associated with the transition to bulbil formation in Titanotrichum oldhamii[J]. Dev. Genes Evol.， 2004， 214： 122-127.

[9]Abraham-Juarez， M.J.， A. Martinez-Hernandez， M.A. Leyva-Gonzalez， et al. ClassI KNOX genes are associated with organogenesis during bulbil formation in Agave tequilana[J]. J. Exp. Bot.， 2010， 61： 4055-4067.

[10]Sandoval， S.C.D.， M.J. Abraham Juárez， J. Simpson. Agave tequilana MADS genes show novel expression patterns in meristems， developing bulbils and floral organs[J]. Sex. Plant Reprod， 2012， 25： 11-26.

[11]Li， C.， Y.N. Zhang， K. Zhang， et al. Promoting owering， lateral shoot outgrowth， leaf development， and ower abscission in tobacco plants overexpressing cotton FLOWERINGLOCUS T （FT）-like gene Gh FT1[J]. Front. Plant Sci.， 2015， 6： 454.

[12]González-Schain， N.D.， M. Díaz-Mendoza， M. Zurczak， et al. Potato CONSTANS is involved in photoperiodic tuberizationin a graft-transmissible manner[J]. Plant J.， 2012， 70： 678-90.

[13]Lee， R.， S. Baldwin， F. Kenel， et al. FLOWERING LOCUS T genes control onion bulb formation and owering[J]. Nat Commun.， 2013， 4： 2884.

[14]程立寶， 齊曉花， 髙學(xué)雙， 等. 蓮藕根狀莖膨大相關(guān)基因的挖掘與表達(dá)分析[J]. 園藝學(xué)報(bào)， 2012， 39（3）： 834-846.

[15]彭斌. 山藥基源植物分子鑒別體系建立與珠芽相關(guān)基因的轉(zhuǎn)錄組分析[D]. 南京： 南京農(nóng)業(yè)大學(xué)， 2016.

[16]姜福星， 魏丕偉， 吳生， 等. 白花虎眼萬年青葉上珠芽的轉(zhuǎn)錄組分析[J]. 分子植物育種， 2017（2）： 519-531.

[17]葉德， 楊支力， 李東海， 等. 基于高通量測序半夏珠芽轉(zhuǎn)錄組研究[J]. 浙江理工大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2017（2）： 282-288.

[18]Grabherr， M.G.， B.J. Haas， M. Yassour， et al. Trinity： reconstructing a full-length transcriptome without a genome from RNA-Seq data[J]. Nature Biotechnology， 2011， 29： 644-652.

[19]Pertea G.， X.Q. Huang， F. Liang， et al.TIGR Gene Indices clustering tools（TGICL）： a software system for fast clustering of large EST dataste[J]. Bioinformatics， 2003， 19（5）： 651.

[20]Altschul S.F.， W. Gish， E.W. Myers， et al.. Basic local alignment search tool[J]. J Mol Biol，1990， 215（3）：403-410.

[21]Conesa， A.， S. Gotz， J.M. Garcia-Gomez， et al. Blast2GO： a universal tool for annotation， visualization and analysis in functional genomics research[J]. Bioinformatics， 2005， 21： 3674-3676.

[22]Quevillon， E.， V. Silventoinen， S. Pillai， et al. InterProScan： protein domains identifier[J]. Nucleic Acids Res.， 2005， 33： W116-W120.

[23]Mortazavi A， Williams B.A， McCue， K， et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq[J]. Nature methods， 2008， 5（7）：621-628.

[24]Benjamini， Y. ， D. Yekutieli. The control of the 1 discovery rate in multiple testing underdependency[J]. The Annals of Statistics， 2001， 29： 1165-1188.

[25]張禮鳳， 徐冉， 張彥威， 等. 大豆PEBP基因家族的初步分析[J]. 植物遺傳資源學(xué)報(bào)， 2015， 16（1）：151-157.

[26]Kinoshita， T.， N. Ono， Y. Hayashi， et al. FLOWERING LOCUS T Regulates Stomatal Opening[J]. Curr Biol.， 2011， 21（14）： 1232-1238.

[27]Liu， L.， S. Farrona， S. Klemme et al. Post-fertilization expression of FLOWERING LOCUS T suppresses reproductive reversion[J]. Front. Plant Sci.， 2014， 5：1-7.

（責(zé)任編輯：江 龍）