基于高通量測序的半夏變異珠芽轉(zhuǎn)錄組分析

馬琛，汪雷，葉德，徐濤

(浙江理工大學(xué) 生命科學(xué)學(xué)院,浙江 杭州 310018)

半夏Pinelliaternata(Thunb.)Breit，多年生草本植物，天南星科半夏屬，在我國，其塊莖是高利用的中藥之一；據(jù)藥典記錄，其具有燥濕化痰，降逆止嘔，消痞散結(jié)功效，常用于痰多咳喘，風(fēng)痰眩暈，痰厥頭痛，嘔吐反胃，梅核氣等。作為重要的傳統(tǒng)中藥，半夏常以干燥塊莖入藥，最早記載于《神農(nóng)本草經(jīng)》，列為下品，其性溫、味辛、略有毒，歸脾、胃、肺；具有外用止痛消腫及抗腫瘤潛力[1-4]。珠芽繁殖是半夏繁殖的重要繁殖方式，其具有快速繁殖，產(chǎn)量高等優(yōu)點(diǎn)，在自然界和農(nóng)業(yè)生產(chǎn)中能增加半夏的存活率。因此，珠芽品質(zhì)的優(yōu)劣是決定半夏質(zhì)量與產(chǎn)量的重要因子，也極大影響了半夏的種植。由于野生的半夏資源已越來越稀少，而半夏的需求卻日漸增大，因此，為了培育優(yōu)質(zhì)高產(chǎn)的半夏品種，利用基因工程手段就成為緩解半夏資源緊缺的解決手段之一。目前的半夏研究主要集中在優(yōu)質(zhì)品種的篩選和高產(chǎn)量、藥理藥效及臨床等方面，對半夏珠芽功能基因的研究較少、進(jìn)展緩慢。

篩選差異表達(dá)基因(differential expressed gene, DEG)可針對性地分析功能基因，利用這些基因?qū)ξ锓N進(jìn)行定向遺傳育種改良，提高抗脅迫能力及其分子機(jī)制的深入研究[5-6]?；虿町惐磉_(dá)是指植株在特定環(huán)境下，發(fā)生生理應(yīng)激反應(yīng)引起功能蛋白的差異表達(dá)，導(dǎo)致生長表型出現(xiàn)差異的現(xiàn)象。轉(zhuǎn)錄組(transcriptome)是指一定狀態(tài)下，細(xì)胞內(nèi)所有轉(zhuǎn)錄產(chǎn)物的集合，轉(zhuǎn)錄組測序能夠全面快速地獲取某一物種特定器官或組織在某一狀態(tài)下的幾乎所有轉(zhuǎn)錄體，并反映其表達(dá)水平[7-9]。利用轉(zhuǎn)錄組分析能夠批量地篩選發(fā)掘無參考基因組物種的功能基因，以及探索在不同樣品中基因的表達(dá)差異。本研究借助高通量測序技術(shù)[10]，比較分析半夏葉基部長珠芽和不長珠芽的DEG，并篩選出影響半夏珠芽激素生物合成和信號傳導(dǎo)的DEG，為深入研究半夏珠芽相關(guān)功能基因的克隆及功能分析提供基礎(chǔ)。

1 材料與方法

1.1 材料

