三葉木通雌雄花芽轉(zhuǎn)錄組測(cè)序及分析

田 天 項(xiàng) 明 陳發(fā)菊

(三峽大學(xué) 生物技術(shù)研究中心，湖北 宜昌 443002)

三葉木通(Akebiatrifoliate)是木通科(Lardizabalaceae)木通屬(Akebia Decne.)落葉木質(zhì)藤本植物，民間又稱(chēng)八月炸、八月瓜.三葉木通花序?yàn)榭偁罨ㄐ?，下部生?～2朵雌花，上部開(kāi)有15～30朵雄花.三葉木通早期均為兩性花結(jié)構(gòu)，在發(fā)育的后期雌蕊和雄蕊會(huì)選擇性退化而形成單性花[1].

高通量測(cè)序技術(shù)是近些年來(lái)發(fā)展出的新技術(shù)，廣泛用于SSR標(biāo)記的開(kāi)發(fā)鑒定[2]，各種生理活動(dòng)的代謝途徑以及基因分析[3]，也被廣泛用于性別相關(guān)基因的挖掘，韋麗君等[4]發(fā)現(xiàn)，生長(zhǎng)素、WRKY和MYB轉(zhuǎn)錄因子相關(guān)的基因與木薯性別決定相關(guān).研究表明簸箕柳的轉(zhuǎn)錄測(cè)序分析發(fā)現(xiàn)33條差異表達(dá)基因定位到十九號(hào)染色體上[5]，王萌通過(guò)轉(zhuǎn)錄組測(cè)序發(fā)現(xiàn)AG與FT基因可能對(duì)核桃的開(kāi)花調(diào)控具有重要作用[6].何曉琳通過(guò)對(duì)蓖麻的不同種系的轉(zhuǎn)錄組對(duì)比分析確定了雜種優(yōu)勢(shì)的基因及代謝通路[7].閆冬春在南瓜轉(zhuǎn)錄組找到了一個(gè)影響雌雄花發(fā)育的關(guān)鍵基因CmACO1[8].在山核桃的雌雄轉(zhuǎn)錄組分析中，研究者聚焦到了MADS-box基因家族，推測(cè)AP1和SEP亞家族基因可能與雌雄花的形成有緊密聯(lián)系.這些研究充分說(shuō)明利用轉(zhuǎn)錄組測(cè)序來(lái)獲得植物性別分化的相關(guān)基因是可行的.

本實(shí)驗(yàn)通過(guò)三葉木通的雌雄花芽進(jìn)行轉(zhuǎn)錄組測(cè)序，利用現(xiàn)有的數(shù)據(jù)庫(kù)與測(cè)序結(jié)果進(jìn)行比對(duì)，并對(duì)雌雄花芽的差異基因進(jìn)行驗(yàn)證分析，為探究三葉木通性別分化機(jī)制奠定基礎(chǔ)及利用基因工程培育新品種提供方向.

1 材料與方法

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

測(cè)序材料三葉木通雌雄花芽采自湖北省神農(nóng)架自然保護(hù)區(qū).以三葉木通雌花芽為對(duì)照組，記為MT1，以三葉木通雄花芽為實(shí)驗(yàn)組，記為MT2.實(shí)驗(yàn)材料采集之后立刻過(guò)氮處理，并保存于-80℃冰箱，轉(zhuǎn)錄組測(cè)序樣本重復(fù)3次，分別為MT1-1，MT1-2，MT1-3，MT2-1，MT2-2和MT2-3.

1.2 三葉木通總RNA提取、文庫(kù)制備及測(cè)序

采用Trizol試劑盒(Invirtrogen,CA,USA)提取三葉木通雄花芽和雌花芽總RNA，紫外分光光度計(jì)和Nanodrop軟件檢測(cè)RNA濃度與純度，1%瓊脂糖凝膠電泳檢測(cè)總RNA質(zhì)量和完整性.

