新疆兩種亞麻籽轉錄組分析及籽油香氣差異基因挖掘

楊佳瑋，李歡康，林雨霏，劉文玉，魏長慶

（石河子大學食品學院，新疆植物藥資源利用教育部重點實驗室，新疆 石河子 832000）

亞麻（Linum usitatissimumL.）是我國重要的經濟農作物[1]。亞麻籽是亞麻的種子，富含不飽和脂肪酸、木酚素、酚酸、植物甾醇、生育酚、類黃酮等多種營養成分[2-4]，對心血管疾病[5-6]、糖尿病[7]、高血壓[8]、高血脂[9]、抑制癌癥[10-12]等均具有有益作用。亞麻全籽含油量約達40%[13]，其籽油因其特殊濃郁的香氣和豐富的營養成分，備受消費者青睞。

轉錄組學是一門研究細胞中基因轉錄情況以及轉錄調控規律的學科，細胞中的轉錄組因其細胞的生長發育、細胞所處的環境而異[14-15]。轉錄組測序技術是基于Illumina HiSeq等測序平臺對特定組織或細胞某個時期轉錄出來的所有mRNA進行測序[16]，具有通量大、低成本、高分辨率、高靈敏度、檢測范圍廣、不需克隆且操作簡單等優點[17]，是研究生物過程與分子機制的重要手段。近年來，轉錄組測序技術在農作物基因挖掘方面的運用很廣泛，以Illumina HiSeq 2000測序平臺為例，此項技術已在辣椒[18]、核桃[19]、花生籽仁[20]、大豆[21-22]、火龍果[23]中成功應用，為研究農作物基因功能提供大量的轉錄本信息，成為模式和非模式植物功能基因挖掘[24-26]的有效工具。

新疆作為中國亞麻的六大產區之一，培育出了很多品種，包括伊亞-1號、伊亞-3號（YY-3）、伊亞-4號、隴亞-7號、隴亞-8號、天亞-6號（TY-6）等系列品種，各個品種亞麻籽油香氣具有明顯差異。課題組前期對YY-3和TY-6品種的亞麻籽油揮發性香氣進行定性定量測定，發現2個品種揮發性香氣差異較大，其中醛類揮發性化合物含量及種類差異明顯，且TY-6亞麻籽油的醛含量（66.15%）顯著高于YY-3樣品（52.27%）[27]。基于此，本研究以亞麻籽為研究對象，選取新疆2個特色亞麻籽品種YY-3和TY-6進行轉錄組測序和分析，旨在從分子角度挖掘亞麻籽油香氣差異形成關鍵基因，并為新疆亞麻品種育種及籽油的香氣改良提供理論依據和技術支撐。

1 材料與方法

1.1 材料與試劑

2個新疆亞麻籽品種YY-3、TY-6均來源于新疆伊犁雅琪娜公司。于亞麻籽成熟之際采集，用冰袋冷藏運回石河子大學食品學院，經液氮速凍后存于-80 ℃冰箱，用于RNA提取。

1.2 儀器與設備

TG16-WS臺式高速離心機 湘儀離心機儀器有限公司；2100生物分析儀 美國安捷倫公司；NanoDrop 2000超微量分光光度計 美國熱電公司。

1.3 方法

1.3.1 RNA的提取與樣本檢測

將采集的樣品使用植物總RNA提取試劑盒提取亞麻籽總RNA，并通過NanoDrop 2000超微量分光光度計檢測RNA的濃度和純度，采用瓊脂糖凝膠電泳分析RNA質量，采用2100生物分析儀檢測RNA的完整性，以保證后續實驗結果的可靠性。

1.3.2 序列測序組裝

將質檢合格后的RNA利用高通量測序平臺（Illumina HiSeq 2000）進行測序。將下機原始數據（raw reads）進行堿基質量、堿基錯誤率和堿基含量統計評估后，進一步進行質量控制，去掉reads中的接頭序列、低質量的reads，得到clean reads進行后續分析[28]，統計其數量使用Q30等指標對數據進行質量評估。用Trinity軟件對所有clean reads進行從頭組裝得到大量轉錄本，取每條基因中最完整的轉錄本作為unigene，并對其進行優化過濾，且對組裝結果進行質量評估。

1.3.3 unigene的注釋、功能分類和生物學通路分析

將2個亞麻籽品種轉錄組測序獲得的所有unigene利用BLAST軟件將基因注釋到非冗余蛋白（Non-Redundant Protein，NR）、直系同源群集（Cluster of Orthologous Group，COG）、基因本體（Gene Ontology，GO）和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）數據庫，并對獲得的注釋情況及功能信息進行統計。

