花生異型雄蕊轉(zhuǎn)錄組分析

王海霞，王銘倫，丁雨龍

(1. 南京林業(yè)大學(xué)： a． 南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心，b． 生物與環(huán)境學(xué)院，c． 竹類研究所，江蘇 南京 210037；2. 青島農(nóng)業(yè)大學(xué)，山東 青島 266109)

花生(ArachishypogaeaLinn.)為嚴(yán)格的閉花受精植物，每朵花內(nèi)有10枚雄蕊，其中，具4枚長橢球狀花藥的雄蕊為對(duì)萼長藥雄蕊(antesepal stamen，Sts)，藥室較大，花絲較長；具4枚圓球狀花藥的雄蕊為對(duì)瓣圓藥雄蕊(antepetal stamen，Stp)，藥室較小，花絲較短；另2枚雄蕊為退化雄蕊(staminode，Sta)，僅存花絲，未見花藥?；ㄉ哂械倪@類雄蕊為異型雄蕊[1]。異型雄蕊現(xiàn)象從18世紀(jì)開始就引起了科學(xué)家們的關(guān)注[2-3]，一些研究者針對(duì)異型雄蕊的種類分布[4-7]以及花粉活性、花粉組織化學(xué)成分和結(jié)實(shí)率等[8-11]方面開展了廣泛的研究，然而，異型雄蕊的遺傳基礎(chǔ)和分子調(diào)控機(jī)制目前仍未知。

花生的基因組較大、遺傳基礎(chǔ)狹窄，且其花器官的性狀穩(wěn)定、遺傳多態(tài)性極低，因而，花生的花器官突變體很少，限制了花生的花發(fā)育分子生物學(xué)研究。雄蕊發(fā)育在被子植物的花器官進(jìn)化中起關(guān)鍵作用，但是與其他花器官(如花瓣)相比，其分子調(diào)控機(jī)制尚處于探索階段，對(duì)花生異型雄蕊發(fā)育的分子機(jī)制進(jìn)行探討可為雄蕊分化的分子機(jī)制研究提供思路。

作者采用RNA-seq技術(shù)對(duì)花生異型雄蕊進(jìn)行轉(zhuǎn)錄組分析，挖掘和探究與花生雄蕊分化相關(guān)的候選基因和生物學(xué)途徑，以期解釋花生異型雄蕊的形成機(jī)制，為花生新品種選育提供理論基礎(chǔ)。

1 材料和方法

1.1 材料

供試花生品種為青島農(nóng)業(yè)大學(xué)花生研究中心王銘倫教授提供的品種‘Florunner’，該品種由品種‘Early Runner’與‘Florispan’雜交獲得，為普通型蔓生品種，花多且花型較大[12]。于2017年4月末將品種‘Florunner’的種子覆膜種植于青島農(nóng)業(yè)大學(xué)試驗(yàn)田中，并于同年6月中旬(盛花期)選取100株發(fā)育良好的植株，每天約17：00時(shí)于各植株同一位置采集發(fā)育成熟、待開放的花蕾(即花瓣已撐破萼片、微露黃色花瓣)，共采集約400枚花蕾。

將采集的花蕾置于冰盒內(nèi)，在解剖鏡下將雄蕊分為對(duì)萼長藥雄蕊、對(duì)瓣圓藥雄蕊和退化雄蕊3組，每組約800枚雄蕊(每組設(shè)置3個(gè)生物學(xué)重復(fù))；將雄蕊分別置于液氮中速凍，-80 ℃保存，送至上海伯豪生物技術(shù)有限公司進(jìn)行RNA提取及測(cè)序分析。

1.2 方法

1.2.1 總RNA的提取和檢測(cè) 采用RNAqueousTMphenol-free total RNA提取試劑盒(美國Ambion公司)分別提取3組雄蕊的總RNA，獲得的總RNA經(jīng)Agilent Bioanalyzer 2100生物分析儀(美國Agilent Technologies公司)電泳檢測(cè)合格，使用RNAClean XP Kit試劑盒(美國Beckman Coulter公司)和RNase-Free DNase Set試劑盒(德國Qiagen公司)進(jìn)行純化。