按照Illumina公司的文庫(kù)制備實(shí)驗(yàn)流程制備文庫(kù)，樣本總RNA使用連接有Oligo(dT)的磁珠富集真核生物mRNA，然后被隨機(jī)打斷成短片段，用隨機(jī)引物法合成cDNA第一鏈，再合成cDNA第二鏈，最后純化回收雙鏈產(chǎn)物，利用T4 DNA連接酶和KlenowDNA聚合酶進(jìn)行進(jìn)一步處理，最后用USER酶降解cDNA第二鏈進(jìn)行PCR擴(kuò)增獲得測(cè)序文庫(kù).文庫(kù)質(zhì)檢合格后采用Illumina NovaseqTM6000 進(jìn)行測(cè)序，測(cè)序讀長(zhǎng)為雙端2*150 bp(PE150).

1.3 質(zhì)量控制、組裝統(tǒng)計(jì)、功能注釋與表達(dá)分析

首先，使用Cutadapt[9]和perl腳本去除低質(zhì)量堿基和未確定堿基，然后使用FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)通過(guò)有效數(shù)據(jù)的Q20，Q30和GC含量驗(yàn)證序列質(zhì)量，后續(xù)分析均基于高質(zhì)量的有效數(shù)據(jù).使用Trinity 2.4.0[10]進(jìn)行轉(zhuǎn)錄組序列的組裝，將轉(zhuǎn)錄本中的序列稱(chēng)為transcript.選擇其中最長(zhǎng)的轉(zhuǎn)錄本作為Unigene.組裝好的Unigenes通過(guò)DIAMOND[11]軟件與各個(gè)數(shù)據(jù)庫(kù)進(jìn)行對(duì)比和功能注釋.

1.4 三葉木通雌雄花芽差異基因的篩選

根據(jù)Patro[12]的方法計(jì)算TPM[13]衡量所有Unigenes表達(dá)水平，利用package edgeR[14]對(duì)Unigenes進(jìn)行差異表達(dá)分析，然后通過(guò)log2>1(log2 fold change值，差異倍數(shù))或者log2<-1并且具有統(tǒng)計(jì)學(xué)意義(P<0.05)進(jìn)行篩選，獲得三葉木通雌雄花芽中的差異基因.

1.5 三葉木通雌雄花芽差異表達(dá)基因的qRT-PCR驗(yàn)證

使用Trizol reagent試劑盒提取三葉木通總RNA，采用逆轉(zhuǎn)錄法(Takara)獲得cDNA.依據(jù)轉(zhuǎn)錄組測(cè)序結(jié)果，隨機(jī)篩選10條差異表達(dá)明顯的基因，進(jìn)一步分析篩選得到細(xì)胞分裂素以及丙酮酸相關(guān)的基因9個(gè).通過(guò)實(shí)時(shí)熒光定量PCR儀(CFX96，BIO-RAD)以及SYBR qPCR試劑盒(南京諾唯贊生物科技有限公司)進(jìn)行qRT-PCR實(shí)驗(yàn).通過(guò)NCBI BLAST設(shè)計(jì)qPCR引物，由華大基因股份有限公司合成.qRT-PCR引物見(jiàn)表1，產(chǎn)物長(zhǎng)度100～250 bp.選擇RG5(EF-α)作為內(nèi)參基因[15]，上游引物為GTCTCCCTCTTCAGGATGTTTAC，下游引物為GACAACCATACCAGGCTTGA.其中前10對(duì)引物用來(lái)驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)的可靠性，后9對(duì)引物為代謝通路引物.

表1 實(shí)時(shí)熒光PCR引物

2 結(jié)果與分析

2.1 三葉木通轉(zhuǎn)錄組數(shù)據(jù)測(cè)序及序列組裝

采用Illumina測(cè)序平臺(tái)技術(shù)對(duì)三葉木通樣本進(jìn)行轉(zhuǎn)錄組測(cè)序，共得到49.84 Gb的原始數(shù)據(jù)(raw data).過(guò)濾掉不合格的序列，共得到49.06 Gb有效數(shù)據(jù)，占原始數(shù)據(jù)的98.43%，組裝后GC含量與長(zhǎng)度分布如圖1～圖2所示，每組數(shù)據(jù)質(zhì)量見(jiàn)表2.對(duì)過(guò)濾后的數(shù)據(jù)進(jìn)行Trinity組裝，去除冗余后得到17 439個(gè)Unigene，通過(guò)六大數(shù)據(jù)庫(kù)注釋?zhuān)谢蚓灰粋€(gè)或者更多的數(shù)據(jù)庫(kù)注釋到(見(jiàn)表3).

