徐 文,安素妨,王 艷,賈琳琳,魯?shù)さ?張瑩瑩,劉建豐,李保全*

(1.河南省農(nóng)業(yè)科學(xué)院 作物設(shè)計(jì)中心，河南 鄭州 450002; 2.河南省農(nóng)業(yè)科學(xué)院 農(nóng)業(yè)經(jīng)濟(jì)與信息研究所，河南 鄭州 450002)

轉(zhuǎn)錄組是特定時(shí)間特定組織內(nèi)由基因轉(zhuǎn)錄產(chǎn)生的全部RNA轉(zhuǎn)錄本，包括編碼蛋白質(zhì)的RNA和非編碼RNA。轉(zhuǎn)錄組測序(RNA-seq)是借助于近年來發(fā)展迅速的高通量測序技術(shù)對特定細(xì)胞在某一功能狀態(tài)下幾乎所有轉(zhuǎn)錄出來的RNA進(jìn)行測序的一種技術(shù)方法[1]。隨著二代測序技術(shù)通量高、成本低、準(zhǔn)確度和靈敏度高等優(yōu)勢，RNA-seq近年來成為揭示基因轉(zhuǎn)錄和表達(dá)調(diào)控規(guī)律的有效手段，在植物發(fā)育、抗病抗逆、分子標(biāo)記開發(fā)和系統(tǒng)發(fā)生等方面得到了廣泛的應(yīng)用[2-6]。同時(shí)RNA-seq在基因結(jié)構(gòu)變異、低豐度轉(zhuǎn)錄本和未知轉(zhuǎn)錄本的檢測方面也發(fā)揮了重要的作用[7-13]。

甘藍(lán)型油菜(Brasscianapus)是世界四大油料作物之一，是我國重要的食用油來源和植物蛋白飼料資源之一，在國民經(jīng)濟(jì)和人們?nèi)粘Ｉ钪姓加兄匾匚?。甘藍(lán)型油菜基因組(AACC) 是白菜(B.rap)基因組(AA)與甘藍(lán)(B.oleracea)基因組(CC)雜交后經(jīng)過染色體加倍形成的[14]，其基因組參考序列于2014年公開發(fā)布，為油菜功能基因組學(xué)的研究奠定了良好的基礎(chǔ)[15]。測序發(fā)現(xiàn)，油菜基因組大小為849.7 Mb，由101 040個(gè)基因組成，但是由于其基因組極其復(fù)雜和注釋方法的局限性，油菜基因組還存在基因組組裝不完整、注釋信息不完善、新轉(zhuǎn)錄本的遺漏等問題。目前，利用RNA-seq研究甘藍(lán)型油菜發(fā)育、抗逆等的報(bào)道較多[16-18]，但是甘藍(lán)型油菜中新轉(zhuǎn)錄本的鑒定尚缺乏深入的研究。鑒于RNA-seq技術(shù)的優(yōu)點(diǎn)和甘藍(lán)型油菜基因組注釋的不完整性，利用RNA-seq和生物信息分析技術(shù)預(yù)測甘藍(lán)型油菜基因組中尚未注釋的轉(zhuǎn)錄本，并通過RT-PCR(Reverse transcription-PCR)和克隆測序的方法來驗(yàn)證，以期進(jìn)一步完善油菜基因組和轉(zhuǎn)錄組注釋信息，為進(jìn)一步挖掘油菜基因組中優(yōu)良的功能基因提供數(shù)據(jù)基礎(chǔ)。

1 材料和方法

1.1 試驗(yàn)材料

本研究以甘藍(lán)型油菜W系為試驗(yàn)材料。選取試驗(yàn)基地生長3周的油菜幼苗，取葉片于液氮中,置于-80 ℃ 冰箱中保存。

1.2 試驗(yàn)方法

1.2.1 RNA提取及檢測 油菜葉片RNA的提取采用Trizol法，按照說明書嚴(yán)格操作。利用Nanodrop檢測RNA的純度，Agilent 2100檢測RNA的完整性。

1.2.2 cDNA文庫構(gòu)建和測序 樣品RNA檢測合格后，進(jìn)行cDNA文庫構(gòu)建，主要流程如下：用帶有Oligo(dT)的beads富集mRNA；加入fragmentation buffer將mRNA打斷成短片段；以mRNA為模板合成第一鏈cDNA，后加入buffer、dNTPs和DNA Polymerase Ⅰ合成第二鏈cDNA；隨后利用AMPure XP beads純化雙鏈cDNA；純化的雙鏈cDNA再進(jìn)行末端修復(fù)、加A尾并連接測序接頭；然后用AMPure XP beads進(jìn)行片段大小選擇；最后通過PCR擴(kuò)增富集得到最終的cDNA文庫。

