紅肉火龍果成熟莖的全長轉錄組測序分析及蔗糖代謝相關基因分離

鄭乾明，王小柯，馬玉華

(貴州省農業科學院果樹科學研究所， 貴陽 550006)

【研究意義】紅肉火龍果(Hylocereuspolyrhizus)是仙人掌科(Cactaceae)量天尺屬(Hylocereus)多年生果樹作物，其葉退化為短刺，肉質莖是主要的源器官。植株具有耐旱的生理特點，已成為貴州喀斯特石漠化地區的特色經濟作物。紅肉火龍果果實含有豐富的可溶性糖、有機酸、維生素C、類黃酮和甜菜苷色素等風味物質和營養成分，深受消費者喜愛[1-6]。果實中積累的風味物質和營養成分均來自于源器官的光合作用合成，蔗糖作為最主要的光合產物，經維管束韌皮部運輸至果實中參與代謝和貯存[7]。研究紅肉火龍果植株源器官(成熟莖)的轉錄組模式，分離光合作用等重要代謝途徑的關鍵基因，有利于闡明紅肉火龍果植株蔗糖合成的分子機制，對改良果實品質，提高果實商品價值具有重要意義。【前人研究進展】紅肉火龍果是近年發展的新興水果，目前尚無基因組序列發布。當前僅有基于Illumina平臺的第二代轉錄組測序，對其果實或根開展系列測序和研究。例如，對紅肉和白肉火龍果果實發育多個時期的轉錄組測序和比較，獲得若干參與甜菜苷色素合成的基因[8]，以及1個潛在的候選MYB轉錄因子[9]；對紅肉火龍果根系在鹽脅迫下的轉錄組測序分析，獲得一系列上調和下調的差異表達基因[10]。針對紅肉火龍果成熟果實與成熟莖的轉錄組測序分析，獲得79 658個Unigene，平均長度為690 bp，約44.11%的轉錄本在公共數據庫中注釋[11]。同時，還獲得若干在成熟莖中表達，參與蔗糖代謝和轉運的Unigene，但均未獲得基因全長。因此，目前針對紅肉火龍果的轉錄組測序數據較少，僅有的轉錄組數據均由Illumina平臺的第二代測序獲得，其序列長度較短，難以獲得目標基因的全長序列。近年來，由Pacific Biosciences公司開發的PacBio測序是應用較為廣泛的第3代測序平臺。在此基礎上發展起來的全長轉錄組測序，具有長讀長，無PCR擴增，無需拼接，較高的隨機錯誤率等特點[12]。目前全長轉錄組測序主要用于構建全長轉錄組數據庫、轉錄本可變剪接分析、轉錄本結構研究和完善基因注釋等[13]。【本研究切入點】提取紅肉火龍果成熟莖的總RNA，進行基于PacBio Sequel平臺的全長轉錄組測序，獲得全長轉錄組數據庫。【擬解決的關鍵問題】獲得紅肉火龍果成熟莖的全長轉錄組序列，根據功能注釋和分類結果篩選參與蔗糖合成和降解等重要代謝途徑的候選功能基因，為后續開展基因功能研究奠定基礎。

1 材料與方法

1.1 材料

于2018年11—12月，在貴州省鎮寧縣的火龍果種植園，選擇紅肉火龍果品種“紫紅龍”。在成年結果植株上選取生長良好、無病蟲害的健壯成熟莖，剪下后立即置于冰盒帶回實驗室。用無水乙醇擦凈后晾干，然后切取莖組織液氮中速凍，置于-80 ℃保存備用。

1.2 核酸提取

利用Trizol試劑提取成熟莖組織的總RNA，使用瓊脂糖凝膠電泳、紫外光分光光度計和Bioanalyszer 2100(Agilent Technologies)檢測其質量與濃度。

1.3 全長轉錄組測序