三葉半夏種子取自甘肅天水，種于浙江理工大學(xué)6號樓317實(shí)驗(yàn)室內(nèi)；選取發(fā)育正常、長勢優(yōu)良且無病蟲害侵襲的植株進(jìn)行采樣，用解剖刀分別刮取半夏葉基部相同部位長珠芽和不長珠芽的組織，用蒸餾水將其清洗干凈，立即用液氮冷凍，放入-80 ℃冰箱保存?zhèn)溆?。樣品與分組信息：半夏葉基部長珠芽為樣品BY，不長珠芽的為樣品ZC。本實(shí)驗(yàn)利用AxyPrep總RNA小量制備試劑盒，按步驟提取總RNA，并進(jìn)行質(zhì)量檢測。經(jīng)質(zhì)檢合格后，使用連接有Oligo(dT)的磁珠富集真核生物mRNA。片段化試劑(fragmentation buffer)將mRNA隨機(jī)打斷成短片段，以片段化的mRNA為模板，利用隨機(jī)引物(random hexamers)合成一鏈cDNA，隨后合成二鏈cDNA。AMPure XP beads純化雙鏈產(chǎn)物，利用T4 DNA聚合酶和Klenow DNA聚合酶將DNA的黏性末端修復(fù)為平末端，3’末端加堿基A并加接頭，AMPureXP beads進(jìn)行片段選擇，最后進(jìn)行PCR擴(kuò)增獲得最終測序文庫。文庫質(zhì)檢合格后用Illumina HiseqTM 4000進(jìn)行測序，開展轉(zhuǎn)錄組研究，共獲得單基因(unigene)28 395條。本研究在此基礎(chǔ)上分析BY和ZC樣品的DEGs。

1.2 DEG的篩選條件

通過RPKM(reads per kilo base of exon model per million mapped reads)法來計算樣品BY和樣品ZC的unigene表達(dá)的豐度值，表示在每百萬測序堿基中每千個轉(zhuǎn)錄子測序堿基中所包含的測序reads數(shù)。計算公式：

RPKM=total exon reads÷mapped reads (million)÷exonlenth(KB)。

本研究根據(jù)公式計算變異珠芽和正常珠芽之間表達(dá)差異的倍數(shù)，用log2fold_change表示樣本基因表達(dá)的相對差異倍數(shù)，計算方法為：log2(Sample_BY FPKM/Sample_ZC FPKM)，即樣本BY除以ZC的基因表達(dá)水平。將RPKM的log2倍數(shù)絕對值>1且FDR(false discovery rate，錯誤發(fā)現(xiàn)率)≤0.05的unigene篩選出，即為DEG。

1.3 DEG的功能注釋

將篩選到的unigene序列分別映射到公共數(shù)據(jù)庫Swiss-Prot、NR、Pfam、KEGG、COG、GO中，通過BLASTX得到與指定unigene具有序列最高相似性(取閾值e<=1e-10)的蛋白(表1)，進(jìn)一步得到相應(yīng)unigene所注釋的基因及其蛋白功能信息。比對到公共數(shù)據(jù)庫NR的結(jié)果顯示，高達(dá)71.22%的unigenes序列找到了相似序列。而近緣物種在相似性序列的匹配中，海棗所占比例最高(44.2%)，其他依次是葡萄(37.5%)、可可樹(8.2%)、無油樟(3.2%)、大桉(2.6%)、桐油樹(2.2%)和其他物種(2.1%)。

表1 BLAST注釋結(jié)果統(tǒng)計情況

將篩選到的半夏unigenes比對到數(shù)據(jù)庫COG里，結(jié)果顯示，共有17 130條unigenes(占總數(shù)60.33%)被注釋到25種COG分類中(表2)。被注釋的unigenes涵蓋大部分生命活動，涉及KOG功能類別也比較全面；其中，最大的分類是“基因功能預(yù)測”，包含2 063條(12.04%)unigenes；其次是“信號傳導(dǎo)機(jī)制”，包含了1 788條(10.44%)unigenes；而“細(xì)胞運(yùn)動”是最少的，僅4條(0.02%)unigenes；在其他類別中，各自的基因表達(dá)豐度都不相同。

表2 COG功能分類

2 結(jié)果與分析

2.1 DEG篩選

De novo測序中最重要的過程就是基因差異表達(dá)分析，通過差異表達(dá)分析可以了解不同樣本間基因的差異表達(dá)情況。根據(jù)條件|Log2(樣品BY RPKM/樣品ZC RPKM)|≥1、且FDR≤0.05，共篩選到13 656條DEG，其中包含6 139條上調(diào)表達(dá)的DEG和7 517條下調(diào)表達(dá)的DEG。

