裂殖壺菌誘變株高產DHA機理的轉錄組學分析

胡愛云 ，龔定芳 ，李 莎

(1.無錫太湖學院 江蘇省物聯網應用技術重點建設實驗室，江蘇 無錫 214064；2.江南大學 生物工程學院，工業生物技術教育部重點實驗室，江蘇 無錫 214122)

隨著科技的發展，人們對DHA的生理功能有了較為詳細的了解[1-3]。DHA的來源主要分為兩類：一類是傳統魚油DHA[4]，主要從DHA含量較高的深海魚類中提取獲得；另一類是微生物發酵法獲取DHA[5]。與其他海洋真菌相比，裂殖壺菌不僅能夠大量積累油脂，還可以積累許多對人體有益的活性物質，具有生長速度快、DHA含量高等優勢，市場前景廣闊。裂殖壺菌藻油的安全性已經得到了美國FDA的認可，在我國衛生部2012年頒布的相關政策中提到允許在奶粉中添加裂殖壺菌藻油DHA[6]。在前期研究中，利用ARTP誘變技術獲取一株高產DHA的裂殖壺菌SR21突變菌株，被命名為NA2[7]，但其高產DHA的分子機理尚不明確。

轉錄組學已經被廣泛用于各項科學研究，目前已經在裂殖壺菌的研究中得到運用。王康[8]使用基于RNA-seq的轉錄組學技術，對添加了萘氧乙酸或茉莉酸培養24 h的裂殖壺菌及正常培養24 h的裂殖壺菌進行轉錄組分析，根據基因差異表達的情況，證實了磷酸戊糖途徑和糖酵解途徑增強，氨基酸合成代謝途徑減弱為油脂合成提供了更多的前體物質乙酰輔酶A，分析了提高DHA產量的機理。陳偉[9]基于比較轉錄組學，對甘油和葡萄糖兩種碳源條件下的差異基因的表達進行了全面的分析，對甘油促進DHA大量積累的機理進行了解析。Ma等[10]基于轉錄組分析考察了低溫情況下基因差異表達情況，結果表明低溫導致脂肪酸合酶和3種多不飽和脂肪酸合酶亞基的表達分別下調約3.7倍和2倍。這些差異基因介導了冷適應的關鍵步驟，如信號轉導、細胞成分生物發生、碳水化合物和脂質代謝等。

本研究將利用RNA-seq技術分析ARTP誘變菌株NA2高產DHA的分子機理，并利用熒光定量PCR對基因的差異表達情況進行驗證。

1 材料與方法

1.1 實驗材料

菌種：裂殖壺菌菌株(SchizochytriumlimacinumSR21)，購自美國菌種保藏中心(ATCC)，-80℃保存于ATCC By+790培養基中。活化培養基：ATCC By+790固體培養基。種子培養基：葡萄糖 30 g/L，酵母粉6 g/L，魚粉蛋白胨6 g/L。搖瓶培養基：參考文獻[11] 。誘變菌株NA2，在前期研究中由ARTP誘變得到[7]。

1.2 轉錄組分析方法

1.2.1 RNA測序及數據處理

RNA的提取：取搖瓶培養96 h的原始菌株SR21和誘變菌株NA2發酵液4 mL，10 000 r/min離心3 min，利用UNIQ-10柱式Trizol試劑盒提取總RNA。

RNA-seq步驟參照文獻[12]。

測序數據的處理：Illumina HiSeq上機測序得到原始測序序列(Raw Reads)，為了確保信息分析準確性，必須對原始測序序列過濾，首先檢查測序錯誤率分布及A/T/G/C 含量分布，然后去除帶接頭的、質量低的以及N(N表示無法確定堿基信息)所占比例大于10%的測序序列，得到過濾后序列。再采用Trinity[13]對過濾后序列進行拼接獲得轉錄本，作為后續分析的參考序列，用Corset[14]程序對轉錄本進行層次聚類，以聚類后得到的每條基因中最長的轉錄本作為該基因的代表，稱為Gene，以此進行后續的分析。

1.2.2 測序數據分析

參照七大數據庫(參照文獻[15] )對測序的結果進行基因功能注釋。

采用DESeq2[16]進行差異基因篩選，篩選閾值為差異顯著性p值小于0.05(padj<0.05)且對差異倍數取以2為底的對數之后的值大于1(|log2FoldChange|>1)。利用Pathway顯著性富集[17]進一步對差異表達基因進行Kegg富集分析。

