三葉木通雌雄花芽轉錄組測序及分析

田 天 項 明 陳發菊

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

三葉木通(Akebiatrifoliate)是木通科(Lardizabalaceae)木通屬(Akebia Decne.)落葉木質藤本植物，民間又稱八月炸、八月瓜.三葉木通花序為總狀花序，下部生有1～2朵雌花，上部開有15～30朵雄花.三葉木通早期均為兩性花結構，在發育的后期雌蕊和雄蕊會選擇性退化而形成單性花[1].

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

本實驗通過三葉木通的雌雄花芽進行轉錄組測序，利用現有的數據庫與測序結果進行比對，并對雌雄花芽的差異基因進行驗證分析，為探究三葉木通性別分化機制奠定基礎及利用基因工程培育新品種提供方向.

1 材料與方法

1.1 實驗材料

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

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

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

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

1.3 質量控制、組裝統計、功能注釋與表達分析

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

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

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

1.5 三葉木通雌雄花芽差異表達基因的qRT-PCR驗證

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

表1 實時熒光PCR引物

2 結果與分析

2.1 三葉木通轉錄組數據測序及序列組裝

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

圖1 Unigene組裝長度統計 圖2 Unigene的GC含量統計

表2 過濾前后的Reads統計

表3 Unigene功能注釋統計

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

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

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

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

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

2.2.2 三葉木通的KEGG功能分析

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

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

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

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

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

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

2.2.4 三葉木通雌雄花芽差異基因的q-PCR驗證

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

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

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

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

2.2.6 三葉木通性別相關基因q-PCR驗證

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

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

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

3 討 論

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

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

圖10 丙酮酸部分代謝通路

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

4 結 論

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