蒺藜苜蓿（Medicago truncatula）全長轉錄組測序及 分析

尚驍堯 周玲芳 尹芊芊 晁躍輝

（北京林業大學草業與草原學院，北京 100083）

蒺藜苜蓿是一種非常重要的豆科模式植物，其具有基因組小，自花授粉，生命周期短，結實率高，具有穩定的遺傳轉化體系，越來越受到科研工作者的青睞。蒺藜苜蓿有2個品種已經被測序成功，其中，A17最早于2011年測序完成，目前已經更新至4.0版本［1-2］；另一個品種為R108，與A17相比，具有更高的遺傳轉化率，更短的生命周期，種子更易萌發，而且目前世界上最大的蒺藜苜蓿突變體庫就是使用的R108，所以被廣泛應用于基因功能研究，基于以上優點，R108蒺藜苜蓿于2014年被測序成功，并于2017年12月進行了版本更新［3-4］。

二代高通量測序（next-generation high-throughput sequencing，NGS）技術的出現，使高通量測序的成本有了大幅度的降低，使得轉錄組分析變得更加便捷。通過該技術，在蒺藜苜蓿中已經開展了大量的工作。通過NGS技術，分析一個MtTdp1α缺失的蒺藜苜蓿突變體，結果顯示該突變體中與DNA修復，DNA損傷以及染色體重組相關的基因轉錄水平發生了顯著的變化［5］。通過NGS技術分析蒺藜苜蓿根瘤的發生過程，鑒定出1 140個差異表達基因，其中很多涉及磷轉移的基因表達水平明顯上調；大量涉及根瘤形成的基因表達水平顯著上調，而抑制根瘤形成的基因表達水平下調［6］。利用該技術，在蒺藜苜蓿長期添加細胞分裂素及細胞分裂素抑制劑的過程中，鑒定出大量與蒺藜苜蓿響應細胞分裂素相關的基因，這些基因為研究蒺藜苜蓿生長發育和激素調節提供了候選［7］。

研究表明，植物在轉錄水平具有更加復雜的過程，如轉錄前水平調控，主要包括一些可變性剪接（alternative splice，AS）和可變多聚腺苷酸化（alternative polyadenylation，APA）［8-9］。比 如 在擬南芥中，80%的含有內含子的基因存在可變性剪接［10］，這些可變性剪接受到多種脅迫條件或者激素水平的調節［11-14］，同樣也發揮著重要的生物學功能。而這些內容難以通過NGS技術實現，因為NGS技術不能提供全長序列［15-16］。SMRT的出現，彌補了NGS技術的不足，因為通過SMRT可以獲得基因全長mRNA序列，準確地區分不同的剪接亞型（splice isoform）及鑒定APA位點等。SMRT在很多物種中有了成功的應用，如通過該技術，進行了人轉錄組的深度分析，發現了大概14 000個可變剪接基因，其中有＞10%是從來沒有注釋的［17］。應用該技術，在小麥中發現了3 026個新基因，在玉米中鑒定大量的新的基因亞型，大大豐富了這些植物基因組信息［18-19］。在一些豆科植物中也應用該技術進行全長轉錄組分析，如在紅三葉中鑒定出5 492條可變性剪接，4 333條lncRNA以及3 762條融合轉錄本，這些信息在以往的基因組測序和二代轉錄組分析中從未出現［20］。雖然蒺藜苜蓿的基因組測序工作已經完成，但是AS、lncRNA、APA位點以及融合轉錄本信息的不足，如目前在蒺藜苜蓿R108中無融合轉錄本數據，以上情況嚴重制約了蒺藜苜?；蚪M深入發掘。

本研究利用SMRT技術對蒺藜苜蓿進行全長轉錄組測序及分析，通過與蒺藜苜蓿參考基因組進行數據比對，鑒定新基因及轉錄本，分析蒺藜苜蓿中的AS類型并統計不同類型的數量，同時預測蒺藜苜蓿中lncRNA及融合轉錄本。這些轉錄組數據能夠有效提高蒺藜苜?；蚪M注釋水平，極大豐富基因序列及結構信息，旨為深入發掘蒺藜苜蓿資源提供有用資源。

1 材料與方法

