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

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

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

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

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

1 材料與方法

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

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

1.2 三葉木通總RNA提取、文庫制備及測序

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

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

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

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

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

根據(jù)Patro[12]的方法計(jì)算TPM[13]衡量所有Unigenes表達(dá)水平，利用package edgeR[14]對Unigenes進(jìn)行差異表達(dá)分析，然后通過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)錄組測序結(jié)果，隨機(jī)篩選10條差異表達(dá)明顯的基因，進(jìn)一步分析篩選得到細(xì)胞分裂素以及丙酮酸相關(guān)的基因9個.通過實(shí)時熒光定量PCR儀(CFX96，BIO-RAD)以及SYBR qPCR試劑盒(南京諾唯贊生物科技有限公司)進(jìn)行qRT-PCR實(shí)驗(yàn).通過NCBI BLAST設(shè)計(jì)qPCR引物，由華大基因股份有限公司合成.qRT-PCR引物見表1，產(chǎn)物長度100～250 bp.選擇RG5(EF-α)作為內(nèi)參基因[15]，上游引物為GTCTCCCTCTTCAGGATGTTTAC，下游引物為GACAACCATACCAGGCTTGA.其中前10對引物用來驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)的可靠性，后9對引物為代謝通路引物.

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

2 結(jié)果與分析

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

采用Illumina測序平臺技術(shù)對三葉木通樣本進(jìn)行轉(zhuǎn)錄組測序，共得到49.84 Gb的原始數(shù)據(jù)(raw data).過濾掉不合格的序列，共得到49.06 Gb有效數(shù)據(jù)，占原始數(shù)據(jù)的98.43%，組裝后GC含量與長度分布如圖1～圖2所示，每組數(shù)據(jù)質(zhì)量見表2.對過濾后的數(shù)據(jù)進(jìn)行Trinity組裝，去除冗余后得到17 439個Unigene，通過六大數(shù)據(jù)庫注釋，所有基因均被一個或者更多的數(shù)據(jù)庫注釋到(見表3).

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

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

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

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

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

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

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

注：1～25：生物過程；26～40：細(xì)胞組分；41～50：分子功能.1.生物過程；2.DNA轉(zhuǎn)錄調(diào)控；3.DNA轉(zhuǎn)錄；4.氧化還原過程；5.蛋白質(zhì)磷酸化；6.防御反應(yīng)；7.對鹽脅迫的反應(yīng)；8.對脫落酸的反應(yīng)；9.轉(zhuǎn)化；10.胚發(fā)育結(jié)束于種子休眠；11.對鎘離子的反應(yīng)；12.對冷脅迫的反應(yīng)；13.蛋白質(zhì)泛素化；14.對水分缺乏的反應(yīng)；15.蛋白質(zhì)折疊；16.蛋白質(zhì)水解；17.對細(xì)菌的防御反應(yīng)；18.磷酸化；19.蛋白酶體介導(dǎo)的泛素依賴性蛋白質(zhì)分解代謝過程；20.多細(xì)胞生物發(fā)育；21.脫落酸激活的信號通路；22.信號轉(zhuǎn)導(dǎo)；23.細(xì)胞內(nèi)蛋白質(zhì)轉(zhuǎn)運(yùn)；24.蛋白質(zhì)轉(zhuǎn)運(yùn)；25.細(xì)胞內(nèi)信號轉(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功能分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討 論

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

三葉木通的差異基因GO分析表明，三葉木通的雌雄花芽的差異基因主要在細(xì)胞組成上，與雌雄花結(jié)構(gòu)差異相關(guān).差異基因的KEGG富集分析表明，三葉木通雌雄花芽的差異基因主要富集于植物信號轉(zhuǎn)導(dǎo)，結(jié)合前人的分析研究推斷，植物激素在三葉木通性別分化的過程中起了關(guān)鍵的調(diào)節(jié)作用.而Unigenes的聚類分析表明丙酮酸、細(xì)胞分裂素和油菜素內(nèi)酯在性別分化過程中扮演了重要的角色，參與了多種代謝過程(kegg通路map04075)，細(xì)胞分裂素能夠與其他激素綜合調(diào)控性別分化.油菜素內(nèi)酯在植物體內(nèi)廣泛存在，通過對油菜素內(nèi)酯物的調(diào)節(jié)，能夠影響花器官形成的生理過程，進(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)物質(zhì)的互相轉(zhuǎn)化.而丙酮酸代謝出現(xiàn)問題，可能會影響絨氈層發(fā)育以及減數(shù)分裂，進(jìn)而對生殖過程以及性別分化產(chǎn)生影響[21].

圖10 丙酮酸部分代謝通路

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

4 結(jié) 論

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