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

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

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

近年來測序技術(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)而計算不同節(jié)點的相異系數(shù)并形成相異矩陣，根據(jù)相異矩陣進(jìn)行層次聚類后依此將基因劃分成以樞紐基因(Hub gene)為核心的基因共表達(dá)模塊(Module)，進(jìn)一步將模塊與性狀特征相關(guān)聯(lián)從而提取出目的信息[13]。這一技術(shù)被廣泛應(yīng)用于動植物性狀的調(diào)控機(jī)制解析上，如對玉米花期基因[16]、胡楊異形葉發(fā)育[17]、蘋果黃皮突變[18]、番茄果實成熟[19]以及雨生紅球藻中蝦青素的積累[20]等過程的相關(guān)分析中。

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

1 材料與方法

1.1 材料培養(yǎng)

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

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

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

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

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

對原始測序序列(Raw reads)進(jìn)行質(zhì)量控制，過濾掉帶接頭以及低質(zhì)量的序列得到clean reads。使用Hisat2[28]軟件將clean reads與龍須菜Survey數(shù)據(jù)[21]比對，將比對成功的reads挑出后用SOAPdenovo-Trans進(jìn)行拼接構(gòu)建參考序列(Reference)，拼接后過濾掉小于200 bp的contig。使用bowtie對每個樣品的reads進(jìn)行比對并使用RSEM統(tǒng)計樣品中基因的read count數(shù)目，edgeR[29-30]進(jìn)行差異基因分析，TMM算法對read count數(shù)換算成FPKM (fragments per kilobase of exon per million fragments mapped, 每千堿基外顯子百萬片段數(shù))值來估算基因表達(dá)水平并進(jìn)行矩陣標(biāo)準(zhǔn)化矯正。對差異表達(dá)基因的表達(dá)矩陣進(jì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á)量表對數(shù)據(jù)進(jìn)行過濾后構(gòu)建表達(dá)矩陣，并對樣品進(jìn)行聚類分析；利用WGCNA 中的pickSoftThreshold函數(shù)計算軟閾值β參數(shù)使得滿足無尺度網(wǎng)絡(luò)分布前提條件；選取R2達(dá)到0.9 時的閾值參數(shù)β[13]，計算所有基因間的相異系數(shù)創(chuàng)建鄰接矩陣，并用TOMsimilarity函數(shù)將鄰接矩陣轉(zhuǎn)換為拓?fù)渚仃嚭螅∧娴玫较喈愋跃仃嚥⒁源藰?gòu)建系統(tǒng)聚類樹。設(shè)定每個模塊最小基因數(shù)為30 個，利用cutreeDynamic動態(tài)剪枝法修剪聚類樹識別初始共表達(dá)模塊。隨后用moduleEigengene函數(shù)計算每個模塊的特征向量值ME (moduleEigengene)，用mergeCloseModules函數(shù)選擇相似性在100%以上的模塊進(jìn)行合并融合(cutHeight=0.0)。

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

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

利用Alldegrees1函數(shù)計算模塊內(nèi)基因間連接數(shù)，篩選模塊內(nèi)連接數(shù)前10%的基因為Hub基因。

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

對模塊內(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ù)處理及去除離群樣本