后續文庫構建、測序和數據分析委托安諾優達基因科技有限公司開展。首先富集含有Poly A的mRNA，反轉錄合成第一鏈cDNA，然后PCR擴增合成cDNA。構建SMRTbell文庫并進行損傷修復和末端修復，片段兩端連接測序接頭，形成具有莖環結構的插入片段測序文庫。利用PacBio Sequel(Pacific Biosciences)測序儀測序，獲得原始數據。

原始數據過濾獲得Polymerase reads，去除接頭序列后獲得Subreads并進行質控，獲得高質量的插入片段序列。對來源于同一條Polymerase read的Subreads片段進行合并和糾錯，獲得環狀一致性序列(Circular Consensus Sequencing，CCS)。分析CCS序列，根據是否含有完整接頭和poly A序列篩選獲得全長非嵌合體(Full-length no chimera，FLNC)序列。使用ICA(isoform-level clustering algorithm)算法，對FLNC進行相同結構間的自我聚類和糾錯，產生去冗余Consensus序列。將未測完整的插入片段序列比對到Consensus序列進行糾錯，獲得準確度大于99%的高質量轉錄本。使用Blastn和Blastx程序在NCBI非冗余核酸和蛋白質數據庫進行檢索注釋。

1.4 序列分析

利用Trans Decoder Release v3.0.1鑒定轉錄本序列中的編碼序列(Coding sequence，CDS)。對獲得的高質量轉錄本序列，在NCBI核酸數據庫(Nucleotide Sequence Database，NT)、NCBI非冗余蛋白數據庫(Non-Redundant Protein Database，NR)、蛋白質真核直系同源簇數據庫(Eukaryotic Orthologous Groups，KOG)、基因本體論數據庫(Gene Ontology，GO)、蛋白質家族域數據庫(Protein Families Database，Pfam)、KEGG(Kyoto Encyclopedia of Genes and Genomes)數據庫進行基因功能注釋。通過與植物轉錄因子數據庫(Plant Transcription Factor Database，Plant TFDB)比對查找轉錄因子。序列一致性計算使用DNAStar軟件的MegAlign程序，序列拼接使用BioEdit軟件。

2 結果與分析

2.1 全長轉錄組測序

從表1看出，通過PacBio Sequel測序，獲得聚合酶Reads數量740 069條，堿基數量為15.8 Gb。聚合酶Reads平均長度為21 366 bp，N50長度為42 865 bp。過濾后獲得插入片段(Subreads)數量共8 483 057條，堿基數據量為15.1 Gb。其平均長度為1775 bp，N50長度為2200 bp。將來源于同一環狀分子的Subreads聚類，獲得環狀一致序列(CCS)序列共計481 124條，其堿基數據量為1.1 Gb。CCS平均長度為2378 bp，N50長度為2623 bp，平均測序次數為13.87次。

表1 紅肉火龍果成熟莖的全長轉錄組測序結果

進一步篩選獲得全長非嵌合體(FLNC)數量為333 226條，堿基數據量為0.7 Gb，平均長度為2071 bp，N50長度為2250 bp。將來源于同一轉錄本的FLNC聚類獲得Consensus序列，并將非全長序列對Consensus序列校正和去冗余，獲得轉錄本30 973條。大部分轉錄本長度分布于1.8 kb，平均長度為1846 bp(圖1)。其中準確度大于99%的高質量轉錄本數量為30 313條，平均長度為1829 bp。

圖1 紅肉火龍果成熟莖的全長轉錄組獲得的轉錄本長度分布

2.2 CDS預測

從表2可知，對30 313條高質量轉錄本進行CDS預測，獲得CDS數量為40 255條。其長度最短為297 bp，最長為5202 bp，平均長度為926 bp，N50長度為1179 bp。統計其長度分布表明，CDS長度在200～600 bp的數量為14 372條，占35.70%；長度在800～1200 bp的數量為15 253條，占37.89%；長度為1200～2000 bp的數量有8456條，占21.01%；長度大于2000 bp 的數量為2172條，占5.40%。

