基于Iso-Seq技術(shù)的野生馬齒莧葉片轉(zhuǎn)錄組分析

葉梅榮，黃守程，王曉鵬，劉愛榮，崔 峰，康 健

(安徽科技學(xué)院 生命與健康科學(xué)學(xué)院， 安徽 鳳陽 233100)

全長(zhǎng)轉(zhuǎn)錄組測(cè)序(isoform sequencing，Iso-Seq)是指利用三代測(cè)序技術(shù)獲得生物體內(nèi)完整的RNA分子的全長(zhǎng)序列，建庫過程中RNA分子無需打斷，測(cè)序可直接獲得完整的cDNA序列[1]。第三代測(cè)序技術(shù)可以獲得更長(zhǎng)乃至全長(zhǎng)的轉(zhuǎn)錄組[2]。三代單分子測(cè)序技術(shù)的應(yīng)用降低了無參考轉(zhuǎn)錄組的分析難度，能夠較為容易地獲取完整基因組和全長(zhǎng)轉(zhuǎn)錄本，對(duì)深入研究生物體的轉(zhuǎn)錄機(jī)制極有利。三代測(cè)序是無需PCR擴(kuò)增、操作更加簡(jiǎn)單、通量高、速度快、讀長(zhǎng)堿基對(duì)的新型測(cè)序技術(shù)，但比二代測(cè)序錯(cuò)誤率高，因此結(jié)合三代和二代測(cè)序技術(shù)更有利于研究生物體的轉(zhuǎn)錄組[3]。三代測(cè)序技術(shù)剛興起，但是發(fā)展很快，如Kuang等[4]利用Iso-seq探討東北紅豆杉(TaxuscuspidateSieb. et Zucc.)的轉(zhuǎn)錄組，揭示了東北紅豆杉生物合成的復(fù)雜性。

馬齒莧(PortulacaoleraceaL.)別稱有馬莧、五行草、長(zhǎng)命菜、五方草、馬齒菜等，為馬齒莧科 (Portulacaceae)馬齒莧屬(Portulaca)[5]，分布地很廣，全球各地均有分布，被列為世界第八大普通的植物[6]。馬齒莧既可入藥，又可食用，是藥食兩用的野生植物。近年的藥理研究表明，馬齒莧具有抗菌、降血脂、抗衰老、松弛肌肉、抗炎、鎮(zhèn)痛，以及促進(jìn)傷口愈合等功效[7-9]。馬齒莧營(yíng)養(yǎng)豐富，富含多種礦質(zhì)元素和多種有機(jī)營(yíng)養(yǎng)保健成分如多糖、總黃酮、α-亞麻酸、β-胡蘿卜素、ω-3脂肪酸、萜類、生物堿、去甲腎上腺素、褪黑激素等，具有改善血液循環(huán)、提高人體免疫力、防治心血管疾病、抑制微生物生長(zhǎng)等多重功效[6-8，10]。馬齒莧適應(yīng)性極強(qiáng)，對(duì)高溫、干旱、高濕、高鹽、重金屬污染等逆境的抵抗能力強(qiáng)大，是優(yōu)良的生態(tài)修復(fù)植物[11-14]。馬齒莧集保健、食用、藥用、生態(tài)價(jià)值于一身，具有進(jìn)一步研究開發(fā)的價(jià)值。

馬齒莧耐高溫和干旱，不耐低溫，因此本研究利用干旱和低溫脅迫分別處理馬齒莧，同時(shí)設(shè)置對(duì)照，目前馬齒莧沒有參考基因組，因此本研究基于三代(Iso-Seq)測(cè)序技術(shù)和二代測(cè)序技術(shù)研究馬齒莧葉片的全長(zhǎng)轉(zhuǎn)錄組，以期為馬齒莧屬的基因研究提供參考。

1 材料與方法

1.1 材料

