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．