宋韡，湯麗群，高帆*

1(山西體育職業(yè)學(xué)院，山西 太原，030006)2(山西大學(xué) 生命科學(xué)學(xué)院，山西 太原，030006)

巴夫杜氏藻(Dunaliellaparva)屬于綠藻門(Chlorophyta)，綠藻綱(Chlorophyceae)，團(tuán)藻目(Volvocales)，杜氏藻科(Dunaliellaceae)，杜氏藻屬(Dunaliella)[1]，多分布于鹽湖、海洋及濕地，為無細(xì)胞壁的單細(xì)胞真核生物，可進(jìn)行無性或有性生殖。其含有豐富的生物活性物質(zhì)和藻類色素，如多糖、甘油、β-胡蘿卜素等，在保健品、食品加工、醫(yī)療和生物柴油等方面具有廣闊的應(yīng)用前景[2]。

目前，國內(nèi)外對(duì)于巴夫杜氏藻的研究多集中于培養(yǎng)條件優(yōu)化[3]、關(guān)鍵基因克隆分析[4]、活性物質(zhì)提取等方面[5]，對(duì)于其轉(zhuǎn)錄組學(xué)方面的研究很少。進(jìn)入21世紀(jì)以來，隨著國內(nèi)外高通量測(cè)序技術(shù)快速發(fā)展，以Ion Torrent、Roche 454、Thermo Fisher、SOLiD、HiSeq為代表的測(cè)序平臺(tái)[6]被廣泛使用。該技術(shù)已成為在不同生長階段、不同脅迫處理等條件下獲取物種關(guān)鍵功能基因，預(yù)測(cè)代謝通路及表征生物學(xué)過程的重要手段之一[7]。YAO等[8]對(duì)Dunaliellatertiolecta進(jìn)行轉(zhuǎn)錄組測(cè)序，獲得10 185條與脂肪酸合成相關(guān)的基因。PUENTE-SNCHEZ等[9]對(duì)Dunaliellaacidophila進(jìn)行測(cè)序，篩選出不同重金屬脅迫下的關(guān)鍵基因。HE等[10]對(duì)D.salina在高鹽度脅迫下的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行分析，獲取了耐鹽關(guān)鍵基因及其代謝通路。上述轉(zhuǎn)錄組測(cè)序分析均基于杜氏藻屬的其他種類，對(duì)D.parva的研究還僅限于2016年尚?；▽?duì)限氮條件下巴夫杜氏藻的轉(zhuǎn)錄組和蛋白組測(cè)序分析[11]。該研究通過Illumina HiSeq 2000測(cè)序平臺(tái)預(yù)測(cè)了D.parva油脂合成途徑，鑒定了油脂合成關(guān)鍵基因wri1和Lip，但并未對(duì)其不同生長時(shí)期的轉(zhuǎn)錄水平及相關(guān)基因進(jìn)行分析。BGISEQ-500是新一代高通量測(cè)序平臺(tái)，利用該平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序研究已在多個(gè)物種中報(bào)道[12]。本研究利用該測(cè)序平臺(tái)對(duì)不同生長時(shí)期巴夫杜氏藻進(jìn)行轉(zhuǎn)錄組測(cè)序分析，獲得大量差異表達(dá)基因，并預(yù)測(cè)其功能及參與的代謝通路。研究結(jié)果為該藻生長發(fā)育調(diào)控機(jī)制及相關(guān)代謝途徑的研究奠定了基礎(chǔ)，對(duì)該藻類資源改良與分子選育具有重要的指導(dǎo)意義。

1 材料與方法

1.1 藻種收集與培養(yǎng)