馬齒莧種子播種于花盆中(盆高14 cm，下口直徑11 cm，上口直徑23 cm)，基質(zhì)是通用型“中禾”營(yíng)養(yǎng)土，光照培養(yǎng)箱培養(yǎng)溫度白天30 ℃，晚上25 ℃。待種子萌發(fā)成幼苗后，每盆留有5株幼苗。待苗長(zhǎng)至5 cm 高度左右，設(shè)置低溫處理(LTS，白天10 ℃，晚上5 ℃)和對(duì)照處理即常溫處理(C，白天30 ℃，晚上25 ℃)，同時(shí)設(shè)置自然干旱處理(DS，白天30 ℃，晚上25 ℃，即一次澆透水后，連續(xù)14 d不澆水)，對(duì)照組(即常溫處理組)每7 d澆透1次水，每個(gè)處理重復(fù)3次。低溫處理和干旱處理14 d后，分別取各處理的馬齒莧成熟葉片用液氮速凍，然后干冰運(yùn)送到武漢貝納科技服務(wù)有限公司進(jìn)行三代和二代轉(zhuǎn)錄組測(cè)序和分析。馬齒莧在培養(yǎng)過程中因處理時(shí)間較長(zhǎng)已進(jìn)入開花結(jié)實(shí)階段。

1.2 總RNA提取、cDNA文庫構(gòu)建和二代技術(shù)轉(zhuǎn)錄組測(cè)序

提取馬齒莧樣本的RNA，用瓊脂糖凝膠電泳檢測(cè)RNA質(zhì)量，用NanoDrop ND-1000超微量分光光度計(jì)檢測(cè)RNA的純度和濃度，采用Agilent2100生物分析儀精確檢測(cè)RNA的完整性，利用Illumina HiSeq 2500測(cè)序平臺(tái)進(jìn)行高通量測(cè)序。

1.3 文庫構(gòu)建和三代技術(shù)轉(zhuǎn)錄組測(cè)序

利用Oligo(dT)富集帶有PolyA尾的RNA，使用SMARTerPCR cDNA Synthesis Kit將RNA反轉(zhuǎn)錄合成cDNA，使用KAPA HiFi PCR Kits進(jìn)行PCR擴(kuò)增合成cDNA；利用BluePippin進(jìn)行片段篩選，用SMRTbell template prep kit 1.0構(gòu)建SMRTbell文庫，進(jìn)行損傷修復(fù)，末端修復(fù)；在DNA片段兩端連接莖環(huán)狀測(cè)序接頭，并利用外切酶去除連接失敗的片段。文庫定量后，利用PacBio Sequel Ⅱ?qū)DNA文庫測(cè)序，獲得原始下機(jī)數(shù)據(jù)(polymerase reads)。

1.4 數(shù)據(jù)質(zhì)控

對(duì)polymerase reads進(jìn)行過濾，去除adaptor序列、長(zhǎng)度小于50 bp和準(zhǔn)確度小于80%的polymerase reads，然后計(jì)算合格的插入片段(subreads)的數(shù)目、長(zhǎng)度分布和N50，利用CCS軟件校正獲得高質(zhì)量的插入片段序列(reads of insert，ROI)即一致性序列(Circular consensus sequencing，CCS)，利用smrtlink 5軟件對(duì)插入片段的接頭序列進(jìn)行識(shí)別和去除，得到初級(jí)轉(zhuǎn)錄本，初級(jí)轉(zhuǎn)錄本存在大量的冗余序列。利用ICE(isoform-level clustering algorithm)算法將冗余序列聚類到一起，得到新的一致性序列，最終得到一致性序列。再利用CD-HIT(version4.6.7)軟件去除冗余序列分別得到isoform級(jí)別序列，將isoform級(jí)別序列繼續(xù)去除冗余序列得到unigene。所有數(shù)據(jù)已上傳美國(guó)國(guó)家生物技術(shù)信息中心(BioProject： PRJNA734447)。

1.5 Unigene功能注釋