表2 預測CDS的數量或長度分布

2.3 序列注釋

從表3看出，將30 313條高質量轉錄本在公共數據庫中進行序列注釋，共計有13 939條轉錄本被注釋，占45.98%。在NCBI非冗余蛋白數據庫(NR)中注釋的轉錄本數量為13 789條，占所有轉錄本的45.49%。在NCBI核酸數據庫(NT)中注釋的轉錄本數量為10 777條，占所有轉錄本的35.55%。在Pfam數據庫中注釋的轉錄本數量為13 136條，占所有轉錄本的43.33%。在KEGG數據庫注釋的轉錄本數量為7582條，占所有轉錄本的25.01%。

從表4可知，統計注釋到的相關物種上的轉錄本數量依次為甜菜(Betavulgaris)、菠菜(Spinaciaoleracea)、葡萄(Vitisvinifera)、克里曼丁橘(Citrusclementina)、麻瘋樹(Jatrophacurcas)和蓮(Nelumbonucifera)，其數量分別為9371條(50.95%)、4579條(24.90%)、436條(2.37%)、258條(1.40%)、205條(1.11%)和204條(1.11%)。

在GO數據庫中共有13 059條轉錄本被注釋到生物學過程、細胞成分和分子功能等三類，占所有轉錄本的43.08%(表3)。生物學過程中注釋到24個類別，其中含有轉錄本最多的三類依次為細胞的過程(66.41%)、代謝過程(60.94%)和單一的生物過程(39.06%)。細胞成分中注釋到22個類別，其中含有轉錄本最多的三類依次為細胞組分(82.03%)、細胞器(54.69%)和細胞器組分(35.94%)。分子功能中注釋到19個類別，其中含有轉錄本最多的三類依次為結合(64.06%)、催化活性(51.56%)和核酸結合轉錄因子活性(6.25%)。

在KOG數據庫注釋的轉錄本數量為10 018條，占所有轉錄本的33.05%(表3)。從圖2可知，注釋為僅通用功能預測的轉錄本最多，數量為1701條，占16.98%；其次為翻譯后修飾、蛋白質周轉和伴侶蛋白，為1203條，占12.01%；其他較多的分別為信號轉導機制(901條，8.99%)，翻譯、核糖體結構和生物發生(656條，6.55%)，功能未知(580條，5.79%)，碳水化合物的運輸和新陳代謝(543條，5.42%)，轉錄(530條，5.29%)，RNA加工和修飾(526條，5.25%)，細胞內運輸、分泌和囊泡運輸(518條，5.17%)。

表3 紅肉火龍果成熟莖全長轉錄組測序獲得的轉錄本注釋

2.4 轉錄因子家族預測

與植物轉錄因子數據庫(Plant Transcription Factor Database，PlantTFDB)比對，共計有11 628條轉錄本注釋為轉錄因子。統計轉錄因子家族表明，共含有56個家族，其含有轉錄本數量較多的家族如表5所示。MYB_related家族的轉錄本最多，數量為1162條，占10%；其次為bHLH家族(basic/helix-loop-helix)，轉錄本數量為1066條，占9.17%。NAC(NAM、ATAF和CUC)家族、WRKY家族和FAR1家族含有的轉錄本數量依次為788、687和669條，分別占所有轉錄因子數量的6.78%、5.91%和5.75%。

表4 NR數據庫注釋的物種數量