為了排除污染和盡可能的利用基因組信息，對測序數(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。并對轉(zhuǎn)錄組reference分別進(jìn)行GO和KEGG代謝通路注釋(見圖1)。KEGG結(jié)果中注釋到一級大類為新陳代謝(Metabolism)、遺傳信息過程(Genetic information processing)、環(huán)境信息過程(Environmental information processing)、細(xì)胞過程(Celluar processing)、生物系統(tǒng)(Organismal systems)。注釋到的基因數(shù)最多的二級分類為翻譯(Translation)過程，有224個基因(見圖1A)。GO結(jié)果中多數(shù)基因富集到代謝過程(Metabolic process)、細(xì)胞過程(Cellular process)、細(xì)胞(Cell)、細(xì)胞部分(Cell part)、催化活性(Catalytic activity)、結(jié)合(Binding)等Go term(見圖1B)。后續(xù)樣本定量分析與功能注釋則采用此reference。

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

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

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

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

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

經(jīng)篩選發(fā)現(xiàn)與品系(Cultivar)相關(guān)的模塊有 6個(圖3B中相關(guān)性系數(shù)絕對值高于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)五個模塊，以及一個備選模塊 Purple模塊(cor=0.68，P=0.02)。相關(guān)性為正值，表示基因在981品系中呈現(xiàn)較高表達(dá)水平(與魯龍1號相比)，相關(guān)性為負(fù)值，表示基因在981品系中呈現(xiàn)較低表達(dá)水平(與魯龍1號相比)；P<0.5顯示相關(guān)性分析可靠。

另一方面，由于981品系表現(xiàn)出四分孢子放散時期的低育特性，因此選擇四分孢子形成的時期即981.II期作為關(guān)注對象，發(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個模塊。相關(guān)性的正負(fù)表示表達(dá)水平的高低。

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

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

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

圖5中有顏色的圓圈表示有基因富集，顏色越深顯著性越高，圓圈越大相應(yīng)富集的基因越多。981中高表達(dá)的基因中(Up組)，主要富集到信號轉(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品系在四分孢子形成、成熟至放散過程中，信號調(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號品系活躍；另外，魯龍1號品系在細(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個模塊的基因表達(dá)量聚類構(gòu)建熱圖(見圖7)，可以看出這3個模塊在981品系的II期普遍上調(diào)表達(dá)，尤其Green和Turquoise兩個模塊基因是特異性上調(diào)表達(dá)。

將3個模塊內(nèi)的基因進(jìn)行KEGG代謝通路富集(見圖8)。3個模塊中均有富集到Cellular Processes相關(guān)的代謝通路，如Magenta模塊富集到溶酶體(Lysosome)(ko04142，P=0.006 899 873)和干細(xì)胞全能型調(diào)節(jié)信號通路(Signaling pathways regulating pluripotency of stem cells)通路(ko04550，P= 0.117 777 5)；green模塊富集到酵母細(xì)胞周期(Cell cycle-yeast)通路(ko04111，P=0.168 894 3)和肌動蛋白細(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)每個方框內(nèi)上方數(shù)值表示相關(guān)性系數(shù)，下方數(shù)值表示顯著性。箭頭標(biāo)識為品系相關(guān)，*標(biāo)識為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個模塊的基因在龍須菜表達(dá)譜12個樣品下的聚類分析

圖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個模塊的基因在表達(dá)譜12個樣品下的聚類分析

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

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

統(tǒng)計這30個基因在KEGG通路注釋中的結(jié)果發(fā)現(xiàn)有21個基因集中富集在卵母細(xì)胞減數(shù)細(xì)胞分裂(Oocyte meiosis, ko04114)、細(xì)胞周期(Cell cycle, ko04110)、酵母細(xì)胞周期(Cell cycle-yeast, ko04111)和酵母減數(shù)分裂(Meiosis-yeast, ko04113)通路中，將這21個基因分別進(jìn)行KO(Keggorthology)注釋和GO注釋(見表2)。這些特異表達(dá)的基因主要集中在細(xì)胞周期的S期和M期。其中，處于S期的組蛋白去乙?；窰DAC1-2(GL84608、K06067)、復(fù)制起始點識別復(fù)合物ORC4(GL79162，K02606)和DNA復(fù)制許可因子MCM2(GL86684、K02540)上調(diào)表達(dá)，表明DNA復(fù)制可能異常進(jìn)行；而細(xì)胞周期檢驗點蛋白RAD24(GL84298、K06662)主要參與S/G2期轉(zhuǎn)換中的DNA損傷檢驗點，該基因上調(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期又可分為前、前中、中、后和末期五個時期。前期中的黏連蛋白亞基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)會和交換重組中發(fā)揮作用，下調(diào)可能導(dǎo)致交換重組出現(xiàn)缺陷。細(xì)胞分裂從中期向后期轉(zhuǎn)換主要依賴于后期促進(jìn)復(fù)合體APC的作用，無活性的APC需要在CDC20的結(jié)合下有活性，在中期染色體排列正確后，一方面降解cyclin B使CDK1失活，解除對分離酶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)，因此對ESP活化作用減弱，即使ESP1(GL86892、K02365)的表達(dá)量上調(diào)，仍有可能是處于無活性狀態(tài)，由此推測在姐妹染色單體分離階段由于分離酶的作用降低，姐妹染色體分離出現(xiàn)異常。末期的細(xì)胞分裂控制蛋白CDC15(GL85666、K06683)參與細(xì)胞分裂的退出調(diào)控，在981.II期中同樣下調(diào)。

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

續(xù)表2

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

3 討論

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

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

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

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

4 結(jié)語

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