1.1 材料

試驗中使用的蒺藜苜蓿R108種子為Samuel Roberts Noble Foundation惠贈，使用草炭灰∶蛭石∶珍珠巖（1∶1∶1）作為營養土，植物材料放置于人工氣候箱中，生長條件為：白天25℃ 16 h，晚上22℃ 8 h；相對濕度為50%-70%。為了獲得更多的轉錄組測序信息，選取不同的植物組織（包括根、莖、葉、花、莢果）；對蒺藜苜蓿進行不同的脅迫處理（100 mmol/L NaCl和10% PEG2000，處理1 d）；對蒺藜苜蓿進行不同的激素誘導（10 μmol/L ABA、10 μmol/L GA3、10 μmol/L IAA和10 μmol/L 6-BA，處理6 h）；分別收集不同的植物材料，混合為一個樣品，進行全長轉錄組測序。

1.2 方法

1.2.1 建庫及SMRT測序 通過植物RNA提取試劑盒分離植物樣品總RNA，去除基因組DNA后，對提取的RNA進行純度及完整性檢測。吸取5 μg總RNA用于建庫反應，基于大小分為2個文庫（＜4 kb與＞4 kb），使用Pacific Bioscience Sequel平臺對這2個文庫進行SMRT測序。

1.2.2 數據處理 測序獲得的原始數據（raw read）通過SMRTlink軟件（version 4.0）進行處理，運行參 數 為：minLength=200，minReadScore=0.75。處理后的數據使用subread BAM files文件處理后，生成環狀一致性序列（circular consensus sequence，CCS），參 數 設 置 為：parameters：min_length 200，max_drop_fraction 0.8，no_polish TRUE，min_zscore -999，min_passes 1，min_predicted_accuracy 0.8，max_length 18000。通過分析5'接頭、3'接頭及ploy（A）結構，鑒定出FLNC序列，經聚類、校正后，最終獲得校正序列。

1.2.3 基因組比對及序列結構分析 SMRT測序結果，通過GMAP軟件（version：2017-01-14）與蒺藜苜蓿R108參考基因組（R108_HM340_v1.0）進行比對。通過TAPIS pipeline軟件（Version 1.2.1）進行基因結構及APA分析［21］，通過SUPPAA軟件（Version：2017-02-07）進行AS分析［22］，通過MEME對轉錄本的poly（A）位點上游50堿基序列進行poly（A）信號分析，同時對轉錄本進行基因組跨區域分析，以鑒定融合轉錄本［19］。

1.2.4 開放閱讀框、轉錄因子及lncRNA鑒定 校正后的序列的開放閱讀框（open reading fragment，ORF）使用ANGEL pipeline進行分析，其中包含有完整ORF、5'非編碼區及3'非編碼區的轉錄本被鑒定為全長轉錄本，并進行編碼氨基酸序列分 析［23］。對 轉 錄 組 數 據 使 用CPC［24］、CNCI［25］、PLEK［26］以及Pfam［27］4種方法進行lncRNA序列預測，并對預測的lncRNA進行分類［28］。下載已知的蒺藜苜蓿lncRNA（https://greenc.sciencedesigners.com/api/db/greenc/species/Medicago_truncatula/fasta），使用BLASTN程序進行序列比對以尋找新lncRNA，參數設置：e-value ≤1e-10，min-identity=90%，mincoverage=85%。

1.2.5 RT-PCR驗證 為了驗證SMRT分析的新基因及融合轉錄本預測結果，通過RT-PCR方法進行驗證。提取蒺藜苜蓿總RNA，去除基因組DNA后，反轉錄為cDNA。基于SMRT分析結果，通過DNA Man（Version 6.0） 和Primer Premier（Version 5.0）軟件進行引物設計。以cDNA為模板，進行PCR反應，待反應結束，取5 μL PCR產物經1%瓊脂糖凝膠電泳分析，選取含有目標大小的PCR產物，送北京睿博生物技術公司測序。

2 結果

2.1 SMRT測序數據統計