A：RNA加工和修飾；B：染色質結構和動力學；C：能源生產和轉換；D：細胞周期控制、細胞分裂和染色體分割；E：氨基酸轉運和代謝；F：核苷酸轉運和代謝；G：碳水化合物的運輸和新陳代謝；H：輔酶轉運和代謝；I：脂質運輸和新陳代謝；J：翻譯、核糖體結構和生物發生；K：轉錄；L：復制、重組和修復；M；細胞壁/膜/包膜生物發生；N：細胞運動；O：翻譯后修飾、蛋白質周轉和伴侶蛋白；P：無機離子轉運和代謝；Q：次生代謝產物的生物合成、轉運和分解代謝；R：僅通用功能預測；S：功能未知；T：信號轉導機制；U：細胞內運輸、分泌和囊泡運輸；V：防御機制；W：細胞外結構；Y：核結構；Z：碳骨架A: RNA processing and modification; B: Chromatin structure and dynamics; C: Energy production and conversion; D: Cell cycle control, cell division, chromosome partitioning; E: Amino acid transport and metabolism; F: Nucleotide transport and metabolism; G: Carbohydrate transport and metabolism; H: Coenzyme transport and metabolism; I: Lipid transport and metabolism; J: Translation, ribosomal structure and biogenesis; K: Transcription; L: Replication, recombination and repair; M: Cell wall/membrane/envelope biogenesis; N: Cell motility; O: Posttranslational modification, protein turnover, chaperones; P: Inorganic ion transport and metabolism; Q: Secondary metabolites biosynthesis, transport and catabolism; R: General function prediction only; S: Function unknown; T: Signal transduction mechanisms; U: Intracellular trafficking, secretion, and vesicular transport; V: Defense mechanisms; W: Extracellular structures; Y: Nuclear structure; Z: Cytoskeleton圖2 轉錄本的KOG注釋 Fig.2 KOG annotations of transcripts

表5 預測的轉錄因子家族和數量

2.5 蔗糖代謝相關基因克隆

根據紅肉火龍果成熟莖的全長轉錄組測序注釋結果，篩選蔗糖代謝相關的基因(表6)。獲得一條與菠菜蔗糖磷酸合酶(sucrose-phosphate synthase，SPS)具有81.03%相似性的轉錄本Transcript 920，其長度為3637 bp，含有3183 bp的CDS，編碼2060個氨基酸。獲得一條與紅莧菜蔗糖合酶(sucrose synthase，SUS)具有93.90%相似性的轉錄本Transcript 3494，其長度為2719 bp，含有2 412 bp的CDS，編碼803個氨基酸。獲得一條與藜麥液泡酸性轉化酶(acid beta-fructofuranosidase，AINV)具有72.61%相似性的轉錄本Transcript 7574，其長度為2229 bp，含有1935 bp的CDS，編碼644個氨基酸。

表6 全長轉錄組測序獲得的蔗糖代謝相關基因及序列相似性

此前針對紅肉火龍果成熟莖和成熟果實果肉，利用基于Illumina平臺的轉錄組測序進行分析，獲得若干在莖中高表達的蔗糖代謝相關Unigene[11]。SPS相關的c49878_g1長度為464 bp，與Transcript 920序列相似性為100.00%，其在成熟莖中的表達量約為成熟果實的3倍(表7)。SUS相關的c21997_g1長度為614 bp，與Transcript 3494序列相似性為99.16%，其在成熟莖中的表達量約為成熟果實的45倍。AINV相關的c34087_g1和c24912_g1，長度分別為665和320 bp。與Transcript 7574序列相似性分別為99.06%和99.08%，在成熟莖中的表達量分別約為成熟果實的27和72倍。

表7 蔗糖代謝基因與此前Illumina測序結果比較

3 討 論

對缺乏基因組序列的非模式作物開展功能基因組學研究，其前提是獲得基因全長序列。轉錄組測序能夠快速獲得大量基因的轉錄本序列，已在非模式作物中廣泛應用。此前普遍利用基于Illumina平臺的第二代高通量測序，其讀長較短，極難獲得基因全長序列。基于PacBio Sequel平臺的第三代高通量測序具有顯著的長讀長優勢，可用于全長轉錄組測序[14-15]。第三代高通量測序的主要缺點是錯誤率較高，表現為插入、缺失或錯配，均為隨機錯誤，可進行校正[16]。通過Subread聚類的自我校正，獲得的CCS可作為轉錄本的參考序列[17]。