1.2.3 實時熒光定量PCR驗證(RT-qPCR)

按照1.2.1所述提取總RNA，再將其反轉錄成cDNA，反轉錄試劑盒為南京諾唯贊公司的HiScript III RT SuperMix for qPCR(+gDNA wiper)。對cDNA進行梯度稀釋，將稀釋后的模板進行RT-qPCR擴增，驗證引物的擴增效率。采用染料法進行RT-qPCR，驗證差異基因的差異表達倍數。其中以裂殖壺菌18S rRNA為內參，引物設計采用線上設計，網址為https://sg.idtdna.com/scitools/Applications/RealTimePCR/。RT-qPCR反應體系如表1所示。檢測儀器為CFX96 Real-Time system，設置PCR的程序為：95℃預熱30 s，95℃變性5 s后60℃退火延伸30 s并循環40次。融解曲線參數參照試劑盒設置。

表1 RT-qPCR反應體系

最終以2-ΔΔCt來表征基因的差異表達倍數，其中Ct值是指在PCR過程中，擴增產物的熒光信號達到設定的閾值時所需要的循環次數[18]。每個樣品重復3次。

2 結果與討論

2.1 總RNA樣品質量評估

用搖瓶將原始菌株SR21和誘變菌株NA2培養96 h，分別取發酵液4 mL，提取總RNA。利用以下方法檢測所獲取的總RNA的質量：①瓊脂糖凝膠電泳檢測RNA降解程度、污染程度；②Nanodrop檢測RNA的純度(OD260/OD280)；③Qubit對RNA濃度進行準確定量；④利用Agilent 2100精確檢測RNA的完整性。進行兩次平行實驗，具體結果見表2。

通常認為，OD260/OD280處于1.8～2.2之間，OD260/OD230大于2時核酸純度較高，RNA完整性(RNA Integrity Number，RIN)越接近10則表示RNA越完整，28S/18S越接近2則表示RNA質量越好[19]。由表2可知，綜合各項考察指標，本研究中所獲取的4個RNA樣品均滿足RNA-seq的要求。

2.2 測序數據質量評估

Illumina HiSeq上機測序得到原始測序序列，檢查堿基錯誤率分布及A/T/G/C 含量分布。過濾原始數據，去除帶接頭的、低質量的以及N(N表示堿基信息沒辦法確定的)所占比例大于10%的序列，最終得到過濾后數據。具體的評估計算結果如表3所示。

注：Q20、Q30表示堿基質量大于20或30所占的比例；GC表示的是檢測片段中G和C的含量。

從表3可以看出，堿基錯誤率低于0.1%，過濾后序列數據量大，Q20和Q30較大，均高于90%。

進一步對過濾后序列進行轉錄本(Transcript)拼接和Corset層次聚類，得到基因序列(Unigene)，最終得到的序列信息以FASTA格式儲存。拼接轉錄本與基因序列長度分布如圖1所示。

由圖1可知，Unigene共包含了19 213個基因片段，長度在2 000 bp以上的共有5 191個。

2.3 基因功能注釋

為得到更加全面的基因功能信息，本研究參照七大數據庫(Nr、Nt、Pfam、Go、Kog、Swiss-prot、Kegg)對最長的轉錄本進行了基因功能注釋。選擇與研究目的相近、注釋成功率高的5個數據庫繪制了維恩圖，見圖2。從圖2可以看出，能夠在5個數據庫中被正確注釋的基因共有1 042個，同時仍有大量未能被正確注釋的基因存在。注釋率較低的原因有以下兩點：①轉錄組數據里面有一部分是非編碼序列；②裂殖壺菌是研究較少的菌種，現有數據庫中該菌種及其相關菌種注釋較少。

圖2 基因功能注釋結果維恩圖

(1)Nr注釋：通過與Nr庫進行比對注釋，既能夠清楚地知道研究物種的基因序列與親緣關系較近的物種基因序列的相似性，也可以得到研究物種基因的功能信息。對比結果為：裂殖壺菌菌株與球孢白僵菌(Beauveriabassiana)、馬尾松毛蟲菌(Spizellomycespunctatus)、裂孢蛙糞霉菌(Basidiobolusmeristosporus)、根內球囊霉(Rhizophagusirregularis)、節水霉菌(Gonapodyaprolifera)以及其他菌種的基因序列相似性分別為14.7%、 5.5%、 5.4%、 3.2%、 3.1%及68.1%。