1.2.2 文庫制備和測(cè)序 對(duì)純化的總RNA進(jìn)行mRNA的富集和片段化，采用SuperScript Ⅱ Reverse Transcriptase試劑盒(美國Invitrogen公司)合成cDNA第1鏈，然后合成cDNA第2鏈，隨之進(jìn)行末端修復(fù)，在3′末端加腺苷酸并連接測(cè)序接頭，富集、純化后構(gòu)建cDNA文庫；采用Illumina HiSeqTM2000高通量測(cè)序平臺(tái)(美國Illumina公司)對(duì)cDNA文庫進(jìn)行雙端測(cè)序。

1.2.3 原始數(shù)據(jù)統(tǒng)計(jì)和功能注釋 原始測(cè)序文件經(jīng)堿基識(shí)別，應(yīng)用Seqtk軟件進(jìn)行誤差過濾；將3組雄蕊樣本的測(cè)序數(shù)據(jù)合并形成pool reads，然后應(yīng)用Trinity軟件[13-14]進(jìn)行denovo拼接；將拼接得到的final unigenes序列與NCBI官方非冗余蛋白質(zhì)序列數(shù)據(jù)庫(Nr)和UniProt數(shù)據(jù)庫進(jìn)行BLASTx比對(duì)(E≤1×10-5)[15]。

應(yīng)用rpstblastn程序(http：∥www.biomedcentral.com/content/supplementary/1471-2105-13-42-s1/Cloud- BioLinux-Package-Documentation/docs/rpstblastn.html)將unigenes序列與真核生物直系同源序列數(shù)據(jù)庫(KOG)進(jìn)行比對(duì)(E≤1×10-5)，取排名前5的比對(duì)結(jié)果進(jìn)行KOG功能分析[16-17]。將unigenes序列與基因功能分類體系數(shù)據(jù)庫(GO)進(jìn)行BLASTx比對(duì)(E≤1×10-5)，取最佳比對(duì)結(jié)果進(jìn)行GO功能分析[18]。

應(yīng)用KEGG KAAS在線pathway比對(duì)分析工具(https：∥www.genome.jp/kaas-bin/kaas_main)對(duì)unigenes進(jìn)行KEGG代謝途徑分析[19]。

2 結(jié)果和分析

2.1 測(cè)序結(jié)果和數(shù)據(jù)組裝結(jié)果分析

測(cè)序結(jié)果顯示：測(cè)序樣品的Q20值均大于96%，數(shù)據(jù)量約31.46 G，clean reads所占比例均在96.0%以上。說明轉(zhuǎn)錄組測(cè)序質(zhì)量較高，可以用于denovo拼接。

對(duì)clean reads進(jìn)行轉(zhuǎn)錄組裝和拼接，得到197 193個(gè)contigs，片段總長度為240 463 900 bp，最長片段為15 810 bp，非冗余unigenes片段的平均長度為1 219 bp，GC含量為39.1%，片段長度N50為1 966 bp。統(tǒng)計(jì)結(jié)果顯示：長度大于等于2 000 bp的unigenes有38 740個(gè)，占unigenes總數(shù)的19.6%；長度小于2 000 bp但大于等于1 000 bp 的unigenes有46 857個(gè)，占unigenes總數(shù)的23.8%；長度小于1 000 bp但大于等于200 bp 的unigenes有111 596個(gè)，占unigenes總數(shù)的56.6%。總體上長度大于等于1 000 bp的unigenes占unigenes總數(shù)的43.4%，表明組裝完整性較高，可用于注釋分析。

2.2 unigenes功能注釋分析