此前利用基于Illumina平臺的轉錄組測序對紅肉火龍果成熟莖和成熟果實果肉進行分析，獲得Unigene平均長度為690 bp，長度大于1000 bp的Unigene僅占18.66%[11]。本研究獲得30 313條高質量轉錄本，平均長度為1829 bp，顯著長于Illumina測序。同時，長度大于1000 bp的轉錄本占比為88.96%，大于2000 bp占比為41.50%。由此可見，全長轉錄組測序獲得的轉錄本長度明顯長于Illumina測序，有利于獲得基因的全長序列。

獲得紅肉火龍果成熟莖的30 313條高質量轉錄本，僅有45.98%的序列在常見公共數據庫中被注釋，注釋到的物種主要為藜科(Chenopodiaceae)的甜菜和菠菜。與近期一些非模式作物的全長轉錄組測序相比，其轉錄本的注釋率明顯偏低。例如，對禾本科(Gramineae)薏苡屬(Coix)薏苡(C.lacryma-jobiLinn.)苗期葉片開展全長轉錄組測序，轉錄本平均長度為2318 bp，約91.50%轉錄本被注釋，物種為高粱、玉米、谷子、水稻和甘蔗等禾本科作物[18]。對薔薇科(Rosaceae)栒子屬(Cotoneaster)作物山東栒子(C.schantungensis)葉片、花和成熟果實進行全長轉錄組測序，約98.86%轉錄本被注釋，物種為蘋果、白梨、桃、梅、野草莓和沙梨等薔薇科作物[19]。對桑科(Moraceae)木波羅屬(Artocarpus)作物波羅蜜(A.heterophyllus)莖和葉進行全長轉錄組測序，獲得轉錄本平均長度為1684 bp，約97.56%轉錄本被注釋，物種為同屬于桑科的川桑(Morusnotabilis)[20]。上述作物在同一科內均有其他作物的基因組序列發布，在其親緣關系較近的情況下注釋成功率極高。紅肉火龍果屬石竹目(Caryophyllales)仙人掌科，該科尚未有物種的基因組序列發布。目前僅石竹目藜科的甜菜和菠菜基因組序列發布，因而相關轉錄本的注釋率較低。在未來仙人掌科作物基因組序列發布后，有必要對紅肉火龍果成熟莖的轉錄本序列重新注釋。

火龍果等仙人掌科作物的葉片退化為刺，其肉質莖是主要的光合器官。研究成熟莖的轉錄組模式，分離光合作用代謝途徑的相關基因全長，有助于研究蔗糖的合成、貯藏和轉運，進而調控果實品質形成。此前利用Illumina平臺的轉錄組測序對紅肉火龍果成熟莖和成熟果實果肉進行分析，僅獲得若干參與成熟莖中蔗糖代謝和轉運相關基因部分序列[11]。本研究獲得蔗糖磷酸合酶SPS、蔗糖合酶SUS、液泡酸性轉化酶AINV等均參與蔗糖合成和降解等代謝[21]。與Illumina轉錄組測序結果綜合分析，其較短的Unigene序列與其相應的全長轉錄本具有99%以上的核苷酸序列相似性，說明來源于同一基因。這些較短Unigene的數字表達模式表明其主要在成熟莖中表達，推測其參與源器官中蔗糖代謝。下一步可根據全長轉錄組測序獲得的基因全長克隆該基因，開展后續的功能驗證等分析。

4 結 論

研究利用基于Pacbio Sequel平臺的轉錄組測序對紅肉火龍果成熟莖進行全長轉錄組測序分析，獲得30 313個高質量轉錄本，其長度顯著長于Illumina測序結果。全長轉錄組測序獲得的高質量轉錄本，能結合Illumina測序獲得的轉錄本序列和基因數字表達結果，快速篩選若干候選基因的全長序列，為開展基因功能研究提供有力基礎。