(2)Go是一套國際標準化的用于描述基因功能的分類系統。Go分為3個大類：生物過程(Biological Process，BP)、細胞組成(Cellular Component， CC )、分子功能(Molecular Function，MF)。Go注釋對基因編碼的產物進行深度解析，分別探究其參與的生物過程、具有的分子功能和所處的細胞環境。對基因進行Go注釋以后，將注釋成功的基因依照下一層級進行分類，結果如圖3所示。從圖3可知：歸類到生物過程的這類基因中，大量的基因參與細胞內過程類別、代謝過程和單組織過程類別；參與細胞組成的基因主要集中在細胞組成、膜組成和細胞器組成3個部分；分子功能類別中，具有結合和催化活性功能的基因較多，遠高于其他功能下的基因數目。

(3)Kegg數據庫：該數據庫是一個系統性分析基因產物功能以及其參與細胞中哪些代謝途徑的數據庫。對基因進行Kocc(Kegg ortholog)注釋后，可以根據其參與的Kegg代謝通路分為以下幾個類別：細胞過程(A，Cellular Processes)、環境信息處理(B，Environmental Information Processing)、遺傳信息處理(C，Genetic Information Processing)、代謝(D，Metabolism)、有機系統(E，Organismal Systems)，結果見圖4。從圖4可知，參與代謝過程的基因所占的比例比較高，主要集中在碳水化合物和氨基酸的代謝過程中。

圖3 Go分類圖

2.4 差異基因表達

2.4.1 基因差異表達分析

根據每個樣品基因表達水平的結果，用DESeq2[16]進行分析，檢測樣品之間的差異表達基因(Differentially Expressed Gene，DEG)，制作火山圖(略)，差異基因篩選閾值為padj<0.05且|log2FoldChange|>1。對于差異基因，如果基因的log2FoldChange>0，則認為該差異基因上調，反之，認為該差異基因下調。以-log10(padj)的數值大小來表征差異的顯著性。與對照組(原始菌株SR21)相比，實驗組(誘變菌株NA2)共有314個基因上調，3個基因下調，且上調部分基因上調的幅度極大。

繪制基因表達的維恩圖，見圖5。

圖5直觀地展現了2個實驗樣品之間獨有和共有的表達基因數目。兩組樣品共有14 278個共同表達基因，實驗組有4 025個特有表達基因，對照組有889個特有表達基因。

圖4 Kegg分類圖

注：以FPLM>0.3為基因表達的標準，每百萬碎片中來自某一基因每千堿基長度的碎片數目。

圖5 基因表達維恩圖

2.4.2 Kegg富集分析

在生物體中，許多基因之間相互協調，共同發揮作用，行使某項生理功能， 通過Pathway顯著性富集能夠明確這些差異表達基因具體參與的是何種生化代謝途徑和信號轉導途徑。富集結果如表4所示。

由表4可知，共有19個上調表達的基因被富集到17個代謝途徑中，主要包括氨基酸代謝、能量代謝等。同時可以看出，同一個基因可能參與幾個代謝途徑， 最終得到9個差異基因，分別是ECHS1、TPIA、SETD1、CHK2、SIR2、MSI、REA1、UBE2Z、PURL。

ECHS1(烯脂酰輔酶A水合酶短鏈1)是脂肪酸分解代謝β氧化通路的關鍵酶[20]，同時參與支鏈氨基酸氧化。從表4可知，ECHS1同時參與賴氨酸、β-谷氨酸、色氨酸、纈氨酸、亮氨酸和異亮氨酸等多種氨基酸的降解。氨基酸的降解可生成更多的乙酰-CoA，為脂質合成提供更多的前體物質。另外ECHS1還參與脂肪酸的延長和降解途徑，β氧化是一個脂肪酸改造的過程，根據機體自身的需要，一些長鏈脂肪酸會分解，釋放能量，然后合成長度適宜的脂肪酸。相比于飽和脂肪酸的β氧化，DHA的β氧化還需要異構酶和還原酶的參與[21]。ECHS1基因上調說明細胞內脂肪酸代謝活躍，脂肪酸合成速度太快，前體物質無法滿足，因此一些脂肪酸被分解以提供更多的能量和原材料供某些特殊脂肪酸的合成[22]，這可能是誘變菌株NA2高產DHA的原因。