將universal gene(unigene)與universal protein(Uniprot)、Pfam protein families(Pfam)、gene ontology(GO)、kyoto encyclopedia of genes and genomes(KEGG)、evolutionary genealogy of genes： non-supervised orthologous groups(eggNOG)、KEGG pathway(pathway)和non-redundant protein(Nr)數(shù)據(jù)庫進(jìn)行同源比對(duì)，獲得所有注釋信息。并對(duì)注釋到各數(shù)據(jù)庫的信息功能注釋和分類。

1.6 簡(jiǎn)單重復(fù)序列標(biāo)記(simple sequence repeats，SSR)信息分析

對(duì)unigene進(jìn)行SSR(Simple Sequence Repeats)位點(diǎn)查找，設(shè)置參數(shù)為：(1/10)(2/6)(3/5)(4/5)(5/4)(6/4)，其中(3/5)代表3個(gè)堿基重復(fù)5次或以上，以此類推。

1.7 差異表達(dá)基因篩選

二代原始測(cè)序數(shù)據(jù)經(jīng)過過濾得到clean reads，因馬齒莧沒有參考基因組，因此與三代測(cè)序得到的全長(zhǎng)非嵌合的轉(zhuǎn)錄本經(jīng)CD-HIT去冗余后得到的isoform作為參考序列進(jìn)行比對(duì)，用RSEM對(duì)比對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)，以FPKM(fragments per kilobase of transcript per million fragments mapped)≥1為標(biāo)準(zhǔn)篩選有效表達(dá)基因，利用DEGseq軟件，以Q value<0.05，|log2FoldChange|>2為標(biāo)準(zhǔn)篩選差異表達(dá)基因(differentially expressed genes，DEGs)。

2 結(jié)果與分析

2.1 初級(jí)轉(zhuǎn)錄本

為了得到盡可能多的轉(zhuǎn)錄本，將低溫處理、干旱處理及對(duì)照組(每個(gè)處理重復(fù)3次)共9個(gè)樣本的葉片混合一起，反轉(zhuǎn)錄成cDNA，總共得到641 732條原始片段，N50為144 732，質(zhì)控后得到插入片段26 758 071，總堿基量為44 G，其N50為2 002，利用CCS軟件質(zhì)控得到高質(zhì)量的插入片段442 760條(表1)，N50為2 397，說明測(cè)序質(zhì)量較好。利用smrtlink軟件對(duì)高質(zhì)量的插入片段進(jìn)行識(shí)別和去除接頭序列，得到全長(zhǎng)的插入片段和全長(zhǎng)的去除嵌合體的插入片段分別為391 258條和301 412條，全長(zhǎng)片段所占比例越高，說明測(cè)序質(zhì)量越好。通過ICE算法去除冗余序列得到188 792條一致序列。利用CD-HIT軟件將相同轉(zhuǎn)錄本聚類到一起，得到103 298條isoform序列。

表1 初級(jí)轉(zhuǎn)錄本統(tǒng)計(jì)

2.2 Unigene級(jí)別聚類分布

根據(jù)CD-HIT軟件進(jìn)一步將來源于一個(gè)基因的isoform聚類到一起，即從103 298條去除冗余得到39 717條unigene序列(圖1)。unigene集中分布在600～3 200 bp，所占百分比為88.34%。

圖1 Unigenes長(zhǎng)度分布圖

其中1 200～1 400 bp的unigene數(shù)量最多，具有3 962條，占百分比為9.98%；其次是1 400～1 600 bp的基因，有3 911條，占百分比為9.85%。

2.3 LncRNA預(yù)測(cè)

LncRNA(長(zhǎng)鏈非編碼RNA)是一類轉(zhuǎn)錄本長(zhǎng)度超過200核苷酸(nt)，不編碼蛋白質(zhì)的RNA分子。根據(jù)核心程序plek采用ncrna_pipeline流程對(duì)得到的去冗余轉(zhuǎn)錄本序列預(yù)測(cè)LncRNA，得到40 952條沒有編碼能力的轉(zhuǎn)錄本序列，長(zhǎng)度主要分布在400～2 200 nt，占比95.10%。其中長(zhǎng)度600～800 nt之間分布數(shù)量最多，為7 636條，所占百分比為18.65%；其次是長(zhǎng)度為800～1 000 nt的基因，有7 514條，占比為18.35%(圖2)。

