姜敏杰， 畢桂萁， 王津果， 張敬宇， 胡依依， 李曉東， 劉 娟， 隋正紅

(中國海洋大學(xué)海洋生物遺傳學(xué)與育種教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266003)

紅藻龍須菜(Gracilariopsislemaneiformis)可提取瓊膠，也可作為鮑魚餌料[1-3]，已經(jīng)成為中國僅次于海帶的第二大栽培海藻[4]。另外龍須菜兼具吸收氮磷，緩解富營(yíng)養(yǎng)化的修復(fù)功能[5]。目前共有3個(gè)中國認(rèn)定的良種：“981”[6]，“2007”[7]和“魯龍1號(hào)”[8]。對(duì)不同良種生理性狀的研究發(fā)現(xiàn)“981”具備與其他良種不同的特性[6,9-12]，其在四分孢子體成熟時(shí)放散孢子數(shù)量低，孢子畸形，即表現(xiàn)出低育的特性[10-12]。此前對(duì)龍須菜的繁殖缺乏分子層面的研究，981品系的低育特性為龍須菜的育性相關(guān)研究提供了良好材料。另外目前龍須菜的栽培產(chǎn)業(yè)仍以營(yíng)養(yǎng)繁殖為主，低育性狀也具有更高的產(chǎn)業(yè)價(jià)值。

近年來測(cè)序技術(shù)發(fā)展和成本降低使得組學(xué)分析手段更加多樣化，轉(zhuǎn)錄組和表達(dá)譜數(shù)據(jù)也成為常見的原始分析數(shù)據(jù)之一。從高通量基因表達(dá)數(shù)據(jù)中篩選挖掘出所關(guān)注的基因信息成為轉(zhuǎn)錄組分析中必要的一環(huán)。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted Gene Co-Expression Network Analysis, WGCNA)側(cè)重于基因之間的相關(guān)性以及基因表達(dá)與所研究性狀之間的關(guān)聯(lián)性[13]，是一種高效、準(zhǔn)確的生物信息學(xué)分析技術(shù)。2005 年，Zhang 和Horvath[14]最早提出WGCNA分析理論。2008 年，Langfelder和Horvath[15]在該研究方法的基礎(chǔ)上開發(fā)了WGCNA的R軟件包，其分析思路是定義基因共表達(dá)相關(guān)矩陣并確定基因網(wǎng)絡(luò)形成的鄰接函數(shù)，進(jìn)而計(jì)算不同節(jié)點(diǎn)的相異系數(shù)并形成相異矩陣，根據(jù)相異矩陣進(jìn)行層次聚類后依此將基因劃分成以樞紐基因(Hub gene)為核心的基因共表達(dá)模塊(Module)，進(jìn)一步將模塊與性狀特征相關(guān)聯(lián)從而提取出目的信息[13]。這一技術(shù)被廣泛應(yīng)用于動(dòng)植物性狀的調(diào)控機(jī)制解析上，如對(duì)玉米花期基因[16]、胡楊異形葉發(fā)育[17]、蘋果黃皮突變[18]、番茄果實(shí)成熟[19]以及雨生紅球藻中蝦青素的積累[20]等過程的相關(guān)分析中。

針對(duì)龍須菜的各類組學(xué)研究也在近幾年相繼展開。Zhou等[21]通過HiSeq 2000測(cè)序預(yù)測(cè)龍須菜基因組大小約為97 Mb，基因數(shù)約為3 490個(gè)以及GC%含量為48%。Hu等[22]在基因組survey結(jié)果的基礎(chǔ)上開發(fā)了豐富的簡(jiǎn)單序列重復(fù)(SSR)標(biāo)記，并用多態(tài)性SSR標(biāo)記對(duì)龍須菜不同地理群體進(jìn)行遺傳多樣性分析。Huang等[23]對(duì)龍須菜野生型和綠色突變型的轉(zhuǎn)錄組數(shù)據(jù)比對(duì)分析后認(rèn)為藻紅蛋白、藻紅膽素和葉綠素光捕捉復(fù)合體的表達(dá)水平改變導(dǎo)致突變的形成。Sun等[24]在對(duì)981品系的全基因組分析中首次報(bào)道了與植物激素如生長(zhǎng)素，水楊酸和茉莉酮酸等相關(guān)的信號(hào)通路。在之前針對(duì)龍須菜四分孢子體發(fā)育放散時(shí)期的生理研究中，認(rèn)為四分孢子體的發(fā)育時(shí)期可劃分成未成熟期(I期)即四分孢子囊未形成時(shí)期、成熟前期(II期)即四分孢子囊開始形成但未有四分孢子放散時(shí)期、成熟期(III)即四分孢子大量放散時(shí)期和成熟后期(IV期)即四分孢子基本放散留下大量孔洞時(shí)期[25]。對(duì)981品系的生理研究發(fā)現(xiàn)其在成熟期(III期)釋放的孢子數(shù)量少且形態(tài)畸形[12]，根據(jù)對(duì)四分孢子形成過程的認(rèn)識(shí)推測(cè)981品系可能在II期經(jīng)減數(shù)分裂形成四分孢子這一過程中出現(xiàn)異常。本文以I、II、IV期不同時(shí)期的藻體為材料，進(jìn)行表達(dá)譜分析后結(jié)合WGCNA技術(shù)，探索龍須菜四分孢子體發(fā)育過程中的相關(guān)調(diào)控機(jī)制，并對(duì)981與魯龍1號(hào)品系之間存在的差異進(jìn)行基因表達(dá)層面上的探究。

