RNA-seq技术在野生动物研究中的应用

肖 慧 吴 尽 赵敏蝶 丁 新 刘学东

(东北林业大学，哈尔滨，150040)

1 RNA-seq技术的发展历史及其在生物学研究中的作用

1964年Holley第一次成功获得一个完整基因的核苷酸序列［6］，此后核酸测序方法就不断地快速发展。随着高通量测序时代的到来，大规模并行测序(massive parallel sequencing，MPS)平台如Roche公司(454 GS-FLX)、Illumina公司(Genome Analyzer II)和 ABI公司(ABSOLiD)彻底变革了测序技术，也改变了转录组的研究方法，产生了RNA测序技术(RNA-seq)。

转录组研究能够从整体水平研究基因功能以及基因结构，揭示特定生物学过程以及疾病发生过程中的分子机理。RNA-seq作为一种新的高效、快捷的研究手段正在改变着人们对转录组的认识。在生物学研究中，RNA-seq可以进行转录本结构研究(基因边界鉴定、可变剪切研究等)、转录本结构变异研究(如基因融合、编码区 SNP研究)、非编码区域功能研究(Non-coding RNA、microRNA前体研究等)、基因表达水平研究以及全新转录本的发现［7］。这些研究为人类认识和了解生物机体和各种功能提供了重要的方法和途径。同样，RNA-seq在野生动物的研究中也得到推广和使用。

2 RNA-seq在野生动物研究中的方法

近年来，随着转录组测序技术的成熟与发展，许多研究人员利用该技术对各种野生动物的转录本进行研究，建立了较为全面的相关野生动物转录组数据库［8－10］，为进一步研究野生动物生长过程中的基因结构和功能的变化及代谢通路的调控等提供实验依据。研究的具体方法通常分为以下几部分。

2.1 样品分离和文库制备

首先，研究者应根据研究对象及研究目的，选择性地进行样品采集和处理。考虑到动物福利，尽量做到低损伤取样甚至无损伤取样，杨晓光［8］为了建立马鹿(Cervus elaphus)鹿茸増生区茸皮和软骨转录组数据库，分别采集3头3～4岁雄性东北马鹿的生长期鹿茸样本，单独提取每个样本中的总RNA后再进行混合，以求尽可能囊括马鹿增生区所有基因。杨秀峰［10］为了获取狼(Canis lupus)和家犬(Canis lupus familiaris)的转录组数据，选取狼和家犬的血液样本进行转录组测序。以上样品相对于动物其他组织器官来说更容易获取，对动物损伤较低且可恢复。在样本采集之后，利用生物公司提供的试剂盒进行总RNA提取，随后将总RNA片段化，连接到特定的接头序列并反转录，得到的cDNA片段进行克隆性扩增，制备文库以备高通量测序使用。

2.2 高通量测序

在上述得到的cDNA片段两端加上接头，使用二代高通量测序仪测序得到足够的多的reads序列。目前各测序平台广泛采用的一种测序方法是双末端测序法，该方法对文库中的每一个cDNA分子的两端均进行序列测定，这样每个cDNA分子就可获得两个经过测序的序列片段(reads)，相对于单端测序增加了物理覆盖度［11－12］，通过对这两个序列片段进行标记就可以在测序得到的数据中识别一对双末端配对序列，显著增强了对数据分析的能力。通常情况下，测序深度在10×～15×以上时覆盖度和测序错误率控制均得以保证。在转录组测序中，双末端测序不仅解决了测序reads不够长的问题，还能发现新的结构变异，区别不同的剪接体，鉴定由染色体重组造成的融合基因［11－14］等问题。

2.3 转录组装配

接下来需要进行转录组的装配。模式生物常具有参考基因组DNA序列信息，可以对测序结果进行基于参考基因组的比对和拼接，再通过进行基因组定位和注释来获得基因组尺度的转录图谱。这种方法虽然方便，但其明显缺陷在于新转录序列的丢失。而非模式生物缺乏参考基因组，只能自体组装(De novo assembly)，这时可以使用一些软件如Velvet与ABySS自体组装表达序列标签(EST)，或在近缘生物数据的帮助下引导装配［15－16］。为提高装配质量，策略之一是尽量增加reads的覆盖度，以及混合使用不同类型的reads。

2.4 生物信息学分析

转录组拼接后得到转录本，为了获取转录本中的信息，接下来需要对转录本进行功能注释和表达定量分析。

2.4.1 转录组功能注释及代谢通路分析

②依法划定饮用水水源性坑塘保护区，禁止一切直接或间接污染水体的行为。制定防止污染和人为破坏的管理办法，从源头上保证水源的可持续利用和农村经济可持续发展。