1.3.4 差異表達基因（differentially expressed genes，DEGs）的篩選及分析

使用TPM（transcripts per million）法估算基因表達水平，然后根據模型進行假設檢驗概率（P值）計算，最后進行多重假設檢驗校正。按照|log2（YY-3/TY-6）|＞1（YY-3/TY-6表示差異倍數，下同）且錯誤發現率（false discovery rate，FDR）＜0.05的原則篩選DEGs。利用比對軟件將篩選出的DEGs進行GO分類以及KEGG通路分析。

1.4 數據處理

兩個品種均進行3 次獨立生物學重復實驗。利用Excel 2007、Origin 8、SPSS 16.0等軟件對實驗數據進行統計分析及圖像處理。FastQC軟件進行測序數據質量控制，CutadaptVersion 1.2.1去除3’端的接頭序列，Trinity軟件進行轉錄組de novo組裝，DEGSeq分析表達差異unigene。

2 結果與分析

2.1 測序數據及比對結果統計

對YY-3、TY-6亞麻籽進行轉錄組測序共得到68 221 320 條raw reads，對其進一步過濾得到61 995 756 條高質量的clean reads（表1）。其中Q20分別為98.80%和98.79%，Q30分別為96.73%和96.68%，Q20和Q30值越高說明測序錯誤越小，GC含量占總堿基數的52%和53%，說明這兩種新疆特色亞麻籽轉錄組的測序質量較高，可用于de novo組裝和后續分析。

表1 兩種亞麻籽轉錄組測序數據統計Table 1 Data statistics of transcriptome sequencing of two flaxseed cultivars

將測序reads與參考數列進行比對得到各個樣品的匹配率（表2）。兩種亞麻籽比對到的reads占比均在90%以上，說明參考基因組選擇較為合適，且數據結果適用于后續基因的功能注釋和基因挖掘。

表2 兩種亞麻籽序列比對統計Table 2 Statistics of sequence alignment of two flaxseed cultivars%

2.2 基因注釋結果分析

利用String Tie軟件結合de novo進行轉錄本重構，將組裝好的轉錄片段與參考基因進行比較合并，將2個亞麻籽品種轉錄組測序獲得的所有unigene通過BLAST軟件與NR、COG、GO、KEGG數據庫進行比對，對所有注釋上的unigene進行分析。由表3可知，獲得注釋的unigene共有33 218 條。其中GO數據庫注釋比例最高，共注釋了17 542 條（52.81%）unigene；COG和NR數據庫其次，分別占unigene總數的51.99%和47.63%，KEGG數據可注釋比例最低，僅有10 713 條（32.25%）unigene被注釋。

表3 亞麻籽轉錄組功能注釋Table 3 Functional annotation of flaxseed transcriptome

2.3 unigene的NR注釋

將33 218個unigene注釋到NR數據庫，共有15 821個unigene注釋到NR數據庫，占總unigene的47.63%。

從匹配序列的相似度分布（圖1a）可見，有28.61%序列的相似度大于80%，26.43%序列的相似度在70%～80%之間，有41.79%序列相似度在40%～70%之間，僅有3.17%的相似度在40%以下。從E-value分布（圖1b）可見，有62.01%的E-value分布在0～10-100之間；有9.12%的E-value分布在10-80～10-60之間。其中E-value越小說明匹配結果的可靠性越高。從注釋匹配的物種分布（圖1c）可見，排在前5 位的物種分別為麻風樹（Jatropha curcas，27.97%）、蓖麻（Ricinuscommunis，13.67%）、毛果楊（Populustrichocarpa，11.20%）、胡楊（Populuseuphratica，10.85%）和可可樹（Theobroma cacao，4.23%）。其中注釋到麻風樹的基因數量最多，說明亞麻籽轉錄組測序組裝的unigene與麻風樹的相似性最高。由于亞麻籽的全部基因組信息尚未被報道，所以本研究中部分unigene在NR數據庫中無法和已知基因匹配，這些unigene可能包含了亞麻籽與其他物種不同的、自身所特有的基因序列。

圖1 亞麻籽unigene的NR注釋分類Fig. 1 Annotation classification of flaxseed unigene via NR

2.4 unigene的COG注釋