通過對2個文庫進行SMRT測序，一共獲得了7 728 183個subread（共計11.77 Gb數據量），通過分析得知平均長度為1 524 bp，N50為2 051 bp （表1）。為了獲得更精確的數據信息，這些數據生成了659 220條CCS，共鑒定出513 951條全長序列，其中509 014條為FLNC序列，平均長度為2 047 bp（表1）。這些FLNC序列被聚類為243 832一致性序列（ICE consensus，表1），聯合非全長序列，共獲得了243 676高質量的全長優化序列（polished consensus），平均長度為2 321 bp，N50為3 380 bp（表1）；這些優化后序列經NGS序列校正后，共生成了236 243條校正序列（corrected consensus），這些序列平均長度為2 425 bp，N50為3 496 bp；對校正后的序列進行分析得知：其中220 233條（占93.22%）長度＞500 bp，180 874條（占76.56%）長度＞1 kb（表1）。

表1 SMRT測序數據統計Table 1 Statistics of sequencing data by SMRT

2.2 新基因及轉錄本鑒定

校正后的序列與蒺藜苜?；蚪M比對后，發現這些校正后的序列可以分為4類（圖1）：（1）未比對到參考基因組上的（unmapped）13 327條（占5.64%）；（2）與參考基因組多個區域比對上的（multiple mapped）18 393條（占7.79%）；（3）與參考基因組正向序列比對上的（mapped to ‘+’）148 971條（占63.06%）；（4）與參考基因組反向序列比對上的（mapped to ‘-’）55 552條（占23.51%）。

圖1 校正后序列與參考基因組比對分析Fig. 1 GMAP analysis of corrected reads to reference genome

與參考基因組比對上的序列經過進一步的數據過濾（去除可能存在的錯誤序列及重復序列）后，共得到了61 125條最終序列（RC61125），這些序列可以分為3種類型（圖2）：（1）8 489條與參考基因組完全一致；（2）43 592條為已知基因的新轉錄本；（3）9 044條來源于新的基因?？傊?，通過SMRT測序共鑒定出7 209新基因以及52 636條新轉錄本。

圖2 RC61125序列的分類Fig. 2 Classification of RC61125 sequence

對SMRT測序鑒定的7 209新基因使用進行功能注釋，結果顯示，共有96.85%的新基因被注釋成功（表2）。

表2 新基因功能注釋Table 2 Functional annotations of novel genes

將這些新基因與NR數據庫比對，發現數量分布最多的前5物種為：蒺藜苜蓿A17（M. truncatula，4 128）、地三葉（Trifolium subterraneum，594）、木豆（Cajanus cajan，100）、鷹嘴豆（Cicer arietinum，70）和大豆（Glycine max，52），這些物種主要為豆科植物（圖3）。KEGG分析顯示：4 532個新基因可以富集到298條KEGG通路中；GO分析結果顯示：1 717個新基因主要可以分為生物進程，分子功能和細胞組分。

圖3 新基因Nr同源物種分布圖Fig. 3 Nr homologous species distribution diagram of novel genes

為了驗證SMRT發現新基因的準確性，隨機選擇5個新基因進行RT-PCR驗證。結果（圖4）顯示，基于SMRT技術發現的新基因序列，設計的引物，均能擴增出目標大小DNA片段，PCR產物經生物技術公司測序也驗證了序列的正確性。以上結果表明，發掘蒺藜苜蓿新基因，SMRT是一種可靠的技術手段。

圖4 新基因RT-PCR驗證Fig. 4 RT-PCR validation of novel genes

2.3 基因結構分析

ORF分析顯示共鑒定出36 235條編碼蛋白質序列，這些序列平均長度為1 027.03 bp，其中13 451條含有完整編碼區（表3）。對序列的5'-UTR及 3'-UTR數量與長度分布統計結果顯示含有11 071個5'-UTR，平均長度為992.86 bp，含有1 990個3'-UTR，平均長度為799.00 bp（表3）。

表3 編碼區及非編碼區鑒定Table 3 CDS and UTR identification