2.2 DEG 功能分析

將篩選所得的13 656條DEGs導(dǎo)入數(shù)據(jù)庫GO中，深入探索其功能定義。結(jié)果表明，在GO數(shù)據(jù)庫中獲得注釋的DEG共有4 595條，占DEG總量的33.65%。其中分別將2 324、592和1 679條DEG歸類至三大主要功能類別，即生物學(xué)途徑(biological process)、細(xì)胞組件(cellular component)和分子生物學(xué)功能(molecular function)。將被注釋的DEG功能定義可進(jìn)一步分在3個主要功能類別的50個亞類中(圖1)，在“生物學(xué)途徑”類別中，“依賴DNA的轉(zhuǎn)錄調(diào)節(jié)(regulation of transcription, DNA-dependent)”是DEG歸入最多的亞類；最少的則是“多細(xì)胞有機(jī)體的發(fā)育(multicellular organismal development)”。在“細(xì)胞組件”類別中，DEG歸入最多的亞類是“膜的完整性(integral to membrane)”；最少的是“葉綠體類囊體膜(chloroplast thylakoid membrane)”。在“分子生物學(xué)功能”類別中，DEG歸入最多的亞類是“ATP結(jié)合(ATP binding)”；最少的是“GTP結(jié)合(GTP binding)”。

圖1 DEG的GO分類

2.3 DEG KEGG富集分析

生物體內(nèi)，生物學(xué)功能的行使依賴于不同基因的相互協(xié)調(diào)，通過分析代謝途徑可以確定DEG的主要代謝途徑及信號通路等。獲得的13 656條DEG比對到KEGG代謝途徑數(shù)據(jù)庫，獲得功能注釋的DEG共3 551條，占DEG總數(shù)的26.00%，包含在239條KEGG代謝途徑中。在KEGG途徑的功能定義中，主要包括產(chǎn)物代謝與合成及細(xì)胞作用等，表3中是富集DEG最多的10條KEGG代謝途徑。

2.4 影響激素合成和信號傳導(dǎo)的DEGs分析

將篩選到的DEGs與Nr、GO、KEGG等蛋白公共數(shù)據(jù)庫分別比對，得到DEGs的功能注釋，并篩選到參與半夏珠芽激素生物合成及信號傳導(dǎo)的關(guān)鍵酶基因17個(表4)。通過分析該基因在半夏變異植株和正常植株中的RPKM值顯示，表達(dá)上調(diào)基因?yàn)?個，即在BY中的表達(dá)量高于在ZC中的表達(dá)量；表達(dá)下調(diào)基因?yàn)?1個，即在ZC中的表達(dá)量高于BY中的表達(dá)量。

表3 半夏DEG富集最多的10條KEGG代謝途徑

3 討論