將unigene注釋到COG數據庫，根據結果對unigene的功能進行分類并統計注釋unigene數目（圖2）。結果表明，有17 269 個unigene注釋到COG數據庫，占總注釋到unigene的51.99%，此數據庫將基因注釋到3 大類，包括細胞過程和信號、信息存儲與處理、代謝。22個各功能類的基因注釋數目有明顯差異，其中翻譯后修飾、蛋白質折疊和分子伴侶最多，占15.48%；其次是轉錄占10.90%，細胞內運輸、分泌和囊泡運輸占比為9.91%，信號轉導機制占比8.87%，翻譯、核糖體結構和生物合成占比8.74%，復制、重組和修復占比6.77%，其余均在5%以下。細胞運動和核結構類別中注釋基因最少，均為0.14%。

圖2 亞麻籽unigene的COG注釋分類Fig. 2 Annotation classification of flaxseed unigene via COG

2.5 DEGs分析

2.5.1 基因表達水平對比

圖3TPM盒形圖Fig. 3 TPM box chart

為進一步分析基因表達差異情況，本研究引入TPM計算基因表達量。圖3為2個樣品的TPM盒形圖，可以看出兩組亞麻籽品種基因表達之間的差異。兩組數據分散程度基本相同，YY-3的中值較大，說明相較TY-6，它的整體表達量較高。

2.5.2 DEGs篩選及聚類分析

根據條件|log2（YY-3/TY-6）|＞1且FDR＜0.05，YY-3亞麻籽樣品相比較TY-6亞麻籽樣品共得到1 658個DEGs，其中上調基因825個，下調基因833個（圖4）。紅色點的縱坐標值均較大，說明這些上調基因的表達量變化差異較為顯著。

圖4 兩種亞麻籽DEGs火山圖Fig. 4 Volcano diagram of DEGs between two flaxseed cultivars

2.5.3 DEGs GO富集分析

對亞麻籽unigene采用interproscan-5.8-49進行注釋。兩個不同品種亞麻籽篩選出的DEGs進行GO富集分析，將其分為3 大類33個分支（圖5）。第一層次的GO Term由細胞組分、生物過程和分子功能三大類構成。其中，歸屬生物過程的GO條目最多（13個）。生物過程中兩組亞麻籽DEGs主要集中在代謝過程、細胞過程和單有機體過程。在細胞組分中，以細胞、細胞部分和膜3個亞類為主。在分子功能中注釋信息最多的兩個功能亞類是結合和催化活性。從圖5可以看出，由于品種不同，導致催化活性等均不同，所以細胞和細胞部件、膜等具有差異，從而導致了兩組亞麻籽代謝過程、細胞過程和單有機體過程的差異性。

圖5 DEGs的GO注釋分類Fig. 5 Annotation classification of differentially expressed genes via GO

2.5.4 DEGs KEGG富集分析

KEGG是基于分子水平信息，特別是大型分子數據集合而生成的基因組測序數據庫和其他高通量實驗得出的數據庫資源，是一個有關通路的主要公共數據庫[29]，可用于進一步研究基因的復雜行為。根據KEGG數據庫對兩個不同品種新疆亞麻籽進行功能分類和通路注釋。380 個DEGs被注釋到108 條代謝途徑，對DEGs進行KEGG通路富集分析，獲得兩種亞麻籽DEGs富集的代謝途徑（圖6）。

圖6 DEGs KEGG富集散點圖Fig. 6 KEGG enrichment scatter plot of DEGs

利用散點圖形式表現KEGG富集程度，計數比為某功能類別中差異表達數與所有功能類別中注釋到的DEGs數的比值，取值范圍為（0，1）。結果表明，兩種亞麻籽DEGs主要富集到脂肪酸代謝，不飽和脂肪酸合成，甘氨酸、絲氨酸和蘇氨酸代謝，抗壞血酸和藻酸鹽代謝，葉酸一碳庫代謝等。

2.5.5 脂質代謝途徑相關DEGs分析

由KEGG富集程度分布圖可知，富集程度最高的前2個KEGG條目為脂肪酸代謝和不飽和脂肪酸合成。其中在脂肪酸代謝中注釋到17個DEGs，在不飽和脂肪酸代謝中注釋到了8個DEGs，對其分別進行聚類分析（圖7、8）。

圖7 脂肪酸代謝相關DEGs熱圖Fig. 7 Heatmap of DEGs related to fatty acid metabolism

圖8 不飽和脂肪酸代謝相關DEGs熱圖Fig. 8 Heatmap of DEGs related to unsaturated fatty acid metabolism

