基于RNA-seq技术的甘蓝型油菜新转录本的鉴定

徐 文,安素妨,王 艳,贾琳琳,鲁丹丹,张莹莹,刘建丰,李保全*

(1.河南省农业科学院 作物设计中心，河南 郑州 450002; 2.河南省农业科学院 农业经济与信息研究所，河南 郑州 450002)

转录组是特定时间特定组织内由基因转录产生的全部RNA转录本，包括编码蛋白质的RNA和非编码RNA。转录组测序(RNA-seq)是借助于近年来发展迅速的高通量测序技术对特定细胞在某一功能状态下几乎所有转录出来的RNA进行测序的一种技术方法[1]。随着二代测序技术通量高、成本低、准确度和灵敏度高等优势，RNA-seq近年来成为揭示基因转录和表达调控规律的有效手段，在植物发育、抗病抗逆、分子标记开发和系统发生等方面得到了广泛的应用[2-6]。同时RNA-seq在基因结构变异、低丰度转录本和未知转录本的检测方面也发挥了重要的作用[7-13]。

甘蓝型油菜(Brasscianapus)是世界四大油料作物之一，是我国重要的食用油来源和植物蛋白饲料资源之一，在国民经济和人们日常生活中占有重要地位。甘蓝型油菜基因组(AACC) 是白菜(B.rap)基因组(AA)与甘蓝(B.oleracea)基因组(CC)杂交后经过染色体加倍形成的[14]，其基因组参考序列于2014年公开发布，为油菜功能基因组学的研究奠定了良好的基础[15]。测序发现，油菜基因组大小为849.7 Mb，由101 040个基因组成，但是由于其基因组极其复杂和注释方法的局限性，油菜基因组还存在基因组组装不完整、注释信息不完善、新转录本的遗漏等问题。目前，利用RNA-seq研究甘蓝型油菜发育、抗逆等的报道较多[16-18]，但是甘蓝型油菜中新转录本的鉴定尚缺乏深入的研究。鉴于RNA-seq技术的优点和甘蓝型油菜基因组注释的不完整性，利用RNA-seq和生物信息分析技术预测甘蓝型油菜基因组中尚未注释的转录本，并通过RT-PCR(Reverse transcription-PCR)和克隆测序的方法来验证，以期进一步完善油菜基因组和转录组注释信息，为进一步挖掘油菜基因组中优良的功能基因提供数据基础。

1 材料和方法

1.1 试验材料

本研究以甘蓝型油菜W系为试验材料。选取试验基地生长3周的油菜幼苗，取叶片于液氮中,置于-80 ℃ 冰箱中保存。

1.2 试验方法

1.2.1 RNA提取及检测 油菜叶片RNA的提取采用Trizol法，按照说明书严格操作。利用Nanodrop检测RNA的纯度，Agilent 2100检测RNA的完整性。

1.2.2 cDNA文库构建和测序 样品RNA检测合格后，进行cDNA文库构建，主要流程如下：用带有Oligo(dT)的beads富集mRNA；加入fragmentation buffer将mRNA打断成短片段；以mRNA为模板合成第一链cDNA，后加入buffer、dNTPs和DNA Polymerase Ⅰ合成第二链cDNA；随后利用AMPure XP beads纯化双链cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头；然后用AMPure XP beads进行片段大小选择；最后通过PCR扩增富集得到最终的cDNA文库。

应用Illumina HiSeq 2500 高通量测序平台对构建好的文库进行测序，测序读长为双端125 bp。

1.2.3 测序数据的质控及比对 对测序得到的raw data，去除含有adapter的reads、去除N碱基比例大于5%的reads、去除低质量的reads得到clean data。从http://plants.ensembl.org/index.html下载已公布的甘蓝型油菜参考基因组序列，用Bowtie 2软件对油菜参考基因组序列建立索引[19]，用Tophat 2软件将得到的clean data与参考基因组序列进行mapping[20],mapped reads以bam格式文件输出。

1.2.4 转录本的表达丰度计算和新转录本的鉴定 转录本的丰度通过计算FPKM(Fragments per kilobase per million mapped reads)值来度量。通过Cufflinks软件包中的cuffdiff命令来实现对每个样品表达丰度的计算[21]。

将每个样品比对得到的bam文件，借助Cufflinks软件与甘蓝型油菜参考基因组的注释文件进行位置信息的比较和整合，初步建立每个样品的转录本文库；通过cuffmerge命令的合并功能对得到的每个样品的转录本文库进行整合，得到1个完整的基因组注释文件merged.gtf；然后用cuffcompare命令与甘蓝型油菜已知的转录本信息进行比较，从而鉴定候选的新转录本。