植物生長中的發(fā)育差異及生理結(jié)構(gòu)的變化甚至病變，細(xì)胞功能多樣性和形態(tài)差異的產(chǎn)生等決定性因子是基因的差異性表達(dá)[11]。因此，分析基因差異性表達(dá)，可有效挖掘影響植物性狀的相關(guān)功能基因，促進(jìn)植物遺傳育種的改良。20世紀(jì)90年代起至今，差減雜交、mRNA差異顯示PCR、抑制性消減雜交、基因表達(dá)的系列分析(SAGE)、表達(dá)序列標(biāo)簽(EST)、cDNA限制性片段長度多態(tài)性分析、基因芯片和熒光定量PCR等經(jīng)典基因表達(dá)差異分析方法已經(jīng)形成[12]。但這些方法都各自存在一定的不足，如差減雜交經(jīng)過層析，所回收的cDNA量很低，嚴(yán)重影響實(shí)驗(yàn)的進(jìn)一步進(jìn)行，且只能檢測2組mRNA之間的表達(dá)差異；mRNA差異顯示逆轉(zhuǎn)錄PCR需反復(fù)操作，進(jìn)行多次電泳及PCR。本研究通過高通量測序技術(shù)挖掘獲得相關(guān)功能基因，并進(jìn)行基因新功能的探討；采用自動化程序，同時分析多組mRNA，檢測出存在差異的單個核苷酸。分析DEG的同時，還能獲得SSR[13]信息及組裝出的完整基因序列。本研究從半夏珠芽轉(zhuǎn)錄組數(shù)據(jù)的28 395條unigenes中，高效批量地篩選出DEGs 13 656條；通過GO歸類，該DEGs被包含在50個功能亞類中，顯著集中在半夏生長發(fā)育過程中的性狀表現(xiàn)及次生代謝等。分析KEGG途徑顯示，有3 551條半夏珠芽差異表達(dá)基因富集在生理代謝過程和次生代謝產(chǎn)物合成途徑中，占DEGs總量的26.00%。本研究獲得的DEGs數(shù)據(jù)豐富了半夏珠芽的基因資源，對研究半夏珠芽的基因調(diào)控機(jī)制意義重大。

表4 半夏珠芽激素合成與信號轉(zhuǎn)導(dǎo)相關(guān)DEG

植物生長發(fā)育過程中，參與調(diào)控功能的內(nèi)源性激素非常重要，其調(diào)節(jié)生物合成及信號傳導(dǎo)的酶基因也成為研究熱點(diǎn)。基于半夏珠芽轉(zhuǎn)錄組DEGs分析，本研究共篩選獲得參與半夏珠芽內(nèi)源性激素生物合成和信號傳導(dǎo)相關(guān)過程的功能基因17個，這些基因在半夏葉基部長珠芽和不長珠芽中呈現(xiàn)出明顯的表達(dá)差異，因此，半夏珠芽生長發(fā)育過程與這些基因及其調(diào)節(jié)的代謝產(chǎn)物密切相關(guān)。