應(yīng)用Illumina HiSeq 2500 高通量測序平臺對構(gòu)建好的文庫進(jìn)行測序，測序讀長為雙端125 bp。

1.2.3 測序數(shù)據(jù)的質(zhì)控及比對 對測序得到的raw data，去除含有adapter的reads、去除N堿基比例大于5%的reads、去除低質(zhì)量的reads得到clean data。從http://plants.ensembl.org/index.html下載已公布的甘藍(lán)型油菜參考基因組序列，用Bowtie 2軟件對油菜參考基因組序列建立索引[19]，用Tophat 2軟件將得到的clean data與參考基因組序列進(jìn)行mapping[20],mapped reads以bam格式文件輸出。

1.2.4 轉(zhuǎn)錄本的表達(dá)豐度計(jì)算和新轉(zhuǎn)錄本的鑒定 轉(zhuǎn)錄本的豐度通過計(jì)算FPKM(Fragments per kilobase per million mapped reads)值來度量。通過Cufflinks軟件包中的cuffdiff命令來實(shí)現(xiàn)對每個(gè)樣品表達(dá)豐度的計(jì)算[21]。

將每個(gè)樣品比對得到的bam文件，借助Cufflinks軟件與甘藍(lán)型油菜參考基因組的注釋文件進(jìn)行位置信息的比較和整合，初步建立每個(gè)樣品的轉(zhuǎn)錄本文庫；通過cuffmerge命令的合并功能對得到的每個(gè)樣品的轉(zhuǎn)錄本文庫進(jìn)行整合，得到1個(gè)完整的基因組注釋文件merged.gtf；然后用cuffcompare命令與甘藍(lán)型油菜已知的轉(zhuǎn)錄本信息進(jìn)行比較，從而鑒定候選的新轉(zhuǎn)錄本。

1.2.5 新轉(zhuǎn)錄本的RT-PCR和測序驗(yàn)證 從鑒定的候選轉(zhuǎn)錄本中選取一部分，以Primer Premier 5.0根據(jù)其序列設(shè)計(jì)引物(表1)；提取甘藍(lán)型油菜葉片RNA，反轉(zhuǎn)錄成cDNA，進(jìn)行PCR擴(kuò)增，將擴(kuò)增產(chǎn)物用瓊脂糖電泳檢測，然后純化回收，連接T載體，轉(zhuǎn)化DH5α感受態(tài)細(xì)胞，挑選陽性克隆，送公司測序。

表1 甘藍(lán)型油菜新轉(zhuǎn)錄本及其PCR擴(kuò)增引物

續(xù)表1 甘藍(lán)型油菜新轉(zhuǎn)錄本及其PCR擴(kuò)增引物

2 結(jié)果與分析

2.1 甘藍(lán)型油菜RNA質(zhì)量的檢測

經(jīng)檢測，17個(gè)甘藍(lán)型油菜葉片RNA均符合cDNA文庫構(gòu)建的標(biāo)準(zhǔn)：OD260/OD280為1.8～2.2、28S rRNA/18S rRNA≥1.5、OD260/OD230≥1.96、RIN值≥9.7，表明RNA純度和完整性較好。

2.2 甘藍(lán)型油菜RNA測序數(shù)據(jù)處理和分析

經(jīng)過質(zhì)量控制后總計(jì)得到848 866 766條clean reads，與甘藍(lán)型油菜參考基因組序列進(jìn)行比對。從表2可以看出，各樣品的reads與其參考基因組序列的比對率在65.00%～85.40%，說明測序數(shù)據(jù)的比對率正常。

表2 clean data與參考基因組序列比對結(jié)果

2.3 甘藍(lán)型油菜新轉(zhuǎn)錄本的識別

通過Tophat/Cufflinks一系列流程的分析，將比對成功的序列進(jìn)行組裝整合。然后與已知甘藍(lán)型油菜參考基因組注釋轉(zhuǎn)錄本信息進(jìn)行比較，最終獲得了由612 085個(gè)外顯子與467 743個(gè)內(nèi)含子組成的103 310個(gè)基因位點(diǎn)(包含158 004個(gè)mRNA)，其中137 756個(gè)轉(zhuǎn)錄本是有多個(gè)外顯子組成的。其中，甘藍(lán)型油菜已知注釋的101 040個(gè)基因位點(diǎn)全部包含在內(nèi)。新鑒定的外顯子有33 811個(gè)，新鑒定的內(nèi)含子有26 839個(gè)。