1 材料與方法

1.1 材料培養(yǎng)

981品系和魯龍1號(hào)品系于2014 年5 月采自青島膠州灣養(yǎng)殖海區(qū)(36°08′N，120°17′E)。使用毛刷盡可能去除藻體表面污染，并用滅菌海水沖洗3遍。每個(gè)品系的一級(jí)分枝切成2 cm小藻段，在含有Pro 培養(yǎng)基[26]的200 mL無菌海水的250 mL三角瓶中放置5 段，分別將981品系和魯龍1號(hào)品系置于其四分孢子放散的最佳條件下培養(yǎng)[27]，每3 d更換一次新鮮培養(yǎng)基。每周在光鏡下對(duì)藻段進(jìn)行觀察以確定四分孢子的形成和放散情況。樣品收取的具體信息見表1，每個(gè)階段樣品設(shè)置3個(gè)生物學(xué)平行。

1.2 總RNA提取與文庫構(gòu)建

采用總RNA 提取試劑盒(天根生化科技有限公司，北京)提取RNA，使用RNase-Free DNase I(天根生化科技有限公司，北京)進(jìn)行柱上DNA消化。提取RNA后進(jìn)行瓊脂糖凝膠電泳檢測(cè)，并用Nanodrop 2000和Agilent 2100檢測(cè)濃度和純度。

使用NEBNext Ultra RNA Library Prep Kit for Illumina(NEB，USA)試劑盒按照說明書步驟進(jìn)行18個(gè)RNA-seq文庫的構(gòu)建，建庫完成后進(jìn)行IlluminaHiSeq 2000 測(cè)序。文庫的構(gòu)建以及測(cè)序均由北京諾禾致源生物信息科技有限公司完成。

1.3 測(cè)序數(shù)據(jù)質(zhì)控與分析

對(duì)原始測(cè)序序列(Raw reads)進(jìn)行質(zhì)量控制，過濾掉帶接頭以及低質(zhì)量的序列得到clean reads。使用Hisat2[28]軟件將clean reads與龍須菜Survey數(shù)據(jù)[21]比對(duì)，將比對(duì)成功的reads挑出后用SOAPdenovo-Trans進(jìn)行拼接構(gòu)建參考序列(Reference)，拼接后過濾掉小于200 bp的contig。使用bowtie對(duì)每個(gè)樣品的reads進(jìn)行比對(duì)并使用RSEM統(tǒng)計(jì)樣品中基因的read count數(shù)目，edgeR[29-30]進(jìn)行差異基因分析，TMM算法對(duì)read count數(shù)換算成FPKM (fragments per kilobase of exon per million fragments mapped, 每千堿基外顯子百萬片段數(shù))值來估算基因表達(dá)水平并進(jìn)行矩陣標(biāo)準(zhǔn)化矯正。對(duì)差異表達(dá)基因的表達(dá)矩陣進(jìn)行提取，聚類檢驗(yàn)樣品重復(fù)性。

表1 龍須菜表達(dá)譜分析的樣品信息

1.4 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