巴夫杜氏藻種(D.parva)，英國CCAP藻種庫，藻細(xì)胞活化后用優(yōu)化杜氏鹽藻培養(yǎng)基進(jìn)行光照靜態(tài)培養(yǎng)，培養(yǎng)基配方：FeC6H5O70.005 g，NaH2PO4·2H2O 0.015 g，CaCl2·2H2O 0.045 g，KCl 0.075 g，NaNO30.45 g，NaHCO30.85 g， MgSO4·7H2O 1.25 g，NaCl 130 g，Co(NO3)2·6H2O 5 mg，CuSO4·5H2O 8 mg，ZnSO4·7H2O 23 mg，維生素B10.5 mg，pH調(diào)至7.5，ddH2O定容至1 000 mL。光照強(qiáng)度20 000 lx，溫度(25±5) ℃，光周期16 h∶8 h，培養(yǎng)室保持持續(xù)無菌通風(fēng)。每隔7 d采樣進(jìn)行顯微鏡檢，每隔21 d進(jìn)行新鮮培養(yǎng)基補(bǔ)充，在培養(yǎng)液總體積不變的前提下確保培養(yǎng)過程無污染。

1.2 生長周期測(cè)定

參考王婷等[13]對(duì)杜氏藻生理指標(biāo)的研究方法，利用紫外-可見光分光光度計(jì)測(cè)定685 nm處不同培養(yǎng)時(shí)間內(nèi)的巴夫杜氏藻培養(yǎng)液吸光值。根據(jù)OD685數(shù)值大小，繪制巴夫杜氏藻生長周期曲線，測(cè)定其生長初期(包括停滯期和對(duì)數(shù)期)、中期(穩(wěn)定期)和后期(衰亡期)的時(shí)間范圍。

1.3 RNA提取，cDNA文庫構(gòu)建及轉(zhuǎn)錄組測(cè)序

分別取3份上述不同生長時(shí)期(OD685平均值所在時(shí)期)的巴夫杜氏藻液各300 mL，4 ℃低溫離心(1 500×g)后用液氮速凍并快速研磨，用TRIzol試劑盒分別提取藻株總RNA，Agilent 2100 Bioanalyzer對(duì)其質(zhì)量和濃度進(jìn)行評(píng)估，利用反轉(zhuǎn)錄試劑盒反轉(zhuǎn)錄為cDNA，SMART cDNA文庫構(gòu)建試劑盒對(duì)合格樣本進(jìn)行文庫構(gòu)建。構(gòu)建好的文庫在BGISEQ-500測(cè)序平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序，上機(jī)測(cè)序操作委托深圳華大基因完成。

1.4 測(cè)序數(shù)據(jù)分析及de novo組裝

利用SOAPnuke[14]對(duì)測(cè)序數(shù)據(jù)進(jìn)行過濾與篩選，具體流程是：刪除5′和3′端含接頭及poly-N(N>5%)的低質(zhì)量reads，篩選獲得的Clean reads(Q20值> 95%，Q30值> 85%)，利用Trinity軟件對(duì)Clean reads進(jìn)行從頭組裝。利用Bowtie2軟件[15]預(yù)測(cè)Clean reads中的UniGenes，利用Tgicl軟件[16]分析UniGenes數(shù)量和長度，利用BUSCO軟件[17]評(píng)估組裝質(zhì)量。

1.5 轉(zhuǎn)錄組數(shù)據(jù)注釋

基于轉(zhuǎn)錄組UniGenes與Gene Ontology(GO)、RefSeq Nonredundant Protein(NR)、Nucleotide Sequence(NT)、SwissPort、Pfam、Kyoto Encyclopedia of Genes and Genomes(KEGG)及Eukaryotic Orthologous Group(KOG)數(shù)據(jù)庫數(shù)據(jù)的比對(duì)相似度(參數(shù)設(shè)置：E值≤1.00E-10，查詢序列長度百分比≥85%, 無冗余Contigs)。基于Pfam和KOG數(shù)據(jù)庫，利用HMMER程序[18]預(yù)測(cè)UniGene潛在功能及抗病毒結(jié)構(gòu)域。

1.6 基因表達(dá)分析

利用Fragments per Kilobase of Transcript per Million Fragments(FPKM)[19]測(cè)算UniGene表達(dá)量，利用RSEM軟件[20]測(cè)算差異表達(dá)基因(differential expression gene，DEG)，利用DEseq2軟件構(gòu)建DEGs分布火山圖(log2fold change≥|2.00|，Q值≤ 0.05)，利用R語言構(gòu)建不同樣本的DEGs聚類熱圖。