圖2 LncRNA長(zhǎng)度分布圖

2.4 ORF預(yù)測(cè)

開放閱讀框(open reading frame，ORF)是蛋白質(zhì)編碼基因的正常核苷酸序列，從起始密碼子到終止密碼子的閱讀框可編碼完整的多肽鏈，期間不存在使翻譯中斷的終止密碼子。利用transDecoder(版本：5.5.0)對(duì)unigene進(jìn)行讀碼框預(yù)測(cè)，鑒定標(biāo)準(zhǔn)為ORF長(zhǎng)度大于300 bp，該序列似然函數(shù)的log值大于0，選取6個(gè)ORF中分?jǐn)?shù)最大的一個(gè)，選擇最長(zhǎng)的一個(gè)，一共預(yù)測(cè)出38 419條CDS序列，其長(zhǎng)度分布圖見圖3。從圖3可以看出，CDS集中分布在200～1 600 bp，占比為83.25%，其中分布最多的是200～400 bp的條帶，有7 464條，占比為19.43%；其次為400～600 bp的條帶，有6 342條，占比為16.51%。

圖3 CDS長(zhǎng)度分布

2.5 Unigene的功能注釋

將樣本unigene與Nr、Uniprot、Pfam、GO、KEGG、eggNOG和Pathway數(shù)據(jù)庫進(jìn)行同源比對(duì)，被注釋的基因有36 381條，注釋率為91.6%；未注釋的基因有3 336條。注釋率最高的為Nr，注釋了36 034條基因，注釋率為90.73%；其次是Uniprot，注釋了29 464條基因，注釋率為74.18%(圖4)。

圖4 Unigene 條目的注釋統(tǒng)計(jì)

2.5.1 GO分類

對(duì)注釋上GO分類的28 875個(gè)基因進(jìn)行功能的預(yù)測(cè)(圖5)，主要分為3大類：分子功能、細(xì)胞組分和生物學(xué)過程，這3大類最多的分別是離子結(jié)合、細(xì)胞組分和生物學(xué)合成過程，分別有11 385、8 272和8 055個(gè)基因。

圖5 GO功能注釋

2.5.2 KEGG pathway分類

對(duì)注釋上KEGG的16 861個(gè)基因進(jìn)行代謝通路分類(圖6)，總共分為5大類，包括有機(jī)系統(tǒng)、代謝、遺傳信息過程、環(huán)境信息過程和細(xì)胞過程；這5大類程基因數(shù)量最多的分別為內(nèi)分泌系統(tǒng)、碳水化合物代謝過程、翻譯、信號(hào)轉(zhuǎn)導(dǎo)及運(yùn)輸和分解代謝，其基因數(shù)量分別為1 057、3 105、1 521、2 442和1 109個(gè)。從圖中還可以看出代謝所占基因數(shù)最多，其占總注釋基因數(shù)55.25%。

圖6 KEGG功能注釋

2.5.3 COG功能分類

利用COG對(duì)基因進(jìn)行直系同源分類，分為26個(gè)組(圖7)，然后對(duì)unigene進(jìn)行功能統(tǒng)計(jì)，結(jié)果發(fā)現(xiàn)信號(hào)轉(zhuǎn)導(dǎo)機(jī)制所占基因數(shù)量最多，為1 318，其次是翻譯、核糖體結(jié)構(gòu)和生物發(fā)生，為1 236個(gè)基因，排第三的是翻譯后修飾、蛋白折疊和蛋白伴侶，為1 234個(gè)基因。

圖7 COG功能注釋分類

2.5.4 Nr注釋分類