利用R 軟件的WGCNA軟件包完成共表達(dá)網(wǎng)絡(luò)構(gòu)建。首先, 導(dǎo)入預(yù)處理的表達(dá)量表對(duì)數(shù)據(jù)進(jìn)行過濾后構(gòu)建表達(dá)矩陣，并對(duì)樣品進(jìn)行聚類分析；利用WGCNA 中的pickSoftThreshold函數(shù)計(jì)算軟閾值β參數(shù)使得滿足無尺度網(wǎng)絡(luò)分布前提條件；選取R2達(dá)到0.9 時(shí)的閾值參數(shù)β[13]，計(jì)算所有基因間的相異系數(shù)創(chuàng)建鄰接矩陣，并用TOMsimilarity函數(shù)將鄰接矩陣轉(zhuǎn)換為拓?fù)渚仃嚭?，取逆得到相異性矩陣并以此?gòu)建系統(tǒng)聚類樹。設(shè)定每個(gè)模塊最小基因數(shù)為30 個(gè)，利用cutreeDynamic動(dòng)態(tài)剪枝法修剪聚類樹識(shí)別初始共表達(dá)模塊。隨后用moduleEigengene函數(shù)計(jì)算每個(gè)模塊的特征向量值ME (moduleEigengene)，用mergeCloseModules函數(shù)選擇相似性在100%以上的模塊進(jìn)行合并融合(cutHeight=0.0)。

1.5 特異性模塊的鑒定及樞紐基因的篩選

重新計(jì)算各模塊特征向量基因ME，利用cor函數(shù)量化模塊ME與樣本特征的相關(guān)系數(shù)并繪制熱圖。相關(guān)系數(shù)越高表示模塊內(nèi)基因與樣本特征關(guān)聯(lián)度越高，相應(yīng)地，基因相對(duì)表達(dá)量越高。選擇與樣本特征具備較高相關(guān)系數(shù)的模塊作為備選特異性模塊。

利用Alldegrees1函數(shù)計(jì)算模塊內(nèi)基因間連接數(shù)，篩選模塊內(nèi)連接數(shù)前10%的基因?yàn)镠ub基因。

1.6 模塊的GO(Gene Ontology)富集分析和基因互作網(wǎng)絡(luò)的構(gòu)建