由圖7、8可知，兩個品種在脂肪酸代謝通路上注釋到的17個DEGs中，有12個基因在TY-6中顯著表達，分別為：MSTRG.8561.1、Lus10040173.BGIv1.0、Lus10036184.BGIv1.0、Lus10012007.BGIv1.0、Lus10029283.BGIv1.0、Lus10006502.BGIv1.0、MSTRG.12234.3、Lus10020160.BGIv1.0、MSTRG.3862.1、Lus10038321.BGIv1.0、MSTRG.117.1和MSTRG.22853.1。這12個基因中有6個基因同時注釋到了不飽和脂肪酸通路上，分別是：Lus10012007.BGIv1.0、Lus10029283.BGIv1.0、Lus10036184.BGIv1.0、MSTRG.117.1、Lus10038321.BGIv1.0和MSTRG.22853.1。由圖7可知，5個基因在YY-3中顯著表達，分別為Lus10020908.BGIv1.0、MSTRG.236.1、Lus10033072.BGIv1.0、Lus10027811.BGIv1.0和Lus10004186.BGIv1.0，其中MSTRG.236.1和Lus10027811.BGIv1.0這2個基因同時注釋到了不飽和脂肪酸通路上，參與調解不飽和脂肪酸的生物合成。

王歡歡等[30]在對番茄果實呈香組分及其代謝途徑的研究中發現，番茄果實中的11 種香氣組分分別來自于以脂肪酸、氨基酸和類胡蘿卜素為前體物質的代謝途徑。在脂肪酸代謝途徑中，不飽和脂肪酸在脂肪氧合酶（lipoxygenase，LOX）催化下生成的氫過氧化物通過氫過氧化物酶（hydroperoxide lyase，HPL）的作用，可轉化成己醛、順-3-己烯醛等C6醛類化合物。Olias等[31]對初榨橄欖油中揮發性香氣成分的研究表明，當植物受到創傷時，LOX變得活躍，導致不飽和脂肪酸形成以己醛為代表的C6醛類化合物。亞麻籽破碎過程中，細胞組織被破壞，亞油酸或亞麻酸在LOX的催化下會生成具有共軛二烯結構的氫過氧化物，進而HPL導致O—O鍵的斷裂，通過進一步裂解迅速生成了大量醛類揮發性物質。杜樂[32]對冷榨亞麻籽油中不飽和脂肪酸氧化產物的研究表明，不飽和脂肪酸自動氧化形成醛類羰基化合物從而影響油的風味和質量，同樣說明了不飽和脂肪酸是揮發性香氣形成的重要前體物質。這些研究均揭示了亞麻籽品種中不飽和脂肪酸含量及種類的差異性與其籽油醛類特征香氣的生成具有重要相關性。此研究中對兩個新疆亞麻籽品種上游轉錄組差異代謝途徑及相關基因的發掘，可對下游兩個品種籽油的香氣差異的進一步探究提供理論依據。

3 結 論

采用Illumina HiSeq測序平臺對兩種新疆特色亞麻籽進行轉錄組測序，兩組數據共獲得68 221 320 條raw reads，過濾得到61 995 756 條clean reads。采用生物信息學軟件對測序結果組裝得到33 218個unigene，將其與NR數據庫進行比對發現，亞麻籽轉錄組測序組裝的unigene與麻風樹的相似性最高。利用COG數據庫進行比對發現可將17 269 個unigene注釋到細胞過程和信號、信息存儲與處理、代謝3 大類22個分支，其中翻譯后修飾、蛋白質折疊和分子伴侶是富集最多的條目，占比15.48%。

按照｜log2（YY-3/TY-6）｜＞1且FDR＜0.05的原則，共篩選出1 658個顯著性DEGs，其中825個基因上調，833個基因下調表達。通過DEGs GO富集分析可知，可將其分為細胞組分、分子功能、生物過程3 大類33個分支，其中歸屬生物過程的GO條目最多（13個）；380 個DEGs被注釋KEGG數據庫的108 條代謝途徑中，兩種亞麻籽DEGs KEGG代謝通路中顯著富集的通路有脂肪酸代謝，不飽和脂肪酸合成，甘氨酸、絲氨酸和蘇氨酸代謝，抗壞血酸和藻酸鹽代謝，葉酸一碳庫代謝等。通過KEGG數據庫對脂質代謝途徑相關DEGs分析可知，在脂肪酸代謝中注釋到17個DEGs，其中有12個基因在TY-6中顯著表達，5個基因在YY-3中顯著表達。在不飽和脂肪酸代謝中注釋到了8個DEGs，其中在TY-6顯著表達的DEGs有6 條，2 條基因在YY-3中顯著表達。研究結果為深入開展不同品種亞麻籽脂肪酸代謝過程研究及籽油香氣改良提供科學依據。