圖1 Unigene組裝長(zhǎng)度統(tǒng)計(jì) 圖2 Unigene的GC含量統(tǒng)計(jì)

表2 過(guò)濾前后的Reads統(tǒng)計(jì)

表3 Unigene功能注釋統(tǒng)計(jì)

去除測(cè)序接頭和過(guò)濾掉不合格序列后得到有效數(shù)據(jù)(clean data)，去除測(cè)序reads接頭序列、截短后長(zhǎng)度小于100 bp的序列、N含量大于5%的序列后統(tǒng)計(jì)GC含量.采取所有樣本混合組裝的策略，將所有樣本進(jìn)行歸一化得到Unigene(即各個(gè)樣本所有基因進(jìn)行去冗余去重復(fù)后得到的唯一序列)，然后使用長(zhǎng)度、GC含量以及N50作為標(biāo)準(zhǔn)對(duì)Unigene組裝質(zhì)量進(jìn)行評(píng)判.N50是將所有組裝結(jié)果按照長(zhǎng)度從高到低排列后，達(dá)到總長(zhǎng)度一半時(shí)的長(zhǎng)度.選取權(quán)威的NCBI_nr、GO、KEGG、Pfam、Swiss-Prot和eggNOG 6種數(shù)據(jù)庫(kù)，使用軟件DIAMOND添加相應(yīng)的功能注釋?zhuān)偻ㄟ^(guò)RSEM軟件計(jì)算各基因的表達(dá)水平，然后通過(guò)R語(yǔ)言對(duì)差異基因進(jìn)行表達(dá)分析.

2.2 三葉木通Unigene的GO和KEGG功能注釋

2.2.1 三葉木通Unigenes的GO功能分析

對(duì)三葉木通的Unigenes進(jìn)行GO功能富集.共注釋到75 631種不同的GO功能，其中生物學(xué)過(guò)程(Biological Process)占比最多，涉及到29 332個(gè)不同的功能，而細(xì)胞組分(Cellular Component)共有個(gè)26090個(gè)功能被注釋到，富集到分子功能(molecular_function)共有23 929個(gè)，同一個(gè)Unigene可以對(duì)應(yīng)多種不同的功能，按照功能注釋富集的Unigenes如圖3所示.

注：1～25：生物過(guò)程；26～40：細(xì)胞組分；41～50：分子功能.1.生物過(guò)程；2.DNA轉(zhuǎn)錄調(diào)控；3.DNA轉(zhuǎn)錄；4.氧化還原過(guò)程；5.蛋白質(zhì)磷酸化；6.防御反應(yīng)；7.對(duì)鹽脅迫的反應(yīng)；8.對(duì)脫落酸的反應(yīng)；9.轉(zhuǎn)化；10.胚發(fā)育結(jié)束于種子休眠；11.對(duì)鎘離子的反應(yīng)；12.對(duì)冷脅迫的反應(yīng)；13.蛋白質(zhì)泛素化；14.對(duì)水分缺乏的反應(yīng)；15.蛋白質(zhì)折疊；16.蛋白質(zhì)水解；17.對(duì)細(xì)菌的防御反應(yīng)；18.磷酸化；19.蛋白酶體介導(dǎo)的泛素依賴(lài)性蛋白質(zhì)分解代謝過(guò)程；20.多細(xì)胞生物發(fā)育；21.脫落酸激活的信號(hào)通路；22.信號(hào)轉(zhuǎn)導(dǎo)；23.細(xì)胞內(nèi)蛋白質(zhì)轉(zhuǎn)運(yùn)；24.蛋白質(zhì)轉(zhuǎn)運(yùn)；25.細(xì)胞內(nèi)信號(hào)轉(zhuǎn)導(dǎo)；26.細(xì)胞核；27.細(xì)胞質(zhì)；28.質(zhì)膜；29.膜的組成成分；30.葉綠體；31.細(xì)胞質(zhì)；32.線粒體；33.膜；34.高爾基體；35.胞外區(qū)；36.胞間連絲；37.內(nèi)質(zhì)網(wǎng)；38.葉綠體基質(zhì)；39.葉綠體包膜；40.液泡；41.分子功能；42.蛋白質(zhì)結(jié)合；43.ATP結(jié)合；44.金屬離子結(jié)合；45.DNA結(jié)合；46.DNA結(jié)合轉(zhuǎn)錄因子活性；47.RNA結(jié)合；48.蛋白質(zhì)絲氨酸/蘇氨酸激酶活性；49.激酶活性；50.鋅離子結(jié)合.