對(duì)模塊內(nèi)基因使用Cytoscape軟件中的binGO功能進(jìn)行GO富集分析，并用Cytoscape[31]軟件輸出網(wǎng)絡(luò)圖像?；贙EGG 生物學(xué)通路數(shù)據(jù)庫(http://www.kegg.jp/)，進(jìn)行KEGG通路富集分析，并借助http://www.omicshare.com的在線工具作圖。

2 結(jié)果與分析

2.1 WGCNA數(shù)據(jù)的預(yù)處理及去除離群樣本

為了排除污染和盡可能的利用基因組信息，對(duì)測(cè)序數(shù)據(jù)向基因組survey拼接結(jié)果上進(jìn)行RNA-seq的可變剪切mapping，將有能map到的RNA-seq數(shù)據(jù)進(jìn)行分離。然后，利用分離出的981品系的RNA-seq reads進(jìn)行de novo拼接，構(gòu)建reference。拼接后，得到13 426條Unigene。并對(duì)轉(zhuǎn)錄組reference分別進(jìn)行GO和KEGG代謝通路注釋(見圖1)。KEGG結(jié)果中注釋到一級(jí)大類為新陳代謝(Metabolism)、遺傳信息過程(Genetic information processing)、環(huán)境信息過程(Environmental information processing)、細(xì)胞過程(Celluar processing)、生物系統(tǒng)(Organismal systems)。注釋到的基因數(shù)最多的二級(jí)分類為翻譯(Translation)過程，有224個(gè)基因(見圖1A)。GO結(jié)果中多數(shù)基因富集到代謝過程(Metabolic process)、細(xì)胞過程(Cellular process)、細(xì)胞(Cell)、細(xì)胞部分(Cell part)、催化活性(Catalytic activity)、結(jié)合(Binding)等Go term(見圖1B)。后續(xù)樣本定量分析與功能注釋則采用此reference。

對(duì)18個(gè)樣品的差異表達(dá)基因的表達(dá)矩陣進(jìn)行提取后聚類。初步聚類分析后，為確保分析結(jié)果的嚴(yán)謹(jǐn)和可靠性，確定采納具有較高相關(guān)性的樣品，去除皮爾遜相關(guān)系數(shù)小于0.95的樣品，每個(gè)取樣階段保留2個(gè)生物學(xué)平行。對(duì)樣品之間相關(guān)性進(jìn)行統(tǒng)計(jì)見圖2，除981.IV組樣品，相關(guān)性為0.89。其他樣品相關(guān)性分別在0.95～0.98之間，顯示相關(guān)性較好。樣品間的聚類分析，顯示重復(fù)性無問題。

圖1 龍須菜轉(zhuǎn)錄本KEGG代謝通路注釋(A)與GO注釋(B)二級(jí)分類統(tǒng)計(jì)

將構(gòu)建好的樣品FPKM表達(dá)矩陣導(dǎo)入WGCNA軟件包，進(jìn)行g(shù)sg數(shù)據(jù)過濾，共去除857個(gè)數(shù)據(jù)不良的基因，剩余12 569個(gè)基因進(jìn)行后續(xù)分析。

2.2 基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

設(shè)置軟閾值β為7，計(jì)算所有基因間的相異系數(shù)創(chuàng)建鄰接矩陣和相異性矩陣并以此構(gòu)建層次聚類樹(見圖3A)，每個(gè)顏色代表一個(gè)模塊，共得到32個(gè)表達(dá)模塊。構(gòu)建每個(gè)模塊與樣本特征的相關(guān)性熱圖和模塊內(nèi)基因數(shù)量的條形圖，熱圖中橫軸為樣品特征，縱軸為32個(gè)模塊，熱圖中顏色越深表示模塊與特征值的相關(guān)性越高(見圖3B)。最大的模塊為Turquoise模塊，包含2 577個(gè)基因，最小基因模塊為Violet模塊，有66個(gè)基因(見圖3B)。

經(jīng)篩選發(fā)現(xiàn)與品系(Cultivar)相關(guān)的模塊有 6個(gè)(圖3B中相關(guān)性系數(shù)絕對(duì)值高于0.7的模塊)： Yellow模塊(cor=-0.75，P=0.005)、Magenta模塊(cor=0.74，P=0.006)、Greenyellow模塊(cor=0.87，P=3e-4)、Cyan模塊(cor=-0.75，P=0.005)、Darkgreen模塊(cor=-0.79，P=0.002)五個(gè)模塊，以及一個(gè)備選模塊 Purple模塊(cor=0.68，P=0.02)。相關(guān)性為正值，表示基因在981品系中呈現(xiàn)較高表達(dá)水平(與魯龍1號(hào)相比)，相關(guān)性為負(fù)值，表示基因在981品系中呈現(xiàn)較低表達(dá)水平(與魯龍1號(hào)相比)；P<0.5顯示相關(guān)性分析可靠。

另一方面，由于981品系表現(xiàn)出四分孢子放散時(shí)期的低育特性，因此選擇四分孢子形成的時(shí)期即981.II期作為關(guān)注對(duì)象，發(fā)現(xiàn)與此特征相關(guān)的模塊(以0.8為閾值)有：Magenta模塊(cor=0.82，P=0.001)、Green模塊(cor=0.97，P=1e-07)、Turquoise模塊(cor=0.99，P=5e-10)3個(gè)模塊。相關(guān)性的正負(fù)表示表達(dá)水平的高低。

2.3 品系相關(guān)模塊的分析與富集

對(duì)與品系(Cultivar)相關(guān)的6個(gè)模塊的基因及模塊特征基因進(jìn)行表達(dá)量分析，構(gòu)建樣品的表達(dá)熱圖(見圖4)。熱圖和柱狀圖為相應(yīng)顏色模塊基因表達(dá)量熱圖，每一行基因都進(jìn)行了相應(yīng)的均一化，顏色越紅表明基因相對(duì)表達(dá)量越高，顏色越綠則越低；柱狀圖為表達(dá)量均一化后的相對(duì)表達(dá)量，可以代表每一列，即每一個(gè)樣品的表達(dá)情況。其中Yellow模塊、Cyan模塊和Darkgreen模塊在981品系的I～I(xiàn)V期樣品中均下調(diào)表達(dá)，而Magenta模塊、Greenyellow模塊和Purple模塊的基因在981品系中均有不同程度的上調(diào)表達(dá)。

進(jìn)一步地將6個(gè)模塊的基因混合，通過FPKM數(shù)據(jù)處理分為981中高表達(dá)基因集與魯龍1號(hào)品系中高表達(dá)基因集。顯示在981品系中，在I～I(xiàn)V期樣品中均保持高表達(dá)水平的基因?yàn)?15個(gè)，設(shè)為Up組；相反，表達(dá)水平比魯龍1號(hào)品系低的基因有817個(gè)，設(shè)為Down組。分別對(duì)Up與Down組使用binGO進(jìn)行 GO富集分析。