研究人员一般会用Blast程序对转录组功能进行注释［17］。常用的参考序列数据库包括NCBI中的非冗余核酸数据库(non-redundant protein database，nr database)、Swiss-Prot等。除了对每条转录本进行功能注释之外，对转录组的注释还有COG(cluster of orthologous group)注释［18］、GO(gene ontology) 注释以及PATHWAY注释。COG注释可以根据基因功能和进化关系对转录组中的序列进行分类，进而可以宏观地认识和比较物种的转录组构成。GO注释是从分子功能、细胞组成以及生物学过程3个方面对基因进行注释。Pathway指代谢通路，Pathway注释主要指构建转录组包含的生物代谢通路和调控关系网络。目前Pathway分析主要有 KEGG［19］和 BioCyc［20］。以 KEGG 为例，KEGG是系统分析基因产物在细胞中的代谢途径以及这些基因产物的功能的数据库，利用KEGG可以进一步研究基因网络在生物学上的复杂功能。

2.4.2 表达定量与表达差异分析

以注释和read映射为基础，基于浏览器对于数据质量可视化和特定事件进行解释非常重要。不过，它们只提供了有关研究的质量画面，大量数据及其相关细节并不能简单地通过这种方式表现出来。因此，大部分RNA-seq第二阶段的内容是关于全基因组转录事件的自动定量研究。其研究内容包括定量已知元件(即，已注释的基因或外显子)，与检测尚未在数据库中注释为外显子的新转录区。通常，定量步骤是进行任何差异表达研究的基础。RNA-seq的基因表达分析技术是基于对reads的计数，对低表达的基因也能够检测，具有灵敏度高、分辨率高、无饱和区等优点［3，21－22］。考虑到样本大小的影响(如 reads 数量随基因长度不同而不同;测序深度不同，测序获得的reads总数也不同)，Mortazavi提出了RPKM(reads per kilobases per million reads)和FPKM(fragments per kilobase per million reads)作为标准化定量指数来消除这两种系统差异。经过这种归一化处理，不同长度、不同测序深度下的基因表达量具有可比性［23］。另外可以在KEGG通路上进行富集分析，之后再进行表达差异分析。

基因表达差异分析是指找出不同时间点、不同组织或者不同处理条件下具有差异表达的基因。而转录组研究的最终目的就是要定量多个样本的表达来获得差异基因表达，确定样本特异性可变剪接异构体和它们差异性丰度。为了分析野生动物基因表达与作用因素之间的关系，可以用统计学的方法对高通量测序得到的基因表达量进行分析比对，找出样本间差异性显著的基因，再进行进一步的分析研究。研究人员通常会采用假设检验算法(P值)对样本之间的差异基因进行筛选，然后再通过错误发现率(false discovery rate，FDR)方法来校准P值。符合以下标准的基因被认为存在基因表达差异:当比对某基因在两样本间的FDR值＜0.001，且log2Ratio值＞1时可认为基因存在表达差异。

3 RNA-seq在野生动物研究中的应用

3.1 构建野生动物基因蓝图

RNA-seq最大的优点是不局限于检测已知基因组序列的转录组，还可以用于全基因组序列未知的非模式生物的转录组测序分析，这一点极大地拓宽了转录组学的研究对象范围，为研究人员提供更多的研究思路。2008年Vera等［24］使用454 GS-20测序技术对非模式生物——庆网蛱蝶(Melitaea cinxia)进行了转录组测序研究，这是第一个利用自体组装进行转录组研究的范例。随后包括野生动物在内的大量生物转录组测序研究不断出现，RNA-seq技术极大地推动野生动物转录组研究，这其中包括国外研究的鲟鱼(Acipenser fulvescens)［25］、虹鳟 (Oncorhynchus mykiss)［26］、响尾蛇(Crotalus adamanteus)［27］、食蟹猴(Macaca fascicu-laris)［28］和国内研究的斑海豹(Phoca largha)［29］、绒山羊［30］、马鹿［8］和狼［10］等。

Hale等通过RNA-seq技术成功获得了鲟鱼的转录组序列，确定超过5000个表达序列标签并鉴定877个候选SNP，大约每460个碱基就有一个SNP。杨秀峰应用RNA-seq技术对2只家犬和3只狼的血液转录组进行测序，该研究通过Illumina HiSeqTM2000平台进行测序，组装后获得了26212个基因，其中新基因1989个;总共鉴定出33229个转录本，其中新转录本1993个。这些研究为构建野生动物基因的蓝图增添了基础数据。

3.2 差异表达分析