2.2.2 三葉木通的KEGG功能分析

對(duì)三葉木通的雌雄花芽進(jìn)行KEGG pathway聚類(lèi)分析，共注釋到10 017個(gè)Unigenes，聚焦到20個(gè)生物途徑上.其中顯著富集的通路(如圖4所示)主要是代謝途徑和遺傳信息處理，具體為碳水化合物，脂質(zhì)，氨基酸的代謝，蛋白質(zhì)折疊、翻譯、分類(lèi)與降解，淀粉和蔗糖代謝等.在對(duì)通路基因富集分析后發(fā)現(xiàn)，許多基因與生長(zhǎng)素、細(xì)胞分裂素相關(guān)的通路密切聯(lián)系.

注：1：有機(jī)體系；2～11：代謝；12～13：人類(lèi)疾病；14～17：遺傳信息處理；18～19：環(huán)境信息處理；20：細(xì)胞過(guò)程.1.環(huán)境適應(yīng)性；2.輔助因子和維生素的代謝；3.其他氨基酸的代謝；4.核苷酸的代謝；5.其他次生代謝產(chǎn)物的生物合成；6.氨基酸的代謝；7.糖類(lèi)的生物合成和代謝；8.脂質(zhì)的代謝；9.碳水化合物的代謝；10.能量代謝；11.萜類(lèi)和聚酮化合物的代謝；12.內(nèi)分泌和代謝疾??；13.抗藥性；14.復(fù)制和修復(fù)；15.轉(zhuǎn)錄；16.折疊、分類(lèi)和降解；17.翻譯；18.膜轉(zhuǎn)運(yùn)；19.信號(hào)傳導(dǎo)；20.轉(zhuǎn)運(yùn)和分解代謝.

2.2.3 三葉木通雌雄花芽差異基因富集分析

對(duì)轉(zhuǎn)錄本進(jìn)行篩選分析后，以雌花為對(duì)照組，發(fā)現(xiàn)雄花中一共有2 579個(gè)Unigenes表達(dá)上調(diào)，2 793個(gè)Unigenes表達(dá)下調(diào).然后利用R語(yǔ)言進(jìn)一步對(duì)這些差異基因繪制火山圖(如圖5所示).通過(guò)差異基因聚類(lèi)分析，可以判斷基因在不同實(shí)驗(yàn)條件下調(diào)控模式的聚類(lèi)模式，根據(jù)樣品基因表達(dá)譜的相近程度繪制了差異基因聚類(lèi)熱圖(如圖6所示)，為了更好地直觀反映聚類(lèi)表達(dá)模式，將差異基因TPM通過(guò)Z值方式進(jìn)行基因表達(dá)展示.其中橫坐標(biāo)為樣本，縱坐標(biāo)為基因，不同的顏色表示不同的基因表達(dá)水平.

圖5 三葉木通差異基因火山圖

圖6 三葉木通差異基因熱圖

2.2.4 三葉木通雌雄花芽差異基因的q-PCR驗(yàn)證

通過(guò)對(duì)差異基因進(jìn)行篩選，選取了10個(gè)基因進(jìn)行 了qRT-PCR驗(yàn)證，這些基因均在雄花芽中產(chǎn)生了高表達(dá)，與測(cè)序結(jié)果一致，說(shuō)明轉(zhuǎn)錄組數(shù)據(jù)真實(shí)可靠.同時(shí)，對(duì)這些差異基因進(jìn)一步挖掘分析表明，與花發(fā)育相關(guān)的MADS-box家族的基因與Unigenes的部分結(jié)果具有較高的重合度，從NCBI上已公布的三葉木通AG、SEP3、AP3序列進(jìn)行了克隆并同時(shí)進(jìn)行了序列比對(duì)，結(jié)果表明堿基重合度很高(如圖7所示).