圖5中有顏色的圓圈表示有基因富集，顏色越深顯著性越高，圓圈越大相應(yīng)富集的基因越多。981中高表達(dá)的基因中(Up組)，主要富集到信號(hào)轉(zhuǎn)導(dǎo)(Signal transduction )(GO: 0007165, FDR= 1.39e-01)、轉(zhuǎn)運(yùn)(Transport) (GO: 0006810, FDR= 1.39e-01)、細(xì)胞分化(Cell differentiation)(GO: 0030154, FDR= 9.66e-02)等終端GO term上，表明981品系在四分孢子形成、成熟至放散過程中，信號(hào)調(diào)節(jié)、細(xì)胞轉(zhuǎn)運(yùn)、細(xì)胞分化等生物學(xué)過程比較活躍。

981品系低表達(dá)組(Down組)主要富集到翻譯(Translation)(GO: 0006412, FDR=6.82e-09)、細(xì)胞穩(wěn)態(tài)(Cellular homeostasis)(GO: 0019725, FDR=7.07e-03)、生物合成過程(Biosynthetic process)(GO: 0009058, FDR= 3.82e-02)等終端GO term上，表明981品系在四分孢子成熟至放散過程中生物合成、細(xì)胞穩(wěn)態(tài)等受到抑制。同樣，核糖體(Ribosome)GO: 0005840, FDR=3.22e-11)、細(xì)胞溶質(zhì)(Cytosol)(GO: 0005829, FDR=2.50e-05)等細(xì)胞成分也有不同程度的減少(見圖6)。以上結(jié)果表明，981品系在四分孢子形成、成熟至放散的過程里，在mRNA翻譯的水平和生物活性大分子的合成上都沒有魯龍1號(hào)品系活躍；另外，魯龍1號(hào)品系在細(xì)胞內(nèi)調(diào)節(jié)維持細(xì)胞穩(wěn)態(tài)的過程中表現(xiàn)較好。

2.4 與981品系II期相關(guān)的模塊基因分析與富集

將與981品系II期相關(guān)的Magenta、Green和Turquoise 3個(gè)模塊的基因表達(dá)量聚類構(gòu)建熱圖(見圖7)，可以看出這3個(gè)模塊在981品系的II期普遍上調(diào)表達(dá)，尤其Green和Turquoise兩個(gè)模塊基因是特異性上調(diào)表達(dá)。

將3個(gè)模塊內(nèi)的基因進(jìn)行KEGG代謝通路富集(見圖8)。3個(gè)模塊中均有富集到Cellular Processes相關(guān)的代謝通路，如Magenta模塊富集到溶酶體(Lysosome)(ko04142，P=0.006 899 873)和干細(xì)胞全能型調(diào)節(jié)信號(hào)通路(Signaling pathways regulating pluripotency of stem cells)通路(ko04550，P= 0.117 777 5)；green模塊富集到酵母細(xì)胞周期(Cell cycle-yeast)通路(ko04111，P=0.168 894 3)和肌動(dòng)蛋白細(xì)胞骨架調(diào)節(jié)通路(Regulation of actin cytoskeleton)(ko04810，P=0.197 946 7)；turquoise模塊富集到卵母細(xì)胞減數(shù)分裂(Oocyte meiosis)(ko04114，P=0.018 493 97)、細(xì)胞周期(Cell cycle)(ko04110，P= 0.122 86)、酵母細(xì)胞周期(Cell cycle-yeast)(ko04111，P=0.142 086 7)、酵母減數(shù)分裂(Meiosis-yeast)(ko04113，P=0.143 526 6)和緊密連接(Tight junction)(ko04530，P=0.157 927 7)等代謝通路。

(B中熱圖內(nèi)每個(gè)方框內(nèi)上方數(shù)值表示相關(guān)性系數(shù)，下方數(shù)值表示顯著性。箭頭標(biāo)識(shí)為品系相關(guān)，*標(biāo)識(shí)為981.II期相關(guān)。紅框所指為與品系相關(guān)的模塊(相關(guān)系數(shù)>0.7)，綠框所指為與981.II期相關(guān)的模塊(相關(guān)系數(shù)>0.8)。The value above each cell in heatmap indicates the correlation coefficient, and the lower value indicates the significance. The arrow indicates the cultivar feature, and the * is marked as the 981.II period feature.The red box is marked as the module related to the strain (correlation coefficient>0.7), and the green box refers to the module related to the 981.II period (correlation coefficient>0.8).)

