不同發育時期藍果忍冬果實的轉錄組分析

孫燕 王金剛 臧丹丹 趙恒田 劉淑華

（1.東北農業大學園藝園林學院，哈爾濱 150030；2.中國科學院東北地理與農業生態研究所黑土區農業生態重點實驗室，哈爾濱 150081）

藍果忍冬（Lonicera caeruleaL.）為忍冬科忍冬屬植物，原產于北半球溫帶地區，現廣泛分布于我國各地，是我國重要的野生果樹資源［1-3］。藍果忍冬果實出汁率高，有香氣，營養價值豐富，具有調節血壓、減緩衰老等作用［4-5］，因其富含花青素、氨基酸、維生素、類黃酮等多種活性物質，被譽為“飲料之王”［6］?！锌扑{1 號’是我國第一個具有自主知識產權和國家備案的藍果忍冬新品種，植株健壯，豐產性強，果實為長橢圓形，口感酸甜無苦味，是優良的鮮果品種，適宜東北寒區廣泛種植［7］。

花青素是植物體內一類重要的次生代謝物質，屬于類黃酮類化合物，在賦予植物色彩時起主要作用［8-9］。在果實生長發育過程中，果皮中色素會發生一系列變化，隨著果實的成熟，葉綠素在水解酶的作用下逐漸分解消失，類胡蘿卜素略有增加；花青素在細胞質中合成，以苯丙氨酸為直接前體，由一系列結構基因編碼的合成酶催化形成，并經過一系列的運輸后貯存于液泡中；漿果的大小、色澤、硬度均會隨著果實的成熟發生顯著變化［10-11］。近年來，關于藍果忍冬果實花青素的研究受到廣泛關注，Zhou 等［12］對藍果忍冬品種‘蓓蕾’果實中的花青素提取并純化，發現提取的花青素ABL-2 的主要成分是花青素-3-葡萄糖苷、花青素-3，5-二葡萄糖苷、花青素-3-蘆丁糖苷和芍藥-3-葡萄糖苷，且能夠顯著抑制細胞生長，誘導DNA 損傷，殺死腫瘤細胞，改善小鼠的生存狀態；Paereckait 等［13］利用分光光度法對11 種藍果忍冬果實中的酚類化合物和花青素含量等進行提取研究并比較不同藍果忍冬品種的抗菌活性。然而，目前有關藍果忍冬基因信息尚不明確，對藍果忍冬果實轉錄組的相關研究少有報道，嚴重限制了藍果忍冬的開發和利用。

本研究以藍果忍冬‘中科藍1 號’果實為試驗材料，使用高通量測序平臺對綠熟期果實和成熟期藍果忍冬果實進行轉錄組測序分析，通過組裝建立藍果忍冬轉錄組數據庫，探究藍果忍冬果實不同發育時期的功能基因和調控網絡，旨在為后續了解藍果忍冬果實發育提供理論參考。

1 材料與方法

1.1 材料

試驗材料為藍果忍冬‘中科藍1 號’品種，由中國科學院東北地理與農業生態研究所培育，分別采取綠熟期（花后28 d，果實成長到體積不再發生變化，果皮顏色為綠色）和成熟期（花后38 d，果皮顏色全部變為藍紫色，果實未軟化）兩個時期的果實（果皮和果肉），各取3 個生物學重復，用液氮速凍后保存在-80℃冰箱中用于后續轉錄組測序和花色苷含量的測定（圖1）。

圖1 不同發育時期的‘中科藍1 號’果實Fig.1 Fruits of ‘Zhongkelan No.1’ at different developmental stages

1.2 方法

1.2.1 不同時期的花色苷總量測定 參照高效液相色譜法（high performance liquid chromatography，HPLC）［14］，利用安捷倫液相色譜系統（Agilent Technologies 1260 Infinity Ⅱ Prime，US）檢測花青素種類及含量，流動相A 為含3% 甲酸的雙蒸水，流動相B 為100% 乙腈。洗脫程序為0 min 90%A，10% B；16 min 86% A，14% B；20 min 90% A，10% B；流速為1 mL/min，柱溫為30℃，檢測波長為520 nm，每個樣品設3 個重復；分別以矢車菊素-3，5-雙葡萄糖苷、飛燕草素-3-O-葡萄糖苷、飛燕草素-3-O-蕓香糖苷、矢車菊素-3-O-葡萄糖苷、矢車菊素-3-O-蕓香糖苷、天竺葵素-3-O-葡萄糖苷、芍藥素-3-O-葡萄糖苷、芍藥素-3-O-蕓香糖苷標準品繪制標準曲線，含量由峰面積計算方法得出，單位為mg/100 g FW。

1.2.2 RNA 提取 采用CTAB 法提取藍果忍冬不同發育時期的總RNA，經PrimeScript One Step RT-PCR Kit Ver.2（TaKaRa）反轉錄為cDNA 并將cDNA 稀釋10 倍用于后續RT-qPCR 實驗。樣品總RNA 經檢測合格后送至北京百邁客生物技術有限公司運用高通量平臺進行RNA-Seq 測序。轉錄組樣品編碼為G：“綠熟期”，B：“成熟期”。

1.2.3 轉錄組測序與分析 為了獲得高質量的轉錄組數據，采用illumina novaseq 6000 平臺進行雙末端測序，截除掉原始序列接頭和引物序列，對低質量的序列進行過濾處理，并利用Trinity 軟件對生成的高質量測序數據進行組裝，得到轉錄本序列Unigene；利用BLAST 軟件在NR、GO、Swiss-Prot、COG、KEGG 和eggNOG4.5 數據庫中進行比較分析并預測Unigene 的氨基酸序列，使用HMMER軟件將Unigene 序列在Pfam 數據庫比對從而獲得Unigene 的功能信息，并使用 TBtools 軟件繪制熱圖。

1.2.4 差異表達基因的篩選和富集分析 用DESeq2軟件對藍果忍冬綠熟期和成熟期果實的轉錄組數據進行組間差異表達分析，篩選標準為FDR（False Discovery Rate）＜0.01 且FC（Fold Change）≥2，并利用Goatools 軟件對篩選出的差異表達基因進行COG 注釋分析、GO 功能富集分析和KEGG 富集分析。

1.2.5 差異表達基因的RT-qPCR 驗證 根據藍果忍冬轉錄組數據隨機篩選出10 個差異基因，在NCBI數據庫中比對并設計定量引物（表1）。選用LcACT1（MT344113）、LcTUB1（MT344114）作為內參基因，利用Agilent AriaMx 實時熒光定量PCR 儀（Agilent Stratagene，USA）對其進行實時熒光定量PCR 反應，程序為94℃預變性30 s，94℃變性15 s，57℃退火30 s，72℃延伸35 s，共44 個循環結束，65-95℃熔解曲線讀取，用2-ΔΔCT算法計算處理。

表1 轉錄組數據RT-qPCR 驗證引物序列Table 1 Primer sequences used in RT-qPCR validation of transcriptome data

2 結果

2.1 藍果忍冬‘中科’兩種時期的花色苷含量變化

利用高效液相色譜（high performance liquid chromatography，HPLC）檢測技術和紫外檢測器對520 nm 波長下進行掃描，檢測到的峰圖如圖2所示，可以看出，藍果忍冬‘中科’綠熟期果實未檢測出花色苷成分（圖2-a），成熟期果實中檢測到的花青素種類共有8 種（A1-A8）（圖2-b），總花色苷含量為598.55 mg/100 g FW。其中，矢車菊素-3-O-葡萄糖苷含量最高，占總花色苷的86.28%，是藍果忍冬‘中科’果實中主要的花色苷（表2），說明不同發育階段藍果忍冬‘中科’果實中花色苷類化合物組成和含量有較大差異。

表2 ‘中科藍1 號’成熟期果實中花色苷含量Table 2 Anthocyanins contents in the ripening stage fruits of ‘Zhongkelan No.1’

圖2 ‘中科藍1 號’不同發育時期果實高效液相色譜Fig.2 HPLC of ‘Zhongkelan No.1’ fruit at different developmental stages

2.2 轉錄組測序數據統計

選擇藍果忍冬綠熟期果實和成熟期果實各3 次生物學重復，共6 個樣品，轉錄組測序共得到40.32 Gb Clean Data，每個樣品Clean Data 均大于6.07 Gb，GC 含量為46.32%-49.38%，Q30 均值在93.89%及以上，表明測序數據較好，可用于后續分析（表3）。

表3 測序數據統計結果Table 3 Statistics of sequencing data

6 個藍果忍冬果實的轉錄組數據經過質量分析處理后共得到236 033 條轉錄本序列，平均長度為1 371.99 bp，N50 長2 186 bp；組裝后獲得81 777個單基因序列，且長度大于1 kb 的Unigene 共有17 339 個（表4），對所獲得的單基因序列與NR、Swiss-Prot、KEGG、COG、KOG、GO 和Pfam 數據庫進行功能注釋比對，注釋成功的Unigene 共有28 992 條，占全部Unigene 的35.45%（表5）。

表4 組裝結果統計Table 4 Assembly result statistics

表5 單基因功能注釋統計Table 5 Single-gene function annotation statistics

2.3 不同發育時期差異表達基因篩選結果

根據評估標準FDR（False Discovery Rate）＜0.01 且FC（Fold Change）≥2，共篩選出3 247 個差異表達基因，其中有1 642 個基因在上調表達，1 605 個基因在下調表達，從圖3可以清楚地看出，在兩個生長時期藍果忍冬果實中基因表達水平存在較大差異。

圖3 差異表達基因的火山圖Fig.3 Volcano map of differentially expressed genes

2.4 不同發育時期差異表達基因COG功能分類

利用COG 數據庫對藍果忍冬果實兩個發育時期的差異表達基因進行對比分析（圖4）發現，注釋到COG 通路中的差異表達基因有1 081 個，分布在23 條途徑中，占全部差異基因的37.88%，其中，碳水化合物運輸和代謝（carbohydrate transport and metabolism）是最大通路，所占數量最多，其次是一般功能預測（general function prediction only），差異基因在脂質運輸和代謝（lipid transport and metabolism）通路和次級代謝物生物合成、運輸和分解代謝（secondary metabolites biosynthesis，transport and catabolism）通路占比也較高，說明在果實成熟發育過程中代謝物質變化較大。

圖4 不同發育時期差異表達基因COG 注釋分類統計圖Fig.4 COG annotation classification of differentially expressed genes at different developmental stages

2.5 不同發育時期差異表達基因GO功能富集結果

對藍果忍冬果實兩個發育時期的差異表達基因進行對比分析，發現注釋到GO 通路中的差異表達基因有1 908 個，分別映射到51 個代謝通路中（圖5），在生物過程（biological process）中，主要富集在代謝過程（metabolic process）和細胞過程（cellular process），其次是單個有機體過程（single-organism process）；在細胞組分（cellula component）通路中，膜部分（membrane）占比最多、其次是細胞（cell）和細胞部分（cell part）；在分子功能（molecular function）通路中，催化活性（catalytic activity）和結合功能（binding）較其他通路顯著富集?；ㄇ嗨厥侵参锎紊x產物，糖苷衍生物。推測可能歸類于代謝過程（metabolic process）和催化活性（catalytic activity）的Unigenes 中，可能有參與藍果忍冬花青素生物合成積累差異相關的良好候選基因。

圖5 不同發育時期差異表達基因GO 功能富集結果Fig.5 Gene ontology annotation of differentially expressed genes at different developmental stages

2.6 不同發育時期差異表達基因KEGG功能富集結果

以q-value＜0.05 作為篩選差異表達基因在KEGG數據庫通路的標準，由圖6可知，2 854 個差異基因中共有1 034 個DEGs 注釋到KEGG 通路中，其中植物激素信號轉導（plant hormone signal transduction）途徑是最大通路，有55 個，占比10.07%；其次是氨基酸生物合成（biosynthesis of amino acid）和碳代謝（carbon metabolism）途徑，均為46 個，占比為8.42%；此外，淀粉和蔗糖代謝（starch and sucrose metabolism）通路與果實的發育也顯著富集，有37 個，占比6.78%。

圖6 不同發育時期差異表達基因KEGG 分類圖Fig.6 KEGG classification of differentially expressed gene at different developmental stages

花青素是類黃酮化合物，且與苯丙氨酸代謝途徑有關。對KEGG 代謝通路進行富集顯著性篩選，得出富集最顯著的20 條代謝通路（圖7），其中，苯丙氨酸，酪氨酸和色氨酸的生物合成（phenylalanine，tyrosine and tryptophan biosynthesis）、類黃酮生物合成（flavonoid biosynthesis）、半胱氨酸和蛋氨酸代謝（cysteine and methionine metabolism）這3 條代謝途徑被顯著富集，說明可能存在大量與藍果忍冬花青素生物合成積累相關的差異基因。

圖7 不同發育時期差異表達基因KEGG 功能富集結果Fig.7 KEGG analysis of differentially expressed genes at different developmental stages

2.7 藍果忍冬轉錄組測序數據的RT-qPCR驗證

為了驗證RNA-seq 的結果，對測序結果隨機挑選了10 個差異表達基因進行RT-qPCR 驗證，分析RT-qPCR 與RNA-seq 結果之間的相關性，每個基因進行3 次生物學重復。結果顯示（圖8），隨機挑選的10 個差異表達基因的表達量變化與轉錄組基因的變化完全一致，說明轉錄組的測序數據和RT-qPCR的數據之間具有良好的一致性，表明本研究RNAseq 數據具有較高的可靠性。

圖8 轉錄組測序數據的RT-qPCR 驗證Fig.8 Validation of transcriptome sequencing data by RT-qPCR

3 討論

近年來，隨著分子生物學技術的發展，新一代高通量測序技術已被廣泛應用于生物分析和醫學研究等領域，因其具有成本低，測序速度快，能夠估計不同發育時期和不同組織中的整體基因的表達量等特點，是當前無參考基因組的非模式生物的經濟有效的選擇［15-18］。目前，已有研究表明在藍果忍冬果實或種子中分離和鑒定出酚類物質、抗氧化成分和揮發性物質等多種生物活性成分，包括黃酮類化合物、花青素、脂肪酸、甘油三酯和總甾醇等化合物，并證實花青素-3-葡萄糖苷為其主要的抗氧化成分［19-22］。本研究基于RNA-seq 測序技術，對‘中科藍1 號’藍果忍冬兩個發育時期的果實進行轉錄組分析，序列經過組裝之后，共獲得81 777 條平均長度為763.35 nt 的Unigene，與前人報道的‘蓓蕾’品種果實發育階段的基因數目相近［23］，比馬尾松、枇杷、荔枝等木本植物組裝的Unigene 多［24-26］，結果表明，組裝的質量以及片段長度可以滿足轉錄組數據分析的需要。有28 992 條Unigene 得到注釋，占比35.45%，雖然比大葉女貞（50.20%）和藍莓（73.60%）的注釋率低，但與文冠果（36.15%）和紅掌（38.87%）相近［27-30］，造成注釋率低的原因可能是缺乏近緣物種的參考基因組信息，未得到注釋的基因也可能是藍果忍冬果實特異性新基因，有待進一步分析研究。

對比兩個測序樣本檢測到的差異表達基因發現，表達量上調的基因數目多于下調基因數目，這表明成熟期果實相對于綠熟期果實的樣本有更多的基因表達量。對這些基因進行KEGG 分析后發現有大量的色素形成相關途徑被顯著富集，包括與花青素代謝直接相關的苯丙氨酸、酪氨酸和色氨酸的生物合成途徑、類黃酮生物合成途徑和氨基酸生物合成途徑，表明藍果忍冬果實色澤調控的復雜性。GO 富集分析發現在果實綠熟期到轉色期果實成熟過程中，參與細胞過程和代謝的基因較多，同時涉及到的功能最多的基因為結合和催化功能，從分子水平說明了這一過程中代謝旺盛。藍果忍冬不同發育時期的差異表達基因在植物激素信號轉導、碳代謝、氨基酸生物合成、脂肪代謝等多條代謝途徑中富集，與類黃酮合成緊密相關，具體受哪些基因控制還需要做進一步信息學分析。

園藝產品品質形成的調控包括多個方面，在植物激素相關研究中，表明GA 對果實成熟軟化［31］、有機酸及糖分積累［32］、果實著色［33］、果實大小和形狀［34］等方面都有一定影響。比如，趙榮華等［35］研究表明GA 水平的調控可同時影響葡萄果實的形態和品質，另有研究發現生長素和赤霉素的相互作用在辣椒的生長發育過程中能平衡細胞分裂和細胞擴張［36］。在轉錄因子相關研究中，Medina-Puche等［37］研究表明YB9/FaMYB11、FabHLH3 和FaTTG1 在調控草莓花青苷合成過程中起關鍵作用，調控果皮花青苷的生物合成［38］。本研究發現藍果忍冬‘中科藍1 號’成熟期果實含有8 種花色苷成分，其中含量最高的花色苷為矢車菊素-3-O-葡萄糖苷，占總花色苷含量的86.28%，而綠熟期果實未檢測出花色苷成分，推測在果實成熟過程中存在調控花色苷含量的轉錄因子，還需進一步研究。

4 結論

本研究對‘中科藍1 號’藍果忍冬綠熟期和成熟期果實進行了轉錄組測序，共有40.32 Gb Clean Data，每個樣品Clean Data 均大于6.07 Gb，Q30 堿基百分比均在93.89%及以上，GC 含量為46.32%-49.38%；共篩選出3 247 個差異表達基因，包含1 642 個上調基因和1 605 個下調基因，為研究藍果忍冬果實轉色過程的的分子調控機制提供了參考。