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

（北京林業(yè)大學草業(yè)與草原學院，北京 100083）

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

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

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

本研究利用SMRT技術(shù)對蒺藜苜蓿進行全長轉(zhuǎn)錄組測序及分析，通過與蒺藜苜蓿參考基因組進行數(shù)據(jù)比對，鑒定新基因及轉(zhuǎn)錄本，分析蒺藜苜蓿中的AS類型并統(tǒng)計不同類型的數(shù)量，同時預測蒺藜苜蓿中l(wèi)ncRNA及融合轉(zhuǎn)錄本。這些轉(zhuǎn)錄組數(shù)據(jù)能夠有效提高蒺藜苜蓿基因組注釋水平，極大豐富基因序列及結(jié)構(gòu)信息，旨為深入發(fā)掘蒺藜苜蓿資源提供有用資源。

1 材料與方法

1.1 材料

試驗中使用的蒺藜苜蓿R108種子為Samuel Roberts Noble Foundation惠贈，使用草炭灰∶蛭石∶珍珠巖（1∶1∶1）作為營養(yǎng)土，植物材料放置于人工氣候箱中，生長條件為：白天25℃ 16 h，晚上22℃ 8 h；相對濕度為50%-70%。為了獲得更多的轉(zhuǎn)錄組測序信息，選取不同的植物組織（包括根、莖、葉、花、莢果）；對蒺藜苜蓿進行不同的脅迫處理（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）；分別收集不同的植物材料，混合為一個樣品，進行全長轉(zhuǎn)錄組測序。

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 數(shù)據(jù)處理 測序獲得的原始數(shù)據(jù)（raw read）通過SMRTlink軟件（version 4.0）進行處理，運行參 數(shù) 為：minLength=200，minReadScore=0.75。處理后的數(shù)據(jù)使用subread BAM files文件處理后，生成環(huán)狀一致性序列（circular consensus sequence，CCS），參 數(shù) 設(shè) 置 為：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）結(jié)構(gòu)，鑒定出FLNC序列，經(jīng)聚類、校正后，最終獲得校正序列。

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

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

1.2.5 RT-PCR驗證 為了驗證SMRT分析的新基因及融合轉(zhuǎn)錄本預測結(jié)果，通過RT-PCR方法進行驗證。提取蒺藜苜?？俁NA，去除基因組DNA后，反轉(zhuǎn)錄為cDNA?；赟MRT分析結(jié)果，通過DNA Man（Version 6.0） 和Primer Premier（Version 5.0）軟件進行引物設(shè)計。以cDNA為模板，進行PCR反應，待反應結(jié)束，取5 μL PCR產(chǎn)物經(jīng)1%瓊脂糖凝膠電泳分析，選取含有目標大小的PCR產(chǎn)物，送北京睿博生物技術(shù)公司測序。

2 結(jié)果

2.1 SMRT測序數(shù)據(jù)統(tǒng)計

通過對2個文庫進行SMRT測序，一共獲得了7 728 183個subread（共計11.77 Gb數(shù)據(jù)量），通過分析得知平均長度為1 524 bp，N50為2 051 bp （表1）。為了獲得更精確的數(shù)據(jù)信息，這些數(shù)據(jù)生成了659 220條CCS，共鑒定出513 951條全長序列，其中509 014條為FLNC序列，平均長度為2 047 bp（表1）。這些FLNC序列被聚類為243 832一致性序列（ICE consensus，表1），聯(lián)合非全長序列，共獲得了243 676高質(zhì)量的全長優(yōu)化序列（polished consensus），平均長度為2 321 bp，N50為3 380 bp（表1）；這些優(yōu)化后序列經(jīng)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測序數(shù)據(jù)統(tǒng)計Table 1 Statistics of sequencing data by SMRT

2.2 新基因及轉(zhuǎn)錄本鑒定

校正后的序列與蒺藜苜?；蚪M比對后，發(fā)現(xiàn)這些校正后的序列可以分為4類（圖1）：（1）未比對到參考基因組上的（unmapped）13 327條（占5.64%）；（2）與參考基因組多個區(qū)域比對上的（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

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

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

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

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

將這些新基因與NR數(shù)據(jù)庫比對，發(fā)現(xiàn)數(shù)量分布最多的前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分析結(jié)果顯示：1 717個新基因主要可以分為生物進程，分子功能和細胞組分。

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

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

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

2.3 基因結(jié)構(gòu)分析

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

表3 編碼區(qū)及非編碼區(qū)鑒定Table 3 CDS and UTR identification

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

圖5 轉(zhuǎn)錄本外顯子分布Fig.5 Distribution of exons in transcripts

2.4 AS及splice isoform分析

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

圖6 AS數(shù)量及類型 Fig.6 Numbers and types of AS events

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

表4 SMRT測序及參考基因組轉(zhuǎn)錄本數(shù)目Table 4 Number of transcripts in SMRT sequence and reference genomes

2.5 lncRNA及融合轉(zhuǎn)錄本鑒定

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

圖7 lncRNA鑒定Fig.7 Identification of lncRNA

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

圖8 融合轉(zhuǎn)錄本RT-PCR驗證Fig. 8 RT-PCR validation of fusion transcripts

2.6 APA分析

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

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

3 討論

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

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

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

4 結(jié)論

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