圖4 與品系相關(guān)的6個(gè)模塊的基因在龍須菜表達(dá)譜12個(gè)樣品下的聚類分析

圖5 Up組生物過程的GO富集網(wǎng)絡(luò)圖

(A：生物過程；B：細(xì)胞成分。A: Biological process; B:Cell component.)圖6 Down組的GO富集網(wǎng)絡(luò)圖

圖7 與981品系II期相關(guān)的3個(gè)模塊的基因在表達(dá)譜12個(gè)樣品下的聚類分析

(Rich factor 是該模塊中位于該通路的基因數(shù)目與所有該通路的基因數(shù)的比值。箭頭所指為屬于細(xì)胞過程的代謝通路。A：Magenta模塊；B：Green模塊；C：Turquoise模塊。Rich factor represent the ratio of the hub genes in this pathway and all the genes of this pathway. Arrows refer to metabolic pathways that are part of the cellular process. A: Magenta；B: Green；C: Turquoise.)

將3個(gè)模塊內(nèi)前30個(gè)富集到Cellular Processes相關(guān)代謝通路的基因整合匯總，共得到30個(gè)基因。將這30個(gè)基因的FPKM值調(diào)出，經(jīng)過均一化處理后構(gòu)建樣品間的表達(dá)熱圖(見圖9)。發(fā)現(xiàn)在981品系的II期樣品內(nèi)有10個(gè)基因表達(dá)量明顯下調(diào)，另外20個(gè)基因則呈現(xiàn)表達(dá)量上升的特征。

圖9 30個(gè)備選基因在表達(dá)譜樣品內(nèi)的表達(dá)熱圖

統(tǒng)計(jì)這30個(gè)基因在KEGG通路注釋中的結(jié)果發(fā)現(xiàn)有21個(gè)基因集中富集在卵母細(xì)胞減數(shù)細(xì)胞分裂(Oocyte meiosis, ko04114)、細(xì)胞周期(Cell cycle, ko04110)、酵母細(xì)胞周期(Cell cycle-yeast, ko04111)和酵母減數(shù)分裂(Meiosis-yeast, ko04113)通路中，將這21個(gè)基因分別進(jìn)行KO(Keggorthology)注釋和GO注釋(見表2)。這些特異表達(dá)的基因主要集中在細(xì)胞周期的S期和M期。其中，處于S期的組蛋白去乙?；窰DAC1-2(GL84608、K06067)、復(fù)制起始點(diǎn)識(shí)別復(fù)合物ORC4(GL79162，K02606)和DNA復(fù)制許可因子MCM2(GL86684、K02540)上調(diào)表達(dá)，表明DNA復(fù)制可能異常進(jìn)行；而細(xì)胞周期檢驗(yàn)點(diǎn)蛋白R(shí)AD24(GL84298、K06662)主要參與S/G2期轉(zhuǎn)換中的DNA損傷檢驗(yàn)點(diǎn)，該基因上調(diào)表達(dá)間接表明S期DNA復(fù)制出現(xiàn)異常。另外在G1/S期中介導(dǎo)泛素化降解途徑的SCF蛋白亞基CUL1(GL85292、K03347)下調(diào)表達(dá)，其參與G1期周期蛋白的降解，SCF改變可能影響到S期的正常調(diào)控。G2/M期的轉(zhuǎn)化通常由CDK1起到關(guān)鍵性調(diào)控作用，而其活性依賴于細(xì)胞周期蛋白cyclin B的積累，但在是981.II期樣品中CCNB1(GL81248，K05868)出現(xiàn)下調(diào)，直接影響CDK1活性，相應(yīng)地，由CDK1介導(dǎo)的底物磷酸化水平降低，可能影響核纖層解聚、核仁解體等等，進(jìn)而影響G2/M期的轉(zhuǎn)化。