生長素即吲哚乙酸，是最早發(fā)現(xiàn)的促進(jìn)植物生長的激素，但對莖、芽、根生長的促進(jìn)作用因濃度而異，調(diào)控細(xì)胞分裂和分化及頂端優(yōu)勢等一系列生理過程[14]。本研究表明，在ZC中，半夏生長素調(diào)控基因低表達(dá)，而在BY中高表達(dá)，推斷在BY中生長素積累高于ZC，這可能是導(dǎo)致半夏葉基部長珠芽的因素之一。赤霉素最顯著的功能是加速細(xì)胞的伸長，促進(jìn)細(xì)胞分裂和擴(kuò)大，且不引起細(xì)胞壁的酸化，除此之外，赤霉素還有著抑制成熟，側(cè)芽休眠，衰老，塊莖形成的生理作用[15]。本研究從半夏珠芽轉(zhuǎn)錄組中篩選到的GA 2-oxidase，是赤霉素氧化酶的編碼基因，該基因主要表達(dá)于ZC 中，表明半夏珠芽赤霉素合成具有特異性。細(xì)胞分裂素(cytokinin kinetin)是一類調(diào)節(jié)植物細(xì)胞生長和發(fā)育的植物激素，與植物生長素有協(xié)同作用；從玉米或其他植物中分離得到，一般在植物根部產(chǎn)生。其有兩種明顯的生理作用：一是促進(jìn)細(xì)胞分裂和調(diào)控其分化，二是延緩蛋白質(zhì)和葉綠素的降解，延遲衰老。本研究從半夏轉(zhuǎn)錄組中鑒定到的乙酰輔酶A羧化酶/生物素羧化酶[16]，在半夏BY中明顯高于半夏ZC中。脫落酸是一種抑制生長的植物激素，廣泛分布于高等植物，除促使葉子脫落外，還可使芽進(jìn)入休眠狀態(tài)，抑制細(xì)胞的分裂；因此可促進(jìn)葉和果實(shí)的衰老和脫落，增強(qiáng)植物抗逆性[17-18]。本研究篩選到9-順式-環(huán)氧類胡蘿卜素雙加氧酶[19]，是一個ABA生物合成的限速關(guān)鍵酶基因，密切影響ABA的積累，且滲透脅迫誘導(dǎo)合成的關(guān)鍵性調(diào)節(jié)位點(diǎn)即為該基因催化的步驟。相比BY中，9-順式-環(huán)氧類胡蘿卜素雙加氧酶在ZC表達(dá)更高，顯示ABA合成在ZC中更為活躍。乙烯在種子萌發(fā)到成熟衰老和植物花、果實(shí)的成熟及衰老等一系列生理過程，起重要的調(diào)節(jié)作用，間接影響植物的性狀差異及品質(zhì)優(yōu)劣[20]。研究挖掘到調(diào)控乙烯生物合成的重要關(guān)鍵酶基因——1-氨基環(huán)丙烷-1-羧酸合酶[21]，是催化乙烯生物合成的限速酶，將S-腺苷甲硫氨酸轉(zhuǎn)變成1-氨基環(huán)丙烷基-1-羧酸和甲基硫腺苷。其在半夏BY中表達(dá)下調(diào)，在半夏ZC中表達(dá)高于在BY中，說明乙烯主要在ZC中合成明顯。油菜素甾醇被稱為第六類植物激素，在植物中含量低微，僅有其他五大激素任意一種的千分之一，但其具有的高生理活性、廣譜、無毒及獨(dú)特作用機(jī)理，對植物的增產(chǎn)和抗逆性具有重大意義[22]。并且生長素和油菜素甾醇通過共同調(diào)控下游基因的轉(zhuǎn)錄、二者信號組分間的互作、合成代謝的調(diào)控以及調(diào)控生長素極性運(yùn)輸?shù)?，在多層次相互作用從而對植物生長發(fā)育實(shí)現(xiàn)精確地調(diào)控[23]。植物中油菜素甾醇的信號轉(zhuǎn)導(dǎo)及其與脫落酸信號的相互調(diào)節(jié)對半夏的生長發(fā)育也具有重要影響。茉莉酸在種子的萌發(fā)、果實(shí)的成熟、根的生長、球莖的形成以及孕育花粉等生長發(fā)育相關(guān)的生理過程，同時在系統(tǒng)防御生物過程起著關(guān)鍵的調(diào)控作用[24]。本研究中茉莉酸生物合成在半夏BY中明顯高于半夏ZC，可以推斷，半夏BY植株在生長過程中有更多的優(yōu)勢，對于半夏產(chǎn)量和質(zhì)量的提高都大有裨益。

本研究利用高通量測序技術(shù)對三葉半夏葉基部長珠芽和不長珠芽進(jìn)行轉(zhuǎn)錄組和基因組分析,豐富了半夏基因組的信息,并挖掘到大量差異表達(dá)基因，為后續(xù)半夏功能基因的分析及克隆提供了豐富的序列數(shù)據(jù)基礎(chǔ)。本研究將篩選到的半夏unigenes比對到數(shù)據(jù)庫COG里，結(jié)果顯示：共有17 130條unigenes被注釋到25種KOG分類中，其中721條基因功能未知。希望未來更多物種可以進(jìn)行轉(zhuǎn)錄組研究，在豐富基因公共數(shù)據(jù)庫的同時，發(fā)掘已知基因的新功能和新的功能基因。本研究篩選得到13 656條顯著性差異表達(dá)基因，需要對其進(jìn)行深入研究，發(fā)掘出半夏生長發(fā)育中控制性狀的關(guān)鍵基因，探索生物學(xué)調(diào)控機(jī)制，為半夏優(yōu)質(zhì)高產(chǎn)品種的改良提供科學(xué)依據(jù)，促進(jìn)中藥半夏的深度開發(fā)。