1.7 GO和KEGG通路富集

利用Blast2GO[21]對(duì)UniGene進(jìn)行GO條目分配和Ontology分類，基于R語言的Phyper功能對(duì)GO條目進(jìn)行富集[錯(cuò)誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)≤0.01]。利用KEGG通路的KAAS程序?qū)niGene進(jìn)行KO條目分配和通路分類，基于R語言的Phyper功能對(duì)KEGG通路進(jìn)行富集(FDR≤0.01)。

2 結(jié)果與分析

2.1 巴夫杜氏藻生長周期鑒定

基于巴夫杜氏藻的3個(gè)生物學(xué)重復(fù)生長曲線擬合，研究發(fā)現(xiàn)其生長周期大體可分為3個(gè)明顯不同的階段，即生長初期(含延滯期和對(duì)數(shù)期)、生長中期(成熟期或穩(wěn)定期)和生長末期(衰亡期或凋亡期)。由圖1-a可知，生長初期的平均OD685值為0.49±0.03，周期為0～45 d；生長中期的平均OD685值約為1.63±0.02，周期約為45～55 d；生長末期的平均OD685值為1.59±0.02，周期約為55～65 d。由圖1-b可知，細(xì)胞色素在生長初期積累較少，整體呈淺綠色，藻細(xì)胞生長速度較快，細(xì)胞呈小球狀；細(xì)胞色素在生長中期積累較多，整體呈深綠色，藻細(xì)胞呈橢球狀，生長速度趨于穩(wěn)定；細(xì)胞色素在生長末期的積累更深，整體呈墨綠色，藻細(xì)胞呈扁球狀，生長速度趨于延滯。

2.2 RNA測(cè)序與de novo組裝

巴夫杜氏藻RNA平均質(zhì)量濃度為56 ng/μL，RNA完整值(RNA integrity number，RIN)平均值為7.1，符合RNA測(cè)序標(biāo)準(zhǔn)(RNA質(zhì)量濃度≥ 30 ng/μL，RIN ≥ 6.0)。3個(gè)生長時(shí)期(30、50、60 d)共9個(gè)生物樣本的原始測(cè)序數(shù)據(jù)共18.89 Gb，其中原始Reads約665.7 Mb，平均每個(gè)樣本約71.76 Mb的Clean reads(總占比約72.23%)(表1)。共鑒定90 153條UniGene，檢測(cè)到272條完整的Contig片段(完整度占比約90.23%，總長約7.65×107bp，平均長度約812 bp)(圖2-a，表1)。N50長1 305 bp(N70長783 bp，N90長345 bp)，GC占比約52.32%。UniGene中共預(yù)測(cè)到49 643條CDS(總長為3.90×107bp，N50為982 bp，GC占比53.65%)，23 528條SSR(總長為20 412 bp)。12個(gè)測(cè)序樣本的Q20和Q30百分比平均值分別為98.56%和89.43%。denovo組裝的BUSCO質(zhì)量評(píng)估結(jié)果顯示，組裝完整的測(cè)序片段、單拷貝及雙拷貝片段總占比>90%(圖2-b)，組裝質(zhì)量較高。

a-巴夫杜氏藻生長周期；b-巴夫杜氏藻形態(tài)學(xué)觀察(標(biāo)尺20 μm)圖1 巴夫杜氏藻生長周期及其細(xì)胞形態(tài)Fig.1 Growth cycle and cell morphology of D. parva

a-Contig長度分布；b-基于BUSCOs的de novo組裝質(zhì)量評(píng)估圖2 轉(zhuǎn)錄組測(cè)序Contig長度及de novo組裝質(zhì)量分析Fig.2 Analysis of transcriptomic Contig length and de novo assembly quality

表1 巴夫杜氏藻轉(zhuǎn)錄組測(cè)序de novo組裝基本信息Table 1 D.parva transcriptome de novo assembly information