注：**代表差異性極顯著，P<0.01

2.2.5 三葉木通的差異基因富集性分析

將三葉木通雌花與雄花的差異表達(dá)基因與GO以及KEGG數(shù)據(jù)庫(kù)比對(duì)分析，結(jié)果與Unigenes差異主要在細(xì)胞組成部分.而差異基因KEGG功能分析(kegg通路 map04626，map04075)顯示，共注釋3160個(gè)差異基因，映射到130條代謝通路，其中植物激素信號(hào)轉(zhuǎn)導(dǎo)與植物病原互作包含的差異基因數(shù)量最多，分別達(dá)到118和113條，而顯著富集的通路還包括淀粉和蔗糖代謝，糖酵解以及苯丙烷生物合成.這些通路在三葉木通性別分化發(fā)揮了重要作用.

2.2.6 三葉木通性別相關(guān)基因q-PCR驗(yàn)證

三葉木通的差異基因聚類(lèi)分析以及GO富集分析均表明，Unigenes主要映射在細(xì)胞分裂素和丙酮酸代謝相關(guān)通路上，因此隨機(jī)選擇9個(gè)與此相關(guān)的差異基因進(jìn)行q-PCR驗(yàn)證.結(jié)果顯示這些基因均在雄花中大量表達(dá)，差異性極強(qiáng)(如圖8～9所示)t檢驗(yàn)表明具有統(tǒng)計(jì)學(xué)意義，表明這些基因可能在雄花的發(fā)育過(guò)程中起著重要作用.

注：**代表差異性極顯著，P<0.01

注：**代表差異性極顯著，P<0.01

3 討 論

植物性別分化是植物生長(zhǎng)發(fā)育過(guò)程中由性別決定基因、環(huán)境因素和激素共同調(diào)控的.現(xiàn)有的研究表明，在部分植物中存在性染色體.日本京都大學(xué)課題組在美味獼猴桃中發(fā)現(xiàn)Y染色體上的特異性基因ShyGirl能夠抑制心皮的發(fā)育[16]，而另外一種特異性基因FriendlyBoy則能夠維持雄花的發(fā)育[17].在黃瓜性別的研究中發(fā)現(xiàn)，自然條件下比人工培育條件下的雌花更少，推測(cè)可能是黃瓜體內(nèi)的miRNA受到環(huán)境的影響選擇性的使心皮敗育，同時(shí)在人工培育條件下施加外源激素乙烯能夠有效促使黃瓜雌花開(kāi)放[18].在對(duì)楊樹(shù)染色體分析后，科學(xué)家發(fā)現(xiàn)XIX染色體上出現(xiàn)了一部分Y染色體連鎖的能夠產(chǎn)生siRNA阻斷對(duì)應(yīng)基因的表達(dá)，從而抑制雌花生成[19].目前，對(duì)植物性別決定的研究主要還處于利用分子標(biāo)記尋找特定的性別遺傳基因，對(duì)于性染色體的劃分以及具體作用還有一定的爭(zhēng)議，雖然近年來(lái)分離鑒定了許多microRNA，發(fā)現(xiàn)甲基化能夠?qū)χ参锏男詣e決定產(chǎn)生影響，但性別分化的機(jī)制仍不清晰.