細(xì)胞分裂期M期又可分為前、前中、中、后和末期五個(gè)時(shí)期。前期中的黏連蛋白亞基SCC1(GL83834、K06670)、凝縮蛋白復(fù)合物亞基CAPH(GL85004、K06676)以及染色體結(jié)構(gòu)維持蛋白SMC1(GL86206、K06636)參與姐妹染色單體黏著和染色體的凝縮，其在981品系的II期中過量表達(dá)可能導(dǎo)致染色體凝縮狀態(tài)緊密。減數(shù)分裂重組蛋白DMC1(GL70705、K10872)在減數(shù)分裂前期I中同源染色體聯(lián)會(huì)和交換重組中發(fā)揮作用，下調(diào)可能導(dǎo)致交換重組出現(xiàn)缺陷。細(xì)胞分裂從中期向后期轉(zhuǎn)換主要依賴于后期促進(jìn)復(fù)合體APC的作用，無活性的APC需要在CDC20的結(jié)合下有活性，在中期染色體排列正確后，一方面降解cyclin B使CDK1失活，解除對(duì)分離酶ESP的磷酸化作用，進(jìn)而分解SCC1使姐妹染色單體分離，另一方面APC利用泛素化降解途徑使ESP的抑制物降解，恢復(fù)ESP活性，進(jìn)而使姐妹染色單體分離[32]。而在981.II期中，盡管APC3(GL86434、K03350)上調(diào)，但作為正調(diào)控因子的CDC20(GL84674、K03363)下調(diào)，因此對(duì)ESP活化作用減弱，即使ESP1(GL86892、K02365)的表達(dá)量上調(diào)，仍有可能是處于無活性狀態(tài)，由此推測(cè)在姐妹染色單體分離階段由于分離酶的作用降低，姐妹染色體分離出現(xiàn)異常。末期的細(xì)胞分裂控制蛋白CDC15(GL85666、K06683)參與細(xì)胞分裂的退出調(diào)控，在981.II期中同樣下調(diào)。

表2 21個(gè)基因的KO和GO注釋信息以及在細(xì)胞分裂中可能發(fā)揮的作用

續(xù)表2

除此之外，還有一些參與信號(hào)轉(zhuǎn)導(dǎo)調(diào)控途徑的基因如絲氨酸/蘇氨酸蛋白激酶PPP2家族(GL75199、K04382；GL85510、K04354；GL82076、K06268)表達(dá)水平也發(fā)生變化，這些都可能是981品系在孢子形成階段出現(xiàn)異常的原因。

3 討論

目前尚未有龍須菜的全基因組測(cè)序數(shù)據(jù)，現(xiàn)有龍須菜基因序列很多是不完整的，僅僅是基因全長(zhǎng)的一部分。在進(jìn)行數(shù)字基因表達(dá)譜分析中，測(cè)序所獲得的序列信息都要以參考基因組序列為模板進(jìn)行注釋及后續(xù)定量，參考序列的質(zhì)量高低直接影響到后續(xù)分析。而經(jīng)過對(duì)數(shù)據(jù)的初步分析，發(fā)現(xiàn)龍須菜轉(zhuǎn)錄組分析可能存在污染以及龍須菜參考基因組數(shù)目太低的問題。龍須菜基因組survey數(shù)據(jù)不夠完整，預(yù)測(cè)基因數(shù)目只有3 490條[21]，嚴(yán)重低于正常藻類基因數(shù)目。采用該預(yù)測(cè)基因集會(huì)因信息缺失而影響準(zhǔn)確性或丟失重要結(jié)果。因此在構(gòu)建參考序列時(shí)采取二次mapping，以survey信息為參考根據(jù)比對(duì)到的RNA-seq的數(shù)據(jù)進(jìn)行de novo拼接，最大化保留轉(zhuǎn)錄本信息。對(duì)轉(zhuǎn)錄本參考序列的質(zhì)檢也表明龍須菜轉(zhuǎn)錄組的完整性應(yīng)該沒有問題，無基因信息的損失。

通過基因的層次聚類，共劃分出32個(gè)共表達(dá)模塊，模塊基因數(shù)從66個(gè)到2 577個(gè)不等。每個(gè)模塊代表了具有相同表達(dá)模式的基因集，模塊的數(shù)量由參數(shù)設(shè)置，模塊數(shù)量多則可能使相似表達(dá)模式的基因分離，給后續(xù)研究帶來工作量，模塊數(shù)量少則可能使不同基因混合，掩蓋基因表達(dá)信息。楊宇昕等[16]對(duì)玉米花期基因共表達(dá)模塊篩選時(shí)劃分出20個(gè)模塊；鮮小華等[33]對(duì)甘藍(lán)型油菜黃籽的研究中對(duì)1 461個(gè)低顯著性關(guān)聯(lián)基因進(jìn)行劃分得到8個(gè)共表達(dá)模塊；鞠正等[19]對(duì)番茄果實(shí)成熟的研究中得到24個(gè)共表達(dá)模塊。