RC61125用于外顯子結構分析，結果（圖5）顯示，17 590條（占28.78%）只含有1個外顯子，8 327條（占13.62%）含有2個外顯子，外顯子數目＞20的序列為1 586條（占2.59%）；基于已經公布的2種蒺藜苜蓿參考基因組序列分析，在R108中，15 865條（占26.00%）轉錄本只含有1個外顯子，而外顯子數目＞20的序列僅為612條（占1.00%），在A17中，13 594條（占23.17%）轉錄本只含有1個外顯子，而外顯子數目＞20的序列僅為696條（占1.19%，圖5）。內含子數目及結構分析，結果顯示，參考基因組R108中內含子數目為184 641個，平均每個基因3.31個，而SMRT共測得269 289個內含子，平均每個基因11.26個。

圖5 轉錄本外顯子分布Fig.5 Distribution of exons in transcripts

2.4 AS及splice isoform分析

通過AS分析顯示，SMRT測序共鑒定了11 761個AS，該數量在參考基因組中僅為1 468個；AS類型分析顯示：SMRT技術鑒定的主要AS類型為RI，該特征與其他植物一致，而參考基因組中主要AS為外顯子跳躍（skipped exon，SE）（圖6）；與參考基因組比較發現，SMRT測序鑒定出9 582個新AS。

圖6 AS數量及類型 Fig.6 Numbers and types of AS events

基因轉錄本分析顯示，在參考基因組中，只有1條轉錄本的基因數目為52 332（占93.94%），而含有＞10條轉錄本的基因數目為5，其中有2個基因（maker-Contig51-snap-gene-51.43 和maker-Contig137-exonerate_est2genome-gene-5.0）的轉錄本數目最多為13（表4）；SMRT數據分析顯示，含有1條轉錄本的基因數目為11 730（占49.03%），含有2條轉錄本的基因數目為4 713（占19.70%），而含有＞10條轉錄本的基因數目為514（占2.15%，表4）。如基于SMRT分析發現，在蒺藜苜蓿 maker-Contig176-snap-gene-26.47基因中含有6條轉錄本，而在參考基因組中僅注釋為1個。

表4 SMRT測序及參考基因組轉錄本數目Table 4 Number of transcripts in SMRT sequence and reference genomes

2.5 lncRNA及融合轉錄本鑒定

通過4種方法，共鑒定出6 595條lncRNA，平均長度為1 250 bp，其中5 265條（占79.83%）僅含有1個外顯子（圖7-A）。根據基因組上的位置關系，將lncRNA分為4種類型（圖7-B）：（1）基因間區的lncRNA（intergenic lncRNA，lincRNA）3 669條（占55.65%）；（2）反義lncRNA（antisense）779條（占11.8%）；（3）內含子型lncRNA（intronic）340條（占5.16%）；（4）與mRNA外顯子有重疊區lncRNA（sense overlapping lncRNA）1 807條（占27.40%）。通過與參考基因組的注釋lncRNA對比分析，發現共鑒定出6 503個新lncRNA。

圖7 lncRNA鑒定Fig.7 Identification of lncRNA

融合轉錄本預測顯示，通過SMRT分析共發現了479條融合轉錄本，而不同染色體間（443條）發生融合轉錄本的概率遠遠大于同一染色體內（36條）。為了驗證SMRT發現融合轉錄本的準確性，隨機選擇了5個融合轉錄本，進行RT-PCR驗證。結果（圖8）顯示，基于SMRT技術發現的融合轉錄本序列，設計的引物，均能擴增出目標大小的DNA片段，PCR產物經生物技術公司測序也驗證了序列的正確性。

圖8 融合轉錄本RT-PCR驗證Fig. 8 RT-PCR validation of fusion transcripts

2.6 APA分析

通過分析發現，SMRT測序共鑒定了23 926個基因，其中12 049個基因（含有295 545轉錄本）至少存在1個poly（A）位點，441個基因存在至少5個poly（A）位點，平均每個基因含有2.14個poly（A）位點（圖9）。對poly（A）剪切位點上、下游各50 bp核酸組成分析顯示，蒺藜苜蓿poly（A）剪切位點上游富含尿嘧啶（uracil，U），下游富含腺嘌呤（adenine，A）；通過MEME進行保守元件分析，結果顯示，在蒺藜苜蓿poly（A）剪切位點上游存在2個保守元件（AAUAAA與UGUA，圖9）。

圖9 poly（A）位點分析Fig. 9 Analysis of poly（A）sites

3 討論