對于鑒定到的158 004個(gè)轉(zhuǎn)錄本，可劃分為6類(表3):與內(nèi)含子鏈匹配的轉(zhuǎn)錄本; 潛在的新轉(zhuǎn)錄本; 與已知外顯子重疊的轉(zhuǎn)錄本;內(nèi)含子與反義鏈上已知內(nèi)含子重疊的轉(zhuǎn)錄本; 未知的基因間隔區(qū)轉(zhuǎn)錄本;外顯子與反義鏈上已知轉(zhuǎn)錄本重疊的轉(zhuǎn)錄本。 其中，未獲得注釋的7 720個(gè)轉(zhuǎn)錄本被劃分為未知的基因間隔區(qū)轉(zhuǎn)錄本。

表3 新鑒定的甘藍(lán)型油菜轉(zhuǎn)錄本和已知轉(zhuǎn)錄本的比較結(jié)果

2.4 甘藍(lán)型油菜新轉(zhuǎn)錄本的鑒定及驗(yàn)證

根據(jù)轉(zhuǎn)RNA-seq預(yù)測的轉(zhuǎn)錄本單元及其表達(dá)豐度，本研究選取了18個(gè)平均表達(dá)豐度在100以上的新鑒定的轉(zhuǎn)錄本進(jìn)行RT-PCR擴(kuò)增并測序驗(yàn)證，結(jié)果發(fā)現(xiàn)15個(gè)新轉(zhuǎn)錄本擴(kuò)增條帶很清晰(圖1)。

為了驗(yàn)證結(jié)果的準(zhǔn)確性和可靠性，將15個(gè)轉(zhuǎn)錄本PCR擴(kuò)增產(chǎn)物進(jìn)行克隆,挑選轉(zhuǎn)化子送公司測序，結(jié)果發(fā)現(xiàn)，測序得到的序列和轉(zhuǎn)錄組數(shù)據(jù)分析組裝出來的序列是一致的(圖2)。

M：Marker；1：XLOC_028977； 2：XLOC_056607；3：XLOC_009511；4：XLOC_042556；5：XLOC_097486；6：XLOC_095366；7:XLOC_097244；8：XLOC_017634；9：XLOC_015102；10：XLOC_066937；11：XLOC_092505；12：XLOC_073270；13：XLOC_046341；14：XLOC_100972；15：XLOC_041903

圖1新轉(zhuǎn)錄本RT-PCR擴(kuò)增電泳檢測結(jié)果

Sequencing表示PCR產(chǎn)物測序序列，XLOC_095366表示通過轉(zhuǎn)錄組分析預(yù)測的轉(zhuǎn)錄本序列圖2 XLOC_095366轉(zhuǎn)錄本比對結(jié)果

3 結(jié)論與討論

本研究利用RNA-seq和生物信息學(xué)技術(shù)對甘藍(lán)型油菜葉片轉(zhuǎn)錄組文庫中的新轉(zhuǎn)錄本進(jìn)行了挖掘和鑒定，并采用RT-PCR和克隆測序?qū)﹁b定到的新轉(zhuǎn)錄本進(jìn)行了驗(yàn)證。研究結(jié)果表明，在甘藍(lán)型油菜基因組已知基因的間隔區(qū)發(fā)現(xiàn)了7 720個(gè)新轉(zhuǎn)錄本，通過RT-PCR擴(kuò)增驗(yàn)證了15個(gè)新鑒定的轉(zhuǎn)錄本，為油菜基因組提供了補(bǔ)充性的注釋信息，為進(jìn)一步挖掘甘藍(lán)型油菜基因組中優(yōu)良的功能基因提供數(shù)據(jù)基礎(chǔ)。

同一個(gè)基因通過可變剪接后形成多種mRNA成熟體即不同的轉(zhuǎn)錄本,經(jīng)翻譯產(chǎn)生結(jié)構(gòu)和功能不同的蛋白質(zhì)或者具有調(diào)控功能的非編碼RNA。從PCR擴(kuò)增檢測的電泳圖中可以看出，新鑒定的轉(zhuǎn)錄本有部分?jǐn)U增的條帶不是單一的，說明這些轉(zhuǎn)錄本可能存在不同的可變剪接體。那么需要進(jìn)一步深入研究這些轉(zhuǎn)錄本的特點(diǎn)，區(qū)分每個(gè)轉(zhuǎn)錄本是編碼蛋白質(zhì)的基因還是具有調(diào)控功能的非編碼RNA。對于具有編碼功能的轉(zhuǎn)錄本，要深入研究其功能并挖掘可能存在的不同剪接體。對于在植物生長發(fā)育、逆境脅迫等方面發(fā)揮著重要作用的非編碼RNA，要克隆這些非編碼RNA并研究它們對功能基因的調(diào)控機(jī)制。因此，深入挖掘和研究這些新轉(zhuǎn)錄本及其作用機(jī)制對完善甘藍(lán)型油菜基因組注釋、為分子育種提供豐富的基因信息具有重要的理論意義。