RNA-seq的另一个优势是它可以捕捉不同组织或状态下的转录组动态变化，通过分析不同因素作用下的基因在RNA水平表达差异性，可以将那些显著差异表达的基因与某些生物学功能联系起来。例如，通过RNA-seq技术对狼和家犬的血液转录组进行分析，GO富集分析发现6个在狼中高表达的基因富集到了狼的先天性免疫系统上，推测这可能与狼在某些病毒的抵抗能力上大于家犬相关［10］。另外，我们研究组对马鹿鹿茸角快速生长期鹿茸生长点茸皮和软骨进行高通量测序，发现茸皮和软骨各自特异性表达6961和2776条基因，通过GO功能富集分析及KEGG通路富集分析，这些差异表达基因主要参与细胞结构、细胞代谢、蛋白质相互作用、催化活性等生物学进程。其中涉及注释了5328个基因的236条通路，这些基因和通路对茸皮和软骨组织特异性生长发育具有重要的调控作用［8］。这些差异基因的发现为进一步研究野生动物体内的分子生物学机制奠定基础。

3.3 非编码区域功能研究

非编码RNA(ncRNA)指的是未被翻译成蛋白质的RNA分子，重要的非编码RNA有转运RNA(tRNAs)、核糖体RNA(rRNAs)以及小RNA如microRNAs(miRNA)、siRNAs等，在一系列与细胞存活有关的活动中发挥重要功能［31］。Berezikov［32］等使用大规模平行测序来比较人类和黑猩猩(Pan troglodytes)大脑的microRNA含量，发现了447个新的miRNA基因，其中许多新的miRNA在灵长类之间并不保守，以此为基础探讨miRNA的进化过程以及人类与黑猩猩的大脑进化和功能的差异。陈艳霞等通过Solexa高通量测序对鹿茸软骨和茸皮组织的miRNA进行了全面的鉴定与特征分析，首次通过同源比对鉴定了鹿茸软骨和茸皮miRNA共684个，其中611个哺乳动物保守miRNA和73个鹿茸新的候选miRNA。通过对这些miRNA靶基因功能注释，发现在快速生长期鹿茸软骨和茸皮中的很多基因参与细胞或细胞器构成、核酸与蛋白质生物合成、催化活性、代谢过程、细胞増殖、配体/受体相互作用及多种信号通路。这些基因和通路在鹿茸快速生长发育过程中起到重要调控作用［9］。

4 展望

自RNA-seq技术问世以来，已经成为生物研究中不可或缺的技术手段，并且随着高通量测序的快速发展，测序成本的不断降低，转录学研究逐渐成为研究热点。RNA-seq技术在野生动物研究中的应用已有数年，NGS测序技术的快速发展又为野生动物转录组学的研究提供了必要的技术手段。尤其对于那些没有基因组序列信息的野生动物，利用RNA-seq技术可以通过比较基因组学对其进行基因组数据注释，这极大拓宽了RNA-seq技术在野生动物领域的研究，在后续的野生动物研究中RNA-seq将发挥越来越大的作用。

［1］ Chu Yongjun，Corey D R．RNA sequencing:platform selection，experimental design，and data interpretation ［J］．Nucleic Acid Therapeutics，2012，22(4):271 －274．

［2］ Maher C A，Kumar-Sinha C，Cao Xuhong A，et al．Transcriptome sequencing to detect gene fusions in cancer［J］．Nature，2009，458(7234):97－101．

［3］ Wang Zhong，Gerstein M，Snyder M．RNA-seq:a revolutionary tool for transcriptomics ［J］．Nature Reviews Genetics，2009，10(1):57－63．

［4］ 刘红亮，郑丽明，刘青青，等．非模式生物转录组研究 ［J］．遗传，2013，35(8):955－970．

［5］ Ingolia N T，Brar G A，Rouskin S A，et al．The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments ［J］．Nature Protocols，2012，7(8):1534－1550．

［6］ Holley R W．Alanine transfer RNA，in Nobel lectures in Adokctdar biology 1933－1975［S］．Elsevier North Holland:New York，NY，USA．1977:285－300．

［7］ 祁云霞，刘永斌，荣威恒．转录组研究新技术:RNA-seq及其应用 ［J］．遗传，2011，33(11):1191－1202．

［8］ 杨晓光．马鹿(Cervus elaphus)鹿茸角顶端茸皮与软骨组织转录组研究［D］．哈尔滨:东北林业大学，2015．

［9］ 陈艳霞．马鹿(Cervus elaphus)鹿茸快速生长期生长点软骨和茸皮组织microRNA表达谱研究 ［D］．哈尔滨:东北林业大学，2015．