將Unigene序列比對(duì)到Nr的蛋白序列，進(jìn)而得到每個(gè)基因的Nr注釋，通過Nr庫比對(duì)注釋的結(jié)果，統(tǒng)計(jì)并繪制比對(duì)上的物種分布圖(圖8)，只統(tǒng)計(jì)比對(duì)最多的前10個(gè)物種，其余的劃到Other species。從圖8中得出比對(duì)最多的前3個(gè)物種分別是栽培甜菜、藜麥和菠菜，說明馬齒莧和這3個(gè)物種蛋白序列相似度高，而劃到其他物種的有9 355個(gè)基因，暗示馬齒莧有自己特有的蛋白序列。且注釋的E-value值在0～1×10-150、1×10-150～1×10-100和1×10-100～1×10-50占比分別為25.65%、18.56%和28.77%，說明可信值較高(圖9)；一致性鑒定發(fā)現(xiàn)一致性非常高，一致性鑒定大于60%的基因數(shù)量為27 490個(gè)(圖10)。

圖8 Nr注釋

圖9 Nr注釋的值

圖10 Nr注釋的一致性分布

2.5.5 Pfam結(jié)構(gòu)域分類注釋

Pfam數(shù)據(jù)庫是一個(gè)蛋白質(zhì)家族大集合，根據(jù)多序列比對(duì)和隱馬爾可夫模型，利用hmmscan進(jìn)行結(jié)構(gòu)域預(yù)測(cè)，圖11顯示注釋最多的前20個(gè)域，其中蛋白激酶(pkinase)最多，有1 080個(gè)基因，其次是酪氨酸蛋白激酶(Pkinase_Tyr)，有1 027個(gè)基因。

圖11 Pfam結(jié)構(gòu)域分類注釋圖

2.5.6 SSR鑒定

圖12顯示unigene總共鑒定了25 898個(gè)SSR位點(diǎn)，其中單個(gè)堿基重復(fù)的最多，為12 772個(gè)，2個(gè)堿基重復(fù)的位列第二，為5 961個(gè)，排在第三的為3個(gè)堿基重復(fù)序列，是5 842個(gè)。

圖12 SSR重復(fù)單元數(shù)

2.6 差異表達(dá)轉(zhuǎn)錄本的篩選

根據(jù)差異表達(dá)分析軟件DESeq2篩選差異表達(dá)基因(DEGs)，篩選標(biāo)準(zhǔn)為：如果Qvalue<0.05，且log2FoldChange>2，則DEGs被認(rèn)為是顯著上調(diào)差異表達(dá)基因；如果Qvalue<0.05，且log2FoldChange<-2，則DEGs被認(rèn)為是顯著下調(diào)差異表達(dá)基因。差異基因火山圖可以顯示不同處理前后的差異基因分布，其中用紅色表示顯著上調(diào)的基因，以綠色表示顯著下調(diào)的基因，以黑色表示無顯著性差異的基因(圖13-A、13-B、13-C)。結(jié)果顯示，C vs DS之間鑒定出16 194個(gè)DEGs，其中包括2 874個(gè)上調(diào)DEGs，13 320個(gè)下調(diào)DEGs；C vs LTS之間獲得了16 093個(gè)差異表達(dá)基因，其中包括6 365個(gè)上調(diào)DEGs，9 728個(gè)下調(diào)基因；DS vs LTS之間發(fā)現(xiàn)了25 498個(gè)DEGs，其中包括15 163個(gè)上調(diào)DEGs，10 335個(gè)下調(diào)DEGs。通過維恩圖比較不同處理比較組的3個(gè)數(shù)據(jù)集(C vs DS、C vs LTS、DS vs LTS)，在3個(gè)比較組中，共有2 430個(gè)DEGs(圖13-D)，除3個(gè)比較組共有的DEGs外，在C vs DS和C vs LTS中還共有2 656個(gè)DEGs，C vs LTS和DS vs LTS中還共有7 919個(gè)DEGs，C vs DS和DS vs LTS比較組中共有8 195個(gè)DEGs，且每個(gè)比較組中都有自己特有的差異表達(dá)基因，C vs DS、C vs LTS、DS vs LTS中特有的差異表達(dá)基因分別有2 913、2 088和6 954個(gè)(圖13-D)。