蒺藜苜蓿作為一種非常重要的豆科模式植物，雖然基因組序列已經公布，但是其轉錄組特征和mRNA結構并沒有深入分析。隨著測序技術的進步，SMRT為該部分工作的進行提供了新的途徑。

全長轉錄組序列對于植物基因組注釋和功能研究都具有極其重要的作用，通過SMRT對蒺藜苜蓿全長轉錄組進行測序，共產生了659 220 CCS，其中509 014被鑒定為FLNC轉錄本。FLNC序列的長度可以用于評估全長轉錄組建庫的優劣，通過長度分析，發現FLNC的序列分布與2個測序文庫長度大小一致（圖1）。同時CCS與FLNC序列不但從長度上具有二代測序不可比擬的優點，同時可以避免二代測序中短序列拼接的帶來的序列錯誤。通過分析蒺藜苜蓿二代測序數據，共獲得了236 243條校正序列，其中76.56%長度＞1 kb，13 451條序列含有完整ORF，該數字遠遠低于三代測序的質量。通過SMRT測序，共鑒定出23 926個基因，平均長度為4 307 bp，不但比二代測序長度要高，而且比參考基因組中還要高2 076 bp。通過SMRT對其他植物測序研究，同樣也支持了類似的觀點，如：使用SMRT對小麥轉錄組測序分析工作，結果顯示，測序平均長度比小麥參考基因組注釋序列要高45 bp［21］。對紫花苜蓿同時進行SMRT及二代測序和分析，其中75.68%的二代測序長度低于1 000 bp，而SMRT測序中僅有3.76%序列長度低于1 000 bp［29］。

SMRT測序在新基因與新轉錄本發現上也具有很大的優勢。通過SMRT測序，在毛竹中發現了8 091個新轉錄本［30］，在穿山甲中鑒定出8 014個新基因［31］。本研究中，鑒定出了新基因7 209個，已知基因的新轉錄本43 592個及新基因新轉錄本9 044個。RT-PCR及PCR產物測序，也進一步證實SMRT測序在新基因及新轉錄本鑒定上的準確性。SMRT測序還能更好的分析AS，使用該技術，在蒺藜苜蓿中共鑒定了11 761 AS，而在參考基因組中該數目僅為1 468。其他植物中公布的數據來看，主要的AS類型為RI，如玉米［19］、毛竹［30］、草莓［32］、無油樟［33］等，而蒺藜苜蓿R108參考基因組的數據顯示，主要的AS類型為SE。通過SMRT測序得知，蒺藜苜蓿R108中AS中發現9 582個新AS，而RI出現頻率最高。在蒺藜苜蓿R108中共鑒定了479條融合轉錄本，融合轉錄本發生在不同染色體間的頻率明顯高于染色體內。蒺藜苜蓿轉錄本的融合特征與在玉米中發現的規律一致［19］。lncRNA作為目前研究的熱點，在植物生長發育中起著重要的調控作用［34-36］，本研究共鑒定了6 595條lncRNA，主要類型為lincRNA。蒺藜苜蓿參考序列中已經注釋的lncRNA為9 676個，經過比對分析，通過SMRT鑒定的lncRNA中92個為已知lncRNA，而6 503個為新lncRNA。SMRT鑒定的lncRNA平均長度為1 250 bp，最長為8 175 bp，而目前已注釋的lncRNA平均長度為372 bp，最大長度為3 890 bp。該結果與SMRT技術在其他植物中的發現一致，在玉米中lncRNA能達到6 kb［19］。通過SMRT技術共鑒定了23 926個基因，其中12 049個基因能夠產生295 545條轉錄本，這些基因至少含有一個poly（A）位點。這些數據豐富了蒺藜苜蓿轉錄組信息，同時也為蒺藜苜?；蚪M版本的升級提供了支持。

4 結論

通過SMRT技術對蒺藜苜蓿進行全長轉錄組測序及分析，共鑒定出7 209新基因，52 636條新轉錄本，9 582個新AS，6 503個新lncRNA及479條融合轉錄本；對基因結構和APA特征進行了分析；對蒺藜苜蓿R108參考基因組進行了數據補充，同時將原有的主要的AS類型修正為RI。