三葉木通的差異基因GO分析表明，三葉木通的雌雄花芽的差異基因主要在細(xì)胞組成上，與雌雄花結(jié)構(gòu)差異相關(guān).差異基因的KEGG富集分析表明，三葉木通雌雄花芽的差異基因主要富集于植物信號(hào)轉(zhuǎn)導(dǎo)，結(jié)合前人的分析研究推斷，植物激素在三葉木通性別分化的過(guò)程中起了關(guān)鍵的調(diào)節(jié)作用.而Unigenes的聚類(lèi)分析表明丙酮酸、細(xì)胞分裂素和油菜素內(nèi)酯在性別分化過(guò)程中扮演了重要的角色，參與了多種代謝過(guò)程(kegg通路map04075)，細(xì)胞分裂素能夠與其他激素綜合調(diào)控性別分化.油菜素內(nèi)酯在植物體內(nèi)廣泛存在，通過(guò)對(duì)油菜素內(nèi)酯物的調(diào)節(jié)，能夠影響花器官形成的生理過(guò)程，進(jìn)而影響性別分化[20].丙酮酸參與的代謝調(diào)控網(wǎng)絡(luò)非常復(fù)雜(如圖10所示)，其中在生物體內(nèi)的物質(zhì)代謝中起著關(guān)鍵橋梁作用，在糖酵解的途徑中氧化生成乙酰CoA，參與三羧酸循環(huán)，從而實(shí)現(xiàn)生物體內(nèi)三大營(yíng)養(yǎng)物質(zhì)的互相轉(zhuǎn)化.而丙酮酸代謝出現(xiàn)問(wèn)題，可能會(huì)影響絨氈層發(fā)育以及減數(shù)分裂，進(jìn)而對(duì)生殖過(guò)程以及性別分化產(chǎn)生影響[21].

圖10 丙酮酸部分代謝通路

GO分析以及qRT-PCR的驗(yàn)證發(fā)現(xiàn)，三葉木通雌雄花芽的差異基因囊括了幾乎已公布的三葉木通MADS-BOX基因家族中的B類(lèi)以及C類(lèi)基因，這與花發(fā)育模型理論中，B類(lèi)以及C類(lèi)基因決定雄蕊和心皮的發(fā)育[22]是一致的，推測(cè)在植物性別分化的過(guò)程中，MADS-box家族中的基因除了直接參與花器官形態(tài)建成之外，可能也在早期參與了性別分化調(diào)控的過(guò)程.pathway的分析表明，激素在雌雄花芽的發(fā)育也起到了不可或缺的關(guān)鍵性作用，通過(guò)進(jìn)一步檢索發(fā)現(xiàn)，在激素相關(guān)的通路中，差異基因更多的富集在細(xì)胞分裂素相關(guān)的通路上，油菜素內(nèi)酯和丙酮酸也發(fā)揮了重要調(diào)控作用.已有研究現(xiàn)性別分化是一個(gè)由多個(gè)因素共同調(diào)控決定的復(fù)雜的生理過(guò)程，三葉木通的性別分化機(jī)制，還有待于進(jìn)一步研究.

4 結(jié) 論

本實(shí)驗(yàn)通過(guò)轉(zhuǎn)錄組測(cè)序組裝，在雌雄花芽各3個(gè)轉(zhuǎn)錄本中一共獲得了17 439個(gè)Unigenes，平均長(zhǎng)度860 bp，N50為1 471，Q20大于98%，Q30大于93%，大部分Unigenes長(zhǎng)度位介于200～2 000 bp之間，表明三葉木通轉(zhuǎn)錄組序列組裝質(zhì)量良好.通過(guò)實(shí)時(shí)熒光定量PCR驗(yàn)證，結(jié)果表明轉(zhuǎn)錄組數(shù)據(jù)真實(shí)可靠.將組裝好的通過(guò)對(duì)比公共數(shù)據(jù)庫(kù)發(fā)現(xiàn)，三葉木通在NR數(shù)據(jù)庫(kù)比對(duì)的序列數(shù)量最多有14 327條，但在KEGG的注釋只有10 017條，比例僅為57.44%.差異基因以及綜合數(shù)據(jù)庫(kù)的對(duì)比分析表明，三葉木通已公布的MADS-box基因均屬于差異基因，同時(shí)KEGG聚類(lèi)分析表明細(xì)胞分裂素和油菜素內(nèi)酯的信號(hào)通路是差異基因富集較多的通路，在雄性花芽的發(fā)育過(guò)程中起著正調(diào)控作用，為后續(xù)對(duì)三葉木通的性別分化機(jī)制研究奠定了基礎(chǔ)并提供了方向.