表4 差異基因Kegg顯著性富集列表

TPIA(磷酸丙糖異構酶)是糖酵解途徑中的關鍵酶之一，在糖酵解和脂肪酸合成等途徑中具有非常重要的作用。TPIA催化磷酸二羥丙酮(DHAP)發生異構化形成3-磷酸甘油醛(G-3P)并用于后續代謝過程，在這個過程中NAD+轉化為NADH，為脂質合成提供能量。

SETD1是組蛋白賴氨酸甲基轉移酶SET1家族的成員，主要對賴氨酸和精氨酸的殘基進行組蛋白甲基化修飾。另外SETD1和CHK2都與DNA的修復有關。細胞周期檢測點CHK2基因上調可加強對G2/M期的檢查，防止帶有基因組 DNA 損傷的細胞進入有絲分裂，維持真核細胞基因組的穩定性。

SIR2是一組進化上高度保守的NAD+依賴的組蛋白去乙酰化酶[23]，控制細胞壽命、轉錄調控、能量代謝、抗氧化應激。研究表明該酶參與染色質沉默，染色質沉默對維持染色體結構穩定和基因調控起重要作用[24]。MSI是一種RNA結合蛋白，參與細胞不對稱分裂及基因轉錄調控，抑制細胞凋亡和分化。REA1 是核糖體合成因子，與MSI共同作用實現了基因的高效表達。

UBE2Z：泛素介導的蛋白質水解過程，特異性強，效率高，需要消耗能量。待降解的蛋白質首先被泛素標記，然后進入細胞的蛋白酶復合的活性結合位點，進而被降解成短肽。PURL：調控嘌呤代謝途徑關鍵酶基因表達水平，影響經AMP→ADP→ATP代謝途徑，促進能量的產生。UBE2Z和PURL基因上調，說明細胞內代謝活躍，表明細胞增殖迅速，這可能與生物量的提高有關。

2.5 RT-qPCR驗證

轉錄組測序分析的結果往往需要通過RT-qPCR進一步驗證，從轉錄組測序結果中篩選出9個基因進行RT-qPCR驗證，以裂殖壺菌18S rRNA為內參，設計引物如表5所示。

表5 引物設計序列

續表5

基因名稱引物序列(5'-3')SIR2F:GATAA CGGACCTAC TGCCATGR:CTTTCCTCAGCACGATCTACGMSIF:ACCTCCCTACAATCAGTCTCGR:TCCAAATTCCAGAGCTGTCAGREA1F:TGGTAATCGGTATCTGTGAGT-TGR:ACCAGTCGTTTCATTGTCCTCUBE2ZF:TCGCCAAGTACCAAGTAAAGAGR:CCCTGAGTATCCCAATTGTCCPURLF:TGGTAATCGGTATCTGTGAGT-TGR:ACCAGTCGTTTCATTGTCCTC

將RT-qPCR的結果整理后得到基因相對表達差異倍數如圖6所示。

圖6 RT-qPCR基因相對表達差異倍數

為了驗證轉錄組測序結果的準確性，以RT-qPCR中-ΔΔCt值與轉錄組分析的結果log2ratio對比，結果如圖7所示。由圖7可知，差異基因的表達量倍數上有所差異，但是在變化趨勢上是一致的。說明轉錄組測序的結果是可信的。

圖7 相關酶基因RNA-seq和RT-qPCR基因相對表達量對比

RT-qPCR的結果進一步證明了脂肪酸代謝、氨基酸代謝途徑中一些關鍵基因的突變，引起轉錄水平的變化。維持基因組穩定性和實現高水平轉錄等的相關基因上調共同促進了DHA的積累。以上結果表明ARTP這種新型的誘變方式是簡單有效的，產生了較多的突變位點，同時利用丙二酸平板成功篩選獲得高產菌株，在維持葡萄糖濃度恒定，間歇補加氮源的補料分批發酵策略下提高了DHA的產量。

3 結 論

采用轉錄組分析的方法揭示了誘變菌NA2中DHA產量提高的分子機理：與原始菌株相比，誘變菌NA2共有314個基因上調，3個基因下調；上調的基因主要參與氨基酸代謝和能量代謝，為多不飽和脂肪酸的積累提供更多的乙酰-CoA和NADPH。