1.2.5 新转录本的RT-PCR和测序验证 从鉴定的候选转录本中选取一部分，以Primer Premier 5.0根据其序列设计引物(表1)；提取甘蓝型油菜叶片RNA，反转录成cDNA，进行PCR扩增，将扩增产物用琼脂糖电泳检测，然后纯化回收，连接T载体，转化DH5α感受态细胞，挑选阳性克隆，送公司测序。

表1 甘蓝型油菜新转录本及其PCR扩增引物

续表1 甘蓝型油菜新转录本及其PCR扩增引物

2 结果与分析

2.1 甘蓝型油菜RNA质量的检测

经检测，17个甘蓝型油菜叶片RNA均符合cDNA文库构建的标准：OD260/OD280为1.8～2.2、28S rRNA/18S rRNA≥1.5、OD260/OD230≥1.96、RIN值≥9.7，表明RNA纯度和完整性较好。

2.2 甘蓝型油菜RNA测序数据处理和分析

经过质量控制后总计得到848 866 766条clean reads，与甘蓝型油菜参考基因组序列进行比对。从表2可以看出，各样品的reads与其参考基因组序列的比对率在65.00%～85.40%，说明测序数据的比对率正常。

表2 clean data与参考基因组序列比对结果

2.3 甘蓝型油菜新转录本的识别

通过Tophat/Cufflinks一系列流程的分析，将比对成功的序列进行组装整合。然后与已知甘蓝型油菜参考基因组注释转录本信息进行比较，最终获得了由612 085个外显子与467 743个内含子组成的103 310个基因位点(包含158 004个mRNA)，其中137 756个转录本是有多个外显子组成的。其中，甘蓝型油菜已知注释的101 040个基因位点全部包含在内。新鉴定的外显子有33 811个，新鉴定的内含子有26 839个。

对于鉴定到的158 004个转录本，可划分为6类(表3):与内含子链匹配的转录本; 潜在的新转录本; 与已知外显子重叠的转录本;内含子与反义链上已知内含子重叠的转录本; 未知的基因间隔区转录本;外显子与反义链上已知转录本重叠的转录本。 其中，未获得注释的7 720个转录本被划分为未知的基因间隔区转录本。

表3 新鉴定的甘蓝型油菜转录本和已知转录本的比较结果

2.4 甘蓝型油菜新转录本的鉴定及验证

根据转RNA-seq预测的转录本单元及其表达丰度，本研究选取了18个平均表达丰度在100以上的新鉴定的转录本进行RT-PCR扩增并测序验证，结果发现15个新转录本扩增条带很清晰(图1)。

为了验证结果的准确性和可靠性，将15个转录本PCR扩增产物进行克隆,挑选转化子送公司测序，结果发现，测序得到的序列和转录组数据分析组装出来的序列是一致的(图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新转录本RT-PCR扩增电泳检测结果

Sequencing表示PCR产物测序序列，XLOC_095366表示通过转录组分析预测的转录本序列图2 XLOC_095366转录本比对结果

3 结论与讨论

本研究利用RNA-seq和生物信息学技术对甘蓝型油菜叶片转录组文库中的新转录本进行了挖掘和鉴定，并采用RT-PCR和克隆测序对鉴定到的新转录本进行了验证。研究结果表明，在甘蓝型油菜基因组已知基因的间隔区发现了7 720个新转录本，通过RT-PCR扩增验证了15个新鉴定的转录本，为油菜基因组提供了补充性的注释信息，为进一步挖掘甘蓝型油菜基因组中优良的功能基因提供数据基础。

同一个基因通过可变剪接后形成多种mRNA成熟体即不同的转录本,经翻译产生结构和功能不同的蛋白质或者具有调控功能的非编码RNA。从PCR扩增检测的电泳图中可以看出，新鉴定的转录本有部分扩增的条带不是单一的，说明这些转录本可能存在不同的可变剪接体。那么需要进一步深入研究这些转录本的特点，区分每个转录本是编码蛋白质的基因还是具有调控功能的非编码RNA。对于具有编码功能的转录本，要深入研究其功能并挖掘可能存在的不同剪接体。对于在植物生长发育、逆境胁迫等方面发挥着重要作用的非编码RNA，要克隆这些非编码RNA并研究它们对功能基因的调控机制。因此，深入挖掘和研究这些新转录本及其作用机制对完善甘蓝型油菜基因组注释、为分子育种提供丰富的基因信息具有重要的理论意义。

参考文献:

[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.