［10］ 杨秀峰．基于高通量测序的狼和家犬血液转录组研究［D］．曲阜:曲阜师范大学，2016．

［11］ Ozsolak F，Milos P M．RNA sequencing:advances，challenges and opportunities［J］．Nature Reviews Genetics，2011，12(2):87－98．

［12］ Maher C A，Palanisamy N，Brenner J C，et al．Chimeric transcript discovery by paired-end transcriptome sequencing［J］．Proceedings of the National Academy of Sciences of the United States of America，2009，106(30):12353－12358．

［13］ Au K F，Jiang Hui，Lin Lan，et al．Detection of splice junctions from paired-end RNA-seq data by SpliceMap［J］．Nucleic Acids Research，2010，38(14):4570－4578．

［14］ Edgren H，Murumagi A，Kangaspeska S，et al．Identification of fusion genes in breast cancer by paired-end RNA-sequencing ［J］．Genome Biology，2011，12(1):R6．

［15］ Birol I，Jackman S D，Nielsen C B，et al．De novo transcriptome assembly with ABySS ［J］．Bioinformatics，2009，25(21):2872－2877．

［16］ Zerbino D R，Birney E．Velvet:algorithms for de novo short read assembly using de bruijn graphs［J］．Genome Research，2008，18(5):821－829．

［17］ Altschul S F，Gish W，Miller W，et al．Basic local alignment search tool［J］．Journal of Molecular Biology，1990，215(3):403－410．

［18］ Tatusov R L，Fedorova N D，Jackson J D，et al．The COG database:an updated version includes eukaryotes［J］．BMC Bioinformatics，2003，4(1):41．

［19］ Ogata H，Goto S，Fujibuchi W，et al．Computation with the KEGG pathway database［J］．Biosystems，1998，47(1/2):119 －128．

［20］ Karp P D，Ouzounis C A，Moore-Kochlacs C，et al．Expansion of the BioCyc collection of pathway/genome databases to 160 genomes［J］．Nucleic Acids Research，2005，33(19):6083－6089．

［21］ 't Hoen P A C，Ariyurek Y，Thygesen H H，et al．Deep sequencing-based expression analysis shows major advances in robustness，resolution and inter-lab portability over five microarray platforms［J］．Nucleic Acids Research，2008，36(21):e141．

［22］ Shendure J．The beginning of the end for microarrays? ［J］．Nature Methods，2008，5(7):585 －587．

［23］ Mortazavi A，Williams B A，Mccue K，et al．Mapping and quantifying mammalian transcriptomes by RNA-seq ［J］．Nature Methods，2008，5(7):621－628．

［24］ Vera J C，Wheat C W，Fescemyer H W，et al．Rapid transcriptome characterization for a nonmodel organism using 454 pyrosequencing［J］．Molecular Ecology，2008，17(7):1636－1647．

［25］ Hale M C，McCormick C R，Jackson J R，et al．Next-generation pyrosequencing of gonad transcriptomes in the polyploid lake sturgeon(Acipenser fulvescens):the relative merits of normalization and rarefaction in gene discovery ［J］．BMC Genomics，2009，10(1):203．

［26］ Salem M，Rexroad C E，Wang Jiannan，et al．Characterization of the rainbow trout transcriptome using Sanger and 454-pyrosequencing approaches［J］．BMC Genomics，2010，11(1):564．

［27］ Rokyta D R，Wray K P，Lemmon A R，et al．A high-throughput venom-gland transcriptome for the eastern diamondback rattlesnake(Crotalus adamanteus)and evidence for pervasive positive selection across toxin classes［J］．Toxicon，2011，57(5):657 －671．

［28］ Huh J W，Kim Y H，Park S J，et al．Large-scale transcriptome sequencing and gene analyses in the crab-eating macaque(Macaca fascicularis)for biomedical research ［J］．BMC Genomics，2012，13(1):163．

［29］ Gao Xianggang，Han Jiabo，Lu Zhichuang，et al．Characterization of the spotted seal Phoca largha transcriptome using Illumina pairedend sequencing and development of SSR markers［J］．Comparative Biochemistry and Physiology Part D:Genomics and Proteomics，2012，7(3):277－284．

［30］ 刘红亮．基于高通量测序的辽宁绒山羊转录组研究与应用［D］．杨凌:西北农林科技大学，2013．

［31］ Ponting C P，Oliver P L，Reik W．Evolution and functions of long noncoding RNAs［J］．Cell，2009，136(4):629 －641．

［32］ Berezikov E，Thuemmler F，Van Laake L W，et al．Diversity of microRNAs in human and chimpanzee brain ［J］．Nature Genetics，2006，38(12):1375－1377．