2.7 差異表達(dá)基因GO功能分析

圖14分別列出前20個(gè)富集最顯著的差異表達(dá)基因，C vs DS差異基因數(shù)量最多的是赤霉素介導(dǎo)的信號(hào)轉(zhuǎn)導(dǎo)途徑，有77個(gè)；其次是一個(gè)氧原子結(jié)合的氧化還原酶活性，有66個(gè)；排在第三的是鈣離子介導(dǎo)的信號(hào)轉(zhuǎn)導(dǎo)途徑，有54個(gè)(圖14-A)；C vs LTS差異基因主要表現(xiàn)在細(xì)胞對(duì)各種光的反應(yīng)、葉綠素合成過程與丙酮酸磷酸二激酶活性等方面(圖14-B)；DS vs LTS差異基因數(shù)目排列前三的分別為丙酮酸磷酸二激酶活性(222個(gè))、木質(zhì)素生物合成過程(有127個(gè))、葉綠體核質(zhì)體(118個(gè))(圖14-C)。

圖14 C vs DS (A), C vs TS (B), DS vs LTS (C)差異表達(dá)基因GO分類圖

3 結(jié)論與討論

馬齒莧被世界衛(wèi)生組織列為最常用的藥用植物之一，是世界上第八大普通植物，被我國(guó)衛(wèi)生健康委員會(huì)劃定為78種藥食同源的野生植物之一，生態(tài)適應(yīng)能力強(qiáng)，抗旱抗?jié)?、耐鹽堿，在我國(guó)野生資源豐富，但未被充分重視和利用。因此，本文研究常溫(對(duì)照)、干旱脅迫和低溫脅迫下的野生馬齒莧葉片的全長(zhǎng)轉(zhuǎn)錄本，并進(jìn)行功能注釋和GO分類等研究，以期為加強(qiáng)馬齒莧的利用和其基因組的研究提供參考。根據(jù)測(cè)序共得到39 717條unigene，通過與Nr、Uniprot、Pfam、GO、KEGG、eggNOG和Pathway數(shù)據(jù)庫進(jìn)行同源比對(duì)發(fā)現(xiàn)，被注釋的基因?yàn)?6 381條，未被注釋的有3 336條，暗示馬齒莧有未知的轉(zhuǎn)錄本，可為馬齒莧屬植物新轉(zhuǎn)錄本的研究提供參考。本研究發(fā)現(xiàn)單個(gè)SSR位點(diǎn)數(shù)最多，為12 772個(gè)，2個(gè)位點(diǎn)和3個(gè)位點(diǎn)的SSR分別為5 961和5 842個(gè)。SSR標(biāo)記具有多態(tài)性高、信息量大等特點(diǎn)，常被用于物種的遺傳多樣研究、遺傳圖譜的繪制等，因此本轉(zhuǎn)錄組的SSR數(shù)據(jù)可以為后續(xù)馬齒莧屬植物的研究提供參考。差異表達(dá)基因篩選結(jié)果表明，干旱和低溫脅迫處理中上調(diào)和下調(diào)的差異表達(dá)基因數(shù)量較多，且不同比較組中除了有共有差異表達(dá)基因外，還有自己特有的差異表達(dá)基因，這些為以后研究這些基因的功能提供了參考依據(jù)。

A，C vs DS 火山圖；B，C vs LTS 火山圖；C，DSvsLTS 火山圖；D，不同處理組比較的維恩圖。

A， C vs DS Volcano； B， C vs LTS Volcano； C， 15 DSvsLTS Volcano； D， Venn diagram of DEGs in different comparisons.

圖13差異表達(dá)基因分布

Fig.13Distribution of differentially expressed genes