參考文獻(xiàn):

[1] Wang Z,Gerstein M,Snyder M.RNA-Seq:A revolutionary tool for transcriptomics[J].Nature Reviews Genetics,2009,10(1):57-63.

[2] Schlueter U,Denton A K,Braeutigam A.Understanding metabolite transport and metabolism in C4 plants through RNA-seq[J].Current Opinion in Plant Biology,2016,31:83-90.

[3] Martin L B,Fei Z,Giovannoni J J,etal.Catalyzing plant science research with RNA-seq[J].Frontiers in Plant Science,2013,4:66.

[4] Kakumanu A,Ambavaram M M,Klumas C A,etal.Effects of drought on gene expression in maize reproductive and leaf meristem tissue revealed by RNA-Seq[J].Plant Physiology,2012,160(2):846-867.

[5] Wei L,Li S,Liu S,etal.Transcriptome analysis ofHouttuyniacordataThunb.by Illumina paired-end RNA sequencing and SSR marker discovery[J].PLoS One,2014,9(1):e84105.

[6] Wei Z,Sun Z,Cui B,etal.Transcriptome analysis of colored calla lily (ZantedeschiarehmanniiEngl.) by Illumina sequencing:De novo assembly,annotation and EST-SSR marker development[J].Peer J,2016,4:e2378.

[7] Chettoor A M,Givan S A,Cole R A,etal.Discovery of novel transcripts and gametophytic functions via RNA-seq analysis of maize gametophytic transcriptomes[J].Genome Biology,2014,15:414.

[8] Alasoo K,Martinez F O,Hale C,etal.Transcriptional profiling of macrophages derived from monocytes and iPS cells identifies a conserved response to LPS and novel alternative transcription[J].Scientific Reports,2015,5:12524.

[9] Du Z Q,Eisley C J,Onteru S K,etal.Identification of species-specific novel transcripts in pig reproductive tissues using RNA-seq[J].Animal Genetics,2014,45(2):198-204.

[10] Roberts A,Pimentel H,Trapnell C,etal.Identification of novel transcripts in annotated genomes using RNA-Seq[J].Bioinformatics,2011,27(17):2325-2329.

[11] Weirick T,Militello G,Mueller R,etal.The identification and characterization of novel transcripts from RNA-seq data[J].Briefings in Bioinformatics,2016,17(4):678-685.

[12] Zhao J,Song X,Wang K.lncScore:Alignment-free identification of long noncoding RNA from assembled novel transcripts[J].Scientific Reports,2016,6:34838.

[13] Weikard R,Hadlich F,Kuehn C.Identification of novel transcripts and noncoding RNAs in bovine skin by deep next generation sequencing[J].BMC Genomics,2013,14:789.

[14] Allender C J,King G J.Origins of the amphiploid speciesBrassicanapusL.investigated by chloroplast and nuclear molecular markers[J].BMC Plant Biology,2010,10:54.

[15] Chalhoub B,Denoeud F,Liu S,etal.Early allopolyploid evolution in the post-neolithicBrassicanapusoilseed genome[J].Science,2014,345:950-953.

[16] Yan X,Dong C,Yu J,etal.Transcriptome profile analysis of young floral buds of fertile and sterile plants from the self-pollinated offspring of the hybrid between novel restorer line NR1 and Nsa CMS line inBrassicanapus[J].BMC Genomics,2013,14:26.

[17] Gill R A,Ali B,Cui P,etal.Comparative transcriptome profiling of twoBrassicanapuscultivars under chromium toxicity and its alleviation by reduced glutathione[J].BMC Genomics,2016,17:885.

[18] Lu K,Peng L,Zhang C,etal.Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits inBrassicanapus[J].Frontiers in Plant Science,2017,8:206.

[19] Langmead B,Salzberg S L.Fast gapped-read alignment with Bowtie 2[J].Nature Methods,2012,9(4):357-359.

[20] Kim D,Pertea G,Trapnell C,etal.TopHat2:Accurate alignment of transcriptomes in the presence of insertions,deletions and gene fusions[J].Genome Biology,2013,14(4):R36.

[21] Trapnell C,Roberts A,Goff L,etal.Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks[J].Nature Protocols,2012,7(3):562-578.