根據(jù)模塊與特征性狀的相關(guān)性關(guān)聯(lián)分析，找到與品系相關(guān)的6個(gè)模塊，對(duì)模塊內(nèi)的差異表達(dá)基因進(jìn)行GO富集結(jié)果表明981品系的信號(hào)調(diào)節(jié)、細(xì)胞轉(zhuǎn)運(yùn)、細(xì)胞分化等生物學(xué)過程比較活躍，而生物合成、細(xì)胞穩(wěn)態(tài)、代謝過程等受到抑制，胞內(nèi)物質(zhì)有不同程度減少。以上結(jié)果表明，與魯龍1號(hào)品系相比，981品系在四分孢子形成、成熟至放散的過程里，mRNA翻譯水平和生物活性大分子合成都維持在較低水平，并且細(xì)胞內(nèi)穩(wěn)態(tài)紊亂，信號(hào)傳遞活躍，推測(cè)981品系的細(xì)胞活動(dòng)受到干擾而出現(xiàn)異常。在此前對(duì)魯龍1號(hào)品系的生理形狀研究中也表明魯龍1號(hào)品系具有比981品系更快的生長(zhǎng)速度和更高的總蛋白含量[11]。

981品系的低育特性直觀表現(xiàn)在成熟期III期即放散四分孢子時(shí)期，針對(duì)四分孢子形成時(shí)期即成熟前期II期的樣品分析，篩選出3個(gè)高度相關(guān)的模塊：magenta、green和turquoise模塊。從模塊內(nèi)基因的KEGG代謝通路富集結(jié)果中找出30個(gè)參與細(xì)胞過程的基因?；虻谋磉_(dá)熱圖表明有10個(gè)基因在981的II期下調(diào)而另外20個(gè)基因表現(xiàn)上調(diào)，這種表達(dá)水平上的異常變化符合要尋找的目的基因特征。對(duì)30個(gè)基因進(jìn)一步篩選得到21個(gè)參與細(xì)胞周期的基因，分析各自的功能發(fā)現(xiàn)981品系在四分孢子發(fā)育階段可能存在DNA復(fù)制、染色體凝縮狀態(tài)、交換重組缺陷、姐妹染色單體分離、細(xì)胞分裂的退出調(diào)控等異常問題。對(duì)981品系四分孢子形成過程的形態(tài)學(xué)觀察結(jié)果表明981品系形成和釋放的四分孢子多表現(xiàn)出形態(tài)畸形和發(fā)育不完全[12]，這些基因的調(diào)控異常表明981品系可能在孢子形成時(shí)期的細(xì)胞減數(shù)分裂出現(xiàn)異常而導(dǎo)致其形態(tài)變化，或可作為深入研究的備選基因。

4 結(jié)語

對(duì)前期獲得的龍須菜不同育性品系(981和魯龍1號(hào))在不同發(fā)育階段的表達(dá)譜數(shù)據(jù)，采用WGCNA構(gòu)建基因共表達(dá)網(wǎng)絡(luò)，共獲得32個(gè)模塊。其中，模塊最大基因數(shù)量為2 577個(gè)，最小基因數(shù)為66個(gè)。經(jīng)過篩選發(fā)現(xiàn)與品系相關(guān)的模塊(以0.7為閾值)共5個(gè)，涉及基因共2 301個(gè)，對(duì)模塊基因進(jìn)行GO功能富集分析，結(jié)果表明主要富集以下功能上：信號(hào)轉(zhuǎn)導(dǎo)、轉(zhuǎn)運(yùn)、細(xì)胞分化、翻譯、細(xì)胞穩(wěn)態(tài)、生物合成過程等。981品系與魯龍1號(hào)品系的富集結(jié)果對(duì)比表明981品系在孢子放散過程中生物活性不高，生物合成水平較低，內(nèi)環(huán)境穩(wěn)態(tài)較差。篩選到21個(gè)與981品系II期四分孢子形成相關(guān)的備選基因，注釋到細(xì)胞周期的不同階段，主要集中在DNA復(fù)制和染色體形態(tài)變化上。對(duì)這些基因的進(jìn)一步研究將有助于深入揭示龍須菜四分孢子發(fā)育過程中的分子調(diào)控機(jī)制。