2.2.1 序列相似性分析 花生異型雄蕊轉(zhuǎn)錄組unigenes序列與Nr數(shù)據(jù)庫的比對(duì)結(jié)果(圖1)顯示：在花生異型雄蕊轉(zhuǎn)錄組中共有129 223個(gè)unigenes注釋到Nr數(shù)據(jù)庫，占unigenes總數(shù)的65.5%；其中，與Nr數(shù)據(jù)庫中的同源序列相比，E值較低(1×10-30

圖1 花生異型雄蕊轉(zhuǎn)錄組unigenes序列的E值分布圖

在Nr數(shù)據(jù)庫中對(duì)花生異型雄蕊轉(zhuǎn)錄組unigenes與其他植物種類轉(zhuǎn)錄組的注釋unigenes進(jìn)行同源序列匹配，結(jié)果見圖2。結(jié)果顯示：在花生異型雄蕊轉(zhuǎn)錄組中共有106 900個(gè)unigenes注釋到已探明的同源序列中，占注釋unigenes總數(shù)的82.7%。其中，與大豆〔Glycinemax(Linn.) Merr.〕的32 466個(gè)unigenes序列相似，占注釋unigenes總數(shù)的25.1%，同源性最高；與蒺藜苜蓿(MedicagotruncatulaGaertn.)的17 343個(gè)unigenes序列相似，占注釋unigenes總數(shù)的13.4%；與菜豆(PhaseolusvulgarisLinn.)的12 710個(gè)unigenes序列相似，占注釋unigenes總數(shù)的9.8%；與小豆〔Vignaangularisvar.angularis(Willd.) Ohwi et H. Ohashi〕的9 804個(gè)unigenes序列相似，占注釋unigenes總數(shù)的7.6%；與赤豆〔Phaseolusangularis(Willd.) Ohwi et H. Ohashi〕的3 958個(gè)unigenes序列相似，占注釋unigenes總數(shù)的3.1%；與花生(Nr數(shù)據(jù)庫)的3 771個(gè)unigenes序列相似，占注釋unigenes總數(shù)的2.9%；與百脈根(LotuscorniculatusLinn.)的3 273個(gè)unigenes序列相似，占注釋unigenes總數(shù)的2.5%；與葡萄(VitisviniferaLinn.)、桃(AmygdaluspersicaLinn.)和可可(TheobromacacaoLinn.)等13種植物的unigenes序列同源性均較低，在已注釋unigenes總數(shù)中所占比例均在2.0%以下。

圖2 花生異型雄蕊轉(zhuǎn)錄組unigenes與其他植物轉(zhuǎn)錄組unigenes的同源序列數(shù)量分布圖

2.2.2 KOG功能分析 在KOG數(shù)據(jù)庫中對(duì)花生異型雄蕊轉(zhuǎn)錄組unigenes序列進(jìn)行功能比對(duì)，結(jié)果見圖3。結(jié)果顯示：在花生異型雄蕊轉(zhuǎn)錄組unigenes序列中，共有133 387個(gè)unigenes注釋到25個(gè)KOG功能類別上，其中有8 776個(gè)unigenes注釋到未知功能(function unknown)，占注釋unigenes總數(shù)的6.6%，這些unigenes的功能仍不明確。

圖3 花生異型雄蕊轉(zhuǎn)錄組unigenes的KOG功能分類

功能明確的unigenes均與配子體發(fā)育功能有關(guān)。其中，有17 135個(gè)unigenes注釋到信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(signal transduction mechanism)功能，占注釋unigenes總數(shù)的12.8%；有16 407個(gè)unigenes注釋到一般功能預(yù)測(cè)(general function prediction only)，占注釋unigenes總數(shù)的12.3%；有9 893個(gè)unigenes注釋到翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(posttranslational modification，protein turnover and chaperone)功能，占注釋unigenes總數(shù)的7.4%；有8 870個(gè)unigenes注釋到轉(zhuǎn)錄(transcription)功能，占注釋unigenes總數(shù)的6.6%；有5 138個(gè)unigenes注釋到細(xì)胞周期調(diào)控、細(xì)胞分裂和染色體分離(cell cycle control，cell division and chromosome partitioning)功能，占注釋unigenes總數(shù)的3.9%；有2 710個(gè)unigenes注釋到染色體結(jié)構(gòu)和活力(chromatin structure and dynamic)功能，占注釋unigenes總數(shù)的2.0%。

2.2.3 GO功能分析 在GO數(shù)據(jù)庫中對(duì)花生異型雄蕊轉(zhuǎn)錄組unigenes序列進(jìn)行功能比對(duì)，結(jié)果見表1。結(jié)果顯示：在花生異型雄蕊轉(zhuǎn)錄組unigenes序列中，注釋到生物過程、細(xì)胞組分和分子功能3大類55個(gè)亞類的unigenes共有48 497個(gè)。

在生物過程的23個(gè)亞類中，分別有31 817、25 479、16 828、6 327、6 192和5 866個(gè)unigenes注釋到代謝過程(metabolic process)、細(xì)胞過程(cellular process)、單一生物過程(single-organism process)、生物調(diào)節(jié)(biological regulation)、應(yīng)激反應(yīng)(response to stimulus)和生物過程調(diào)節(jié)(regulation of biological process)功能，數(shù)量較多，各占注釋unigenes總數(shù)的65.6%、52.5%、34.7%、13.0%、12.8%和12.1%。其中，在核酸代謝過程(GO： 0090304，8334個(gè))、磷酸化(GO： 0016310，4428個(gè))和轉(zhuǎn)錄調(diào)節(jié)(GO：0006355，2322個(gè))等與代謝過程相關(guān)的GO功能類別中unigenes高度富集；而注釋到節(jié)律過程(rhythmic process)、運(yùn)動(dòng)(locomotion)、排毒(detoxification)、生物附著(biological adhesion)和行為(behavior)功能的unigenes數(shù)量則較少。

在細(xì)胞組分的18個(gè)亞類中，分別有13 334、13 334、9 427和8 905個(gè)unigenes注釋到細(xì)胞(cell)、細(xì)胞成分(cell part)、細(xì)胞器(organelle)和細(xì)胞膜(cell membrane)功能，數(shù)量較多，各占注釋unigenes總數(shù)的27.5%、27.5%、19.4%和18.4%。其中，在與核(GO：0005634，4 069個(gè))和核糖體(GO：0005840，639個(gè))等細(xì)胞器相關(guān)的GO功能類別中unigenes高度富集；而注釋到類核(nucleoid)、病毒體部分(virion part)、病毒(virion)、膜外區(qū)域部分(extracellular region part)、胞外基質(zhì)(extracellular matrix)和胞外基質(zhì)組分(extracellular matrix component)功能的unigenes數(shù)量則較少。

表1 花生異型雄蕊轉(zhuǎn)錄組unigenes的GO功能分類

在分子功能的14個(gè)亞類中，分別有28 066和26 091個(gè)unigenes注釋到結(jié)合(binding)和催化活性(catalytic activity)功能，數(shù)量較多，各占注釋unigenes總數(shù)的57.9%和53.8%。其中，在ATP結(jié)合(GO：0005524，7 229個(gè))、轉(zhuǎn)錄因子(GO：0003700，1 161個(gè))和鈣離子結(jié)合(GO：0005509，728個(gè))等與結(jié)合相關(guān)的GO功能類別中unigenes高度富集；而注釋到營養(yǎng)儲(chǔ)存活性(nutrient reservoir activity)、轉(zhuǎn)錄調(diào)控因子活性(translation regulator activity)和金屬伴侶蛋白活性(metallochaperone activity)功能的unigenes數(shù)量則較少。

2.2.4 KEGG代謝通路分析 以KEGG代謝庫為參考，對(duì)花生異型雄蕊轉(zhuǎn)錄組unigenes的代謝通路進(jìn)行統(tǒng)計(jì)和分類，結(jié)果見圖4。結(jié)果顯示：在花生異型雄蕊轉(zhuǎn)錄組中共有30 428個(gè)unigenes注釋到2 708個(gè)酶功能，映射了354個(gè)代謝通路，可分為遺傳信息處理(genetic information processing)、細(xì)胞過程(cellular process)、生物系統(tǒng)(organismal system)、環(huán)境信息處理(environmental information processing)和代謝 (metabolism) 5大類32個(gè)亞類，其中與遺傳信息處理、細(xì)胞過程和環(huán)境信息處理相關(guān)的亞類較少，分別僅有4、4和3個(gè)。

在遺傳信息處理代謝通路中，有3 307個(gè)unigenes注釋到翻譯(transcription)亞類，有3 015個(gè)unigenes注釋到折疊、分類和降解(folding，sorting and degradation)亞類。在細(xì)胞過程代謝通路中，有2 713個(gè)unigenes注釋到運(yùn)輸和分解代謝(transport and catabolism)亞類。在環(huán)境信息處理代謝通路中，有4 911個(gè)unigenes注釋到信號(hào)轉(zhuǎn)導(dǎo)(signal transduction)亞類，其中，有1 055個(gè)unigenes注釋到植物激素信號(hào)轉(zhuǎn)導(dǎo)通路(ko04075)，有285個(gè)unigenes注釋到鈣信號(hào)通路(ko01100)。

： 遺傳信息處理Genetic information processing； ： 細(xì)胞過程Cellular process； ： 生物系統(tǒng)Organismal system； ： 環(huán)境信息處理Environmental information processing； ： 代謝Metabolism.

共有5 475個(gè)unigenes注釋到與生物系統(tǒng)相關(guān)的代謝通路中，占注釋unigenes總數(shù)的18.0%，分別注釋到9個(gè)亞類。其中，有1 913個(gè)unigenes注釋到內(nèi)分泌系統(tǒng)(endocrine system)亞類，有1 453個(gè)unigenes注釋到免疫系統(tǒng)(immune system)亞類。

注釋到與代謝相關(guān)通路的unigenes數(shù)量最多，有14 231個(gè)，占注釋unigenes總數(shù)的46.8%，分別注釋到12個(gè)亞類。其中，有6 529個(gè)unigenes注釋到全球和總覽圖(global and overview maps)亞類，有4 635個(gè)unigenes注釋到糖類代謝(carbohydrate metabolism)亞類，有2 802個(gè)unigenes注釋到氨基酸代謝(amino acid metabolism)亞類，有2 538個(gè)unigenes注釋到脂質(zhì)代謝(lipid metabolism)亞類，有2 337個(gè)unigenes注釋到能量代謝(energy metabolism)亞類，有1 278個(gè)unigenes注釋到核苷酸代謝(nucleotide metabolism)亞類。

3 討論和結(jié)論

花生隸屬于豆科(Fabaceae)落花生屬(ArachisLinn.)，從同源序列的匹配結(jié)果看，其異型雄蕊轉(zhuǎn)錄組unigenes與同科的大豆和蒺藜苜蓿的同源匹配率較高，與大豆的32 466個(gè)unigenes序列相似，與蒺藜苜蓿的17 343個(gè)unigenes序列相似，這一結(jié)果佐證了這些種類親緣關(guān)系。而供試花生品種的異型雄蕊轉(zhuǎn)錄組中僅有3 771個(gè)unigenes與Nr數(shù)據(jù)庫中花生的同源unigenes相匹配，一方面可能是由于本實(shí)驗(yàn)是以花生異型雄蕊為實(shí)驗(yàn)材料，而Nr數(shù)據(jù)庫中的花生基因組注釋信息是以果針為研究材料有關(guān)；另一方面，花生栽培種是由2個(gè)二倍體野生種A.duranensis和A.ipa?nsis自然雜交后經(jīng)過染色體加倍形成的異源四倍體[20-21]，導(dǎo)致花生栽培品種的遺傳關(guān)系復(fù)雜。

有花植物的雄蕊和花粉發(fā)育是一個(gè)復(fù)雜的過程，與花被器官相比，生殖器官的特有基因更多[22]。對(duì)花生異型雄蕊轉(zhuǎn)錄組測(cè)序獲得數(shù)據(jù)量約31.46 G，其中，有129 223個(gè)unigenes在Nr數(shù)據(jù)庫比對(duì)出同源序列信息，且許多unigenes功能與雄蕊產(chǎn)生雄配子的功能有關(guān)。許多轉(zhuǎn)錄因子決定著花器官中萼片、花瓣、雄蕊和心皮[23]的發(fā)育時(shí)序與數(shù)量，進(jìn)而參與控制被子植物花發(fā)育。在花生異型雄蕊轉(zhuǎn)錄組中，有8 870個(gè)unigenes注釋到轉(zhuǎn)錄功能，其中5 360個(gè)unigenes具有轉(zhuǎn)錄因子功能，說明轉(zhuǎn)錄因子對(duì)植物基因表達(dá)的調(diào)控有至關(guān)重要的作用[24]。值得注意的是，有6.6%的unigenes基因功能仍不明確，這與花生的基因組大、重復(fù)序列比例高等原因有關(guān)。

雄蕊和花粉的發(fā)育涉及一系列精密的細(xì)胞質(zhì)和核基因相互作用[25]；Ca2+可作為細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)的信使，通過鈣調(diào)蛋白(CaM)將胞外信號(hào)轉(zhuǎn)換成胞內(nèi)的生理生化反應(yīng)，包括細(xì)胞分裂、分化和凋亡，對(duì)花粉發(fā)育和萌發(fā)起著非常重要的作用[26]。在GO數(shù)據(jù)庫中對(duì)花生異型雄蕊轉(zhuǎn)錄組unigenes序列進(jìn)行功能比對(duì)，結(jié)果顯示其unigenes在與核酸代謝過程、磷酸化、轉(zhuǎn)錄調(diào)節(jié)、核、核糖體、ATP結(jié)合、轉(zhuǎn)錄因子和鈣離子結(jié)合相關(guān)的功能類別中高度富集；而KEGG代謝通路分析結(jié)果表明：在花生異型雄蕊轉(zhuǎn)錄組中注釋到代謝相關(guān)通路的unigenes數(shù)量最多，并以注釋到糖類代謝通路、氨基酸代謝通路、脂質(zhì)代謝通路、能量代謝通路和核苷酸代謝亞類的unigenes數(shù)量較多，說明有大量涉及孢粉素合成、氨基酸合成、核糖體組裝、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工和脂質(zhì)轉(zhuǎn)移等功能的unigenes被注釋；此外，與鈣調(diào)蛋白相關(guān)的unigenes也被注釋到無機(jī)離子運(yùn)輸與代謝以及信號(hào)轉(zhuǎn)導(dǎo)機(jī)制這2個(gè)類別上。由此可知，在花生異型雄蕊轉(zhuǎn)錄組中，一部分基因可參與不同的代謝途徑，同時(shí)花藥發(fā)育的精密調(diào)控需要多個(gè)代謝途徑參與，其中，糖類代謝不僅為花藥發(fā)育提供能量，還可作為信號(hào)物質(zhì)影響花藥和花粉的發(fā)育[27]；而在花粉發(fā)育過程中通常會(huì)觀察到高呼吸率和高能量需求[28]，主要代謝途徑的功能障礙將對(duì)花粉粒發(fā)育產(chǎn)生不利影響[29]，雄性不育系中與能量供應(yīng)相關(guān)的ATP含量顯著降低[30-31]。因而，對(duì)于花生異型雄蕊，特別是退化雄蕊與可育雄蕊進(jìn)行差異表達(dá)基因功能分析，可為花生雄蕊發(fā)育的分子調(diào)控提供更多啟示。

花生異型雄蕊的發(fā)育調(diào)控是一個(gè)復(fù)雜的生物過程，應(yīng)結(jié)合花生基因組測(cè)序成果開展大量的基因組分析，以期挖掘更多有價(jià)值的遺傳信息，從而更全面地了解花生異型雄蕊分化的分子調(diào)控機(jī)制。