2.3 基因注釋與功能預(yù)測(cè)

90 153條UniGenes中平均約66.36%的基因可在七大生物信息數(shù)據(jù)庫注釋到(表2)，其中約5.10%的基因在數(shù)據(jù)庫中均有注釋信息。每四大數(shù)據(jù)庫的UniGene功能注釋數(shù)詳見圖3?；贜R的UniGne注釋物種分類結(jié)果顯示，約23.43%的UniGene可在藻類(Chlamydomonaseustigma，Goniumpectorale和Volvoxcarter)和陸生植物(Dorcocerashygrometricum和Ricinuscommunis)中注釋到，剩余76.57%的數(shù)據(jù)未在已報(bào)道物種基因信息庫中檢測(cè)到注釋信息(圖4-a)。注釋信息最多的812條UniGenes中可預(yù)測(cè)到可變的N端病毒結(jié)構(gòu)域，其次為510條UniGenes中預(yù)測(cè)到的受體樣蛋白(receptor-like protein，RLP)結(jié)構(gòu)域，僅僅3條UniGnes功能被預(yù)測(cè)與mol-like結(jié)構(gòu)域形成相關(guān)(圖4-b)。最多的4 650條UniGenes在KOG數(shù)據(jù)庫中的功能注釋為一般性調(diào)控功能，3 300條UniGenes的功能與后轉(zhuǎn)錄調(diào)控與蛋白質(zhì)周轉(zhuǎn)有關(guān)，最少的UniGne功能被預(yù)測(cè)與細(xì)胞運(yùn)動(dòng)與硫的形成相關(guān)(圖4-c)。

表2 巴夫杜氏藻的UniGene注釋信息Table 2 UniGene annotation of D.parva

2.4 基因表達(dá)與聚類分析

由圖5-a可知，基于UniGene表達(dá)FPKM值，巴夫杜氏藻生長中期和末期的基因整體可聚為一類，具有相似的表達(dá)趨勢(shì)，生長初期的基因單獨(dú)聚為一類，所有UniGene大致可劃分為四大簇(Cluster Ⅰ、Ⅱ、Ⅲ和Ⅳ)。以巴夫杜氏藻生長初期UniGene為對(duì)照，共檢測(cè)到86 387條差異表達(dá)基因(圖5-b和圖5-c)。相對(duì)于生長初期，生長中期共檢測(cè)到上調(diào)基因14 556條，下調(diào)基因22 589條；生長末期共檢測(cè)到上調(diào)基因18 940條，下調(diào)基因30 302條。兩對(duì)比較組間共有差異表達(dá)基因有20 579條(圖5-d)。

a-GO-NR-KEGG-KOG注釋維恩圖；b-GO-SwissProt-KEGG-NT 注釋維恩圖；c-GO-Pfam-KEGG-NR注釋維恩圖圖3 巴夫杜氏藻unigene注釋維恩圖Fig.3 Venn diagram of D.parva unigenes

a-巴夫杜氏藻UniGene物種分類注釋；b-巴夫杜氏藻UniGene 抗病毒基因結(jié)構(gòu)域注釋；c-巴夫杜氏藻UniGene KOG功能注釋圖4 巴夫杜氏藻UniGene注釋Fig.4 UniGene annotation of D.parva

2.5 GO分類與富集

由圖6-a和圖6-b可知，在生物學(xué)過程類別中，GO條目最多的是細(xì)胞過程(生長中期vs.生長初期的基因數(shù)為2 250，生長末期vs.生長初期的基因數(shù)為6 300)；在細(xì)胞組分類別中，GO條目最多的是細(xì)胞(生長中期vs.生長初期的基因數(shù)為1 968，生長末期vs.生長初期的基因數(shù)為4 987)；在分子功能類別中，GO條目最多的是催化活性(生長中期vs.生長初期間的基因數(shù)為2 419，生長末期vs.生長初期的基因數(shù)為7 962)。由兩對(duì)巴夫杜氏藻比較組的差異基因GO富集top 20結(jié)果可知(圖6-c和圖6-d)，富集程度最高的GO條目均是膜的整體成分，生長中期vs.生長初期的基因數(shù)最高為1 140，富集比值最高為0.42；生長末期vs.生長初期的基因數(shù)最高為3 469，富集比值最高為0.41。

2.6 KEGG通路分類與富集

由圖7-a和圖7-b可知，2對(duì)比較組間所有的通路均可劃分為細(xì)胞過程、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機(jī)體系統(tǒng)。其中，通路基因條目最多的是整體代謝通路(生長中期vs.生長初期的基因數(shù)為1 200，生長末期vs.生長初期的基因數(shù)為4 300)。由兩對(duì)巴夫杜氏藻比較組的差異基因KEGG通路富集top 20結(jié)果可知(圖7-c和圖7-d)，生長中期vs.生長初期富集程度最高的通路是剪接體(基因數(shù)為546，富集比值為0.16)；生長末期vs.生長初期富集程度最高的通路是RNA轉(zhuǎn)運(yùn)(基因數(shù)最高為1 521，富集比值最高為0.15)。

a-不同生長時(shí)期的基因表達(dá)聚類熱圖；b-生長中期與生長初期間差異基因表達(dá)散點(diǎn)圖；c-生長末期與生長初期間差異基因 表達(dá)散點(diǎn)圖；d-兩對(duì)比較組間的共有差異基因維恩圖圖5 巴夫杜氏藻不同比較組間的差異表達(dá)基因分析Fig.5 Analysis of differentially expressed genes among different D.parva comparison groups

a-生長中期與生長初期間差異基因GO分類；b-生長末期與生長初期間差異基因GO分類；c-生長中期與生長初期間差異基因 GO富集；d-生長末期與生長初期間差異基因GO富集圖6 巴夫杜氏藻不同比較組間的差異基因GO分類與富集Fig.6 GO classification and enrichment of DEGs among different D.parva comparison groups

a-生長中期與生長初期間差異基因KEGG通路分類；b-生長末期與生長初期間差異基因KEGG通路分類；c-生長中期與生長初期間 差異基因KEGG通路富集；d-生長末期與生長初期間差異基因KEGG通路富集圖7 巴夫杜氏藻不同比較組間的差異基因KEGG通路分類與富集Fig.7 KEGG pathway classification and enrichment of DEGs among different D.parva comparison groups

3 討論

目前，轉(zhuǎn)錄組測(cè)序技術(shù)已成為闡明藻類功能基因的重要手段之一。本研究通過第二代高通量測(cè)序平臺(tái)BGISEQ-500對(duì)不同生長時(shí)期的巴夫杜氏藻進(jìn)行轉(zhuǎn)錄組測(cè)序，獲得的單樣本Clean reads(71.76 Mb)及GC百分比(52.32%)均高于PANAHI等[22]利用Illumina GAIIx測(cè)序平臺(tái)獲得的鹽生杜氏藻(D.salina)數(shù)據(jù)(Clean reads 31.2 Mb，GC百分比50%)。巴夫杜氏藻cDNA測(cè)序文庫的Q20和Q30平均值高于POLLE等[23]報(bào)道的D.salinaCCAP 19/18轉(zhuǎn)錄組測(cè)序的Q20和Q30值。但本研究預(yù)測(cè)的CDS數(shù)(49 643條)卻低于HE等[10]在Illumina HiSeq 2 000測(cè)序平臺(tái)獲得的D.salinaCCAP 19/30CDS數(shù)(49 751條)。巴夫杜氏藻轉(zhuǎn)錄組測(cè)序數(shù)據(jù)denovo組裝后的N50長度(1 305 bp)小于HONG等[7]從D.salina435轉(zhuǎn)錄組測(cè)序獲得的N50長度(1 727 bp)。由此可知，基于BGISEQ-500測(cè)序平臺(tái)獲取的巴夫杜氏藻轉(zhuǎn)錄組測(cè)序數(shù)據(jù)整體是可靠的。

僅23.43%的巴夫杜氏藻UniGene可在植物物種中注釋到，相較于其他藻類物種注釋信息，比例較低(一般>40%)[24]。筆者推測(cè)其原因可能是：(1)該情況可能與巴夫杜氏藻原始系統(tǒng)進(jìn)化地位相關(guān)。巴夫杜氏藻是綠藻的一種，研究顯示綠藻起源較早，相比高等陸生植物而言，其與藍(lán)細(xì)菌甚至與某些原生動(dòng)物系統(tǒng)進(jìn)化關(guān)系較近[25]，而本研究未將UniGene的物種注釋范圍擴(kuò)大到細(xì)菌及原生動(dòng)物。(2)本研究測(cè)序獲得UniGene數(shù)據(jù)未經(jīng)開放型閱讀框(open reading box，ORF)預(yù)測(cè)處理，由此冗余數(shù)據(jù)及非功能基因的片段對(duì)物種注釋信息造成了干擾。因此，在原UniGene數(shù)據(jù)基礎(chǔ)上，擴(kuò)大其物種基因注釋的數(shù)據(jù)庫范圍，并進(jìn)一步處理搜集含ORF框的功能基因數(shù)據(jù)，將可能提高其基因的物種注釋比例。

以生長初期為對(duì)照，巴夫杜氏藻生長中期差異表達(dá)基因數(shù)低于生長末期，這表明到了藻類的生長末期，隨著細(xì)胞的凋亡、代謝廢物及有毒物質(zhì)的積累，用于調(diào)控生長危機(jī)的基因數(shù)明顯增加，以此應(yīng)答不利于藻類營養(yǎng)物質(zhì)積累及細(xì)胞增殖的細(xì)胞內(nèi)環(huán)境[26]。此外，相對(duì)于生長初期，巴夫杜氏藻生長中期和生長末期的下調(diào)基因數(shù)均高于上調(diào)基因數(shù)，這暗示了隨著生長周期延長，對(duì)藻株生長發(fā)育和新陳代謝等生物學(xué)過程起重要作用的基因來說，其總體調(diào)控方式是負(fù)向的，這與杜氏藻在不同非生物脅迫下的差異基因統(tǒng)計(jì)結(jié)果有所差異(一般上調(diào)基因數(shù)高于下調(diào)基因)[8-10]。

不同比較組間基因GO條目與功能富集多數(shù)與細(xì)胞的完整性及其功能組分相關(guān)，這表明隨著生長周期延長，巴夫杜氏藻胞內(nèi)的甘油、多糖、β-胡蘿卜素等生物活性內(nèi)含物逐漸積累，細(xì)胞體積不斷增大[27]，相關(guān)基因的功能多與細(xì)胞周期調(diào)控、細(xì)胞組分合成及催化活性緊密相關(guān)。不同比較組間基因參與的KEGG通路與多種代謝過程相關(guān)，但其詳細(xì)的代謝途徑仍不明確。究其原因，可能與已報(bào)道藻類代謝通路信息較少有關(guān)。

4 結(jié)論

不同生長時(shí)期巴夫杜氏藻轉(zhuǎn)錄組測(cè)序結(jié)果及組裝質(zhì)量較高，上調(diào)表達(dá)基因數(shù)高于下調(diào)表達(dá)基因，GO富集條目最多的是細(xì)胞過程、細(xì)胞整體成分及催化活性，KEGG通路富集程度最高的是剪接體及RNA轉(zhuǎn)運(yùn)代謝通路。下一步，對(duì)該藻株生長發(fā)育關(guān)鍵基因的挖掘及其表達(dá)驗(yàn)證，蛋白組與代謝組學(xué)的聯(lián)合分析其參與的關(guān)鍵代謝通路很有必要，這對(duì)于深入揭示巴夫杜氏藻的生長發(fā)育分子機(jī)制，分子水平調(diào)控其活性物質(zhì)含量，遺傳水平改良其資源品質(zhì)具有重要的意義。