全轉(zhuǎn)錄組學(xué)在畜牧業(yè)中的應(yīng)用

石田培，張莉

?

全轉(zhuǎn)錄組學(xué)在畜牧業(yè)中的應(yīng)用

石田培，張莉

中國(guó)農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193

RNA作為一種大分子參與基因編碼、解碼、調(diào)控、表達(dá)等多種生物學(xué)過(guò)程。目前，對(duì)RNA的功能研究主要通過(guò)全轉(zhuǎn)錄組測(cè)序方法來(lái)完成。全轉(zhuǎn)錄組研究可以對(duì)基因結(jié)構(gòu)與功能進(jìn)行更深層次地分析和探究，揭示基因表達(dá)與生命現(xiàn)象之間的內(nèi)在聯(lián)系?，F(xiàn)階段，基于高通量測(cè)序技術(shù)的轉(zhuǎn)錄本結(jié)構(gòu)研究、基因表達(dá)水平研究及非編碼區(qū)域功能研究在模式動(dòng)物、豬、禽類中已大量開展，但在羊上卻鮮有報(bào)道。本文介紹了利用RNA-seq及Small RNA-seq 技術(shù)研究全轉(zhuǎn)錄組的一般流程及常用策略，綜述了全轉(zhuǎn)錄組學(xué)技術(shù)在畜牧業(yè)領(lǐng)域中的研究進(jìn)展。

基因表達(dá)；全轉(zhuǎn)錄組；RNA測(cè)序技術(shù)；高通量測(cè)序

在后基因組學(xué)研究中，轉(zhuǎn)錄組學(xué)是解讀基因組功能元件和揭示細(xì)胞及組織分子機(jī)制的基礎(chǔ)，在生物表型和基因表達(dá)研究中占據(jù)了重要地位[1]。轉(zhuǎn)錄組一般指狹義轉(zhuǎn)錄組，即能夠編碼蛋白的RNA， 而生物學(xué)中的轉(zhuǎn)錄本概念則為生物體特定細(xì)胞或組織在某一特定條件下的所有轉(zhuǎn)錄產(chǎn)物，涵蓋了編碼RNA (mRNA)和非編碼RNA (non-coding RNA, ncRNA)[2]。同相對(duì)穩(wěn)定的基因組比較，轉(zhuǎn)錄組會(huì)根據(jù)生理狀態(tài)、生長(zhǎng)階段、生長(zhǎng)環(huán)境的變化而變化[3]，具有高度動(dòng)態(tài)可變的特性。除此之外，轉(zhuǎn)錄組圖譜中含有豐富的生物學(xué)信息，包括基因表達(dá)豐度及差異、基因結(jié)構(gòu)、反義轉(zhuǎn)錄本、選擇性剪切、單核苷酸多態(tài)性和基因融合等。

生物體是一個(gè)錯(cuò)綜復(fù)雜的網(wǎng)絡(luò)，構(gòu)成復(fù)雜有機(jī)體的任何簡(jiǎn)單分子都不能孤立存在或行使功能，如轉(zhuǎn)錄組中的各類RNA。這些不同類型的RNA生物學(xué)功能各異，參與眾多信號(hào)通路的調(diào)控。單一的mRNA或ncRNA研究缺乏關(guān)聯(lián)性，對(duì)多種RNA 信息進(jìn)行整合分析，探索潛在的網(wǎng)絡(luò)調(diào)控機(jī)制則成為轉(zhuǎn)錄組學(xué)研究的趨勢(shì)。在傳統(tǒng)轉(zhuǎn)錄組學(xué)(以mRNA[4]為主)研究中，樣品總RNA的提取通常利用Oligo (dT)磁珠富集法，該方法能有效發(fā)掘富含poly (A)尾巴的RNA[5]，包括mRNA和一小部分長(zhǎng)鏈非編碼RNA (long non-coding RNA, lncRNA)，但過(guò)濾掉了其他poly (A) minusRNA組分，即無(wú)poly (A)尾RNAs，因此導(dǎo)致circRNA (circular RNA)和部分lncRNA的丟失，不能全面地反映轉(zhuǎn)錄本的真實(shí)情況。Yang等[6]首次提出了基于poly (A) minus富集的RNA-seq方法，可將sno-lncRNAs (small nucleolar RNAs)、ciRNAs (circular intronic RNA)和circRNAs等這類非poly (A)結(jié)尾的新型RNA分子富集并保留下來(lái)。至此，包括miRNA[7]、lncRNA[8,9]和circRNA[10,11]等非編碼RNA的全轉(zhuǎn)錄組測(cè)序研究應(yīng)運(yùn)而生。本文綜述了基于RNA-seq 和Small RNA-seq技術(shù)的全轉(zhuǎn)錄組研究概況，總結(jié)了全轉(zhuǎn)錄組測(cè)序的一般流程和常用策略，最后對(duì)畜牧業(yè)未來(lái)的發(fā)展趨勢(shì)進(jìn)行展望。

1 全轉(zhuǎn)錄組測(cè)序概述

全轉(zhuǎn)錄組是指在一定時(shí)空狀態(tài)下有機(jī)體特定細(xì)胞或組織能夠轉(zhuǎn)錄出的所有轉(zhuǎn)錄本[12]，其中蘊(yùn)含著重要的生物學(xué)調(diào)控規(guī)律。豐富的RNA組學(xué)信息可用于開展差異表達(dá)ncRNA的靶基因預(yù)測(cè)、差異表達(dá)ncRNA和mRNA的正負(fù)相關(guān)性分析、ceRNA (com-peting endogenous RNAs)對(duì)靶向mRNA的功能富集分析及關(guān)鍵基因挖掘等研究。全轉(zhuǎn)錄測(cè)序是通過(guò)高通量RNA-seq測(cè)序技術(shù)對(duì)RNA序列進(jìn)行測(cè)序，并通過(guò)一些定量平臺(tái)和技術(shù)反應(yīng)其表達(dá)量和表達(dá)差異，從而形成表達(dá)圖譜。全轉(zhuǎn)錄組學(xué)分析就是以表達(dá)譜為基礎(chǔ)，對(duì)RNA進(jìn)行鑒定和注釋，預(yù)測(cè)相應(yīng)靶細(xì)胞或編碼潛能，并基于GO(Gene Ontology)、COG (Clusters of Orthologous Groups)、KEGG (Kyoto Encyclopedia of Genes and Genomes)等數(shù)據(jù)庫(kù)進(jìn)行功能富集、聚類分析、信息挖掘和通路探究等。不同于傳統(tǒng)轉(zhuǎn)錄組學(xué)，全轉(zhuǎn)錄組學(xué)在mRNA研究基礎(chǔ)上，涵蓋了lncRNA、circRNA、miRNA、siRNA(small interfering RNA)和piRNA(piwi-interacting RNA)等多種非編碼RNA及其之間的調(diào)控網(wǎng)絡(luò)分析。目前，轉(zhuǎn)錄組測(cè)序主要集中在單一RNA的定量表達(dá)和注釋分析上，而全轉(zhuǎn)錄組學(xué)測(cè)序可同時(shí)對(duì)多種RNA進(jìn)行鑒定及關(guān)聯(lián)分析(圖1)。

2 全轉(zhuǎn)錄組測(cè)序研究的一般技術(shù)流程和方法

2.1 測(cè)序樣品的準(zhǔn)備

在樣品采集過(guò)程中，組織中的RNA 極易受內(nèi)源或外源RNA酶作用而降解，同樣也容易受到蛋白質(zhì)、DNA、同源和外源酚類等物質(zhì)的污染，因此，樣品質(zhì)量和保存條件是決定試驗(yàn)結(jié)果質(zhì)量的關(guān)鍵。不同物種不同組織部位的RNA得率大不相同，在動(dòng)物中，肝臟、脾臟和胰臟組織的得率較大，而肌肉、腦和心臟等組織得率較小。在皮膚、骨骼和毛發(fā)中RNA的提取難度較大，效果也較差。在準(zhǔn)備樣品時(shí)，首選新鮮組織，剔除非研究所需的組織類型。若組織體積較大，多管分組置于-80℃或液氮中長(zhǎng)期保存。

2.2 建庫(kù)測(cè)序

2.2.1 cDNA文庫(kù)構(gòu)建和RNA測(cè)序

總RNA質(zhì)檢后便可建庫(kù)，由于miRNA與其他3類RNA序列長(zhǎng)度差異較大，需要使用SE (single- end)和PE (paired-end)兩種測(cè)序策略同時(shí)建立兩個(gè)文庫(kù)。鏈特異性建庫(kù)的具體步驟如下(圖2)：(1)通過(guò)試劑盒在樣品總RNA中去除rRNA；(2)利用六堿基隨機(jī)引物合成第一條cDNA鏈；(3)第二條cDNA鏈合成時(shí)，用dUTP代替dTTP，使鏈上布滿U位點(diǎn)；(4)在3￠端加A，加接頭；(5)使用USER酶在尿嘧啶位置上產(chǎn)生一個(gè)單核苷酸缺口，借此消化掉第二條鏈，只保留第一條鏈；(6) PCR擴(kuò)增。小RNA建庫(kù)的具體步驟如下(圖3)：(1)在3￠端加接頭；(2)在5￠端加接頭；(3)反轉(zhuǎn)錄擴(kuò)增；(4)用凝膠純化富集擴(kuò)增的cDNA后測(cè)序。在測(cè)序時(shí)，理論上數(shù)據(jù)量越大越利于后續(xù)低豐度基因的完整組裝，但實(shí)際上并非數(shù)據(jù)量越大越好，需要根據(jù)物種情況及相關(guān)研究決定數(shù)據(jù)量的大小。

圖1 全轉(zhuǎn)錄組測(cè)序一般流程

圖2 鏈特異性建庫(kù)流程

2.2.2 原始數(shù)據(jù)預(yù)處理

測(cè)序產(chǎn)生的數(shù)據(jù)是一系列不能直接使用的原始數(shù)據(jù)，主要是FASTQ格式的讀段(reads)。其中除實(shí)驗(yàn)所需的堿基質(zhì)量信息之外，還包括測(cè)序儀器名稱、上機(jī)次數(shù)、試劑型號(hào)等信息。原始數(shù)據(jù)在建庫(kù)過(guò)程或者測(cè)序過(guò)程中會(huì)產(chǎn)生大量低質(zhì)量讀段、較短的序列、含N序列甚至是一些污染序列(如細(xì)菌基因序列)[13]、接頭序列、載體序列等，因此在數(shù)據(jù)分析之前，應(yīng)對(duì)原始數(shù)據(jù)進(jìn)行過(guò)濾、剪切和校正，以確保后期讀段定位、轉(zhuǎn)錄本組裝、基因定量等流程的順利開展。目前可用的質(zhì)控軟件是FastQC和NGSQC[14]、Qualimap2[15]、HTQC[16]、QCchain[17]、almostSigni-fican[18]、fastq-clean[19]及FaQCs[20]等，最常用的是FastQC和HTQC。FastQC可作為一個(gè)單獨(dú)的JAVA程序，在速度上遠(yuǎn)超其他工具，極短時(shí)間內(nèi)就可以運(yùn)行數(shù)千萬(wàn)的讀段，輸入文件的選擇性也較大，如壓縮或未壓縮的EASTQ或SAM/BAM文件。FastQC除能列出讀段的數(shù)目及質(zhì)量編碼以外，還能可視化報(bào)告堿基內(nèi)容及質(zhì)量、讀段長(zhǎng)度和重復(fù)序列等信息。在實(shí)際操作中，可聯(lián)合使用上述質(zhì)控軟件以達(dá)最佳過(guò)濾效果。

圖3 小RNA建庫(kù)流程

2.2.3 讀段定位及轉(zhuǎn)錄本重組裝

由原始數(shù)據(jù)進(jìn)行深度清理及質(zhì)量控制后獲得待分析數(shù)據(jù)(clean reads)，需要通過(guò)一系列軟件將讀段比對(duì)到參考基因組或者轉(zhuǎn)錄本上，并根據(jù)實(shí)際定位情況進(jìn)行轉(zhuǎn)錄本組裝[21]。在全轉(zhuǎn)錄組數(shù)據(jù)分析時(shí)，選擇比對(duì)程序時(shí)還需考慮剪接比對(duì)(unspliced aligners)情況[22]。當(dāng)生物體沒(méi)有內(nèi)含子或進(jìn)行miRNA測(cè)序時(shí)，可以使用Quality (MAQ)[23]、Burrow-Wheeler Aligner (BWA)[24]和Bowtie[25]等連續(xù)比對(duì)軟件，這些方法用于識(shí)別已知外顯子和接頭，不能識(shí)別由可變剪切所產(chǎn)生的新型外顯子。但當(dāng)讀段定位至外顯子-外顯子的交界接頭處或作圖到有內(nèi)含子的基因組時(shí)，則必須使用剪接比對(duì)程序，如TopHat[26]、Map-Splice[27]、STAR[28]和GSNAP[29]等。可先將讀段分成兩個(gè)短片段再參與定位，同時(shí)還可記錄分割數(shù)據(jù)以供查找后續(xù)轉(zhuǎn)錄本，該方法能夠識(shí)別由可變剪切產(chǎn)生的新轉(zhuǎn)錄本。一般情況下，先利用Bowtie進(jìn)行初步匹配，用匹配成功的reads來(lái)獲取基因組覆蓋區(qū)域，預(yù)測(cè)剪接點(diǎn)(splice junction)局域。然后利用Tophat將未成功匹配的reads劃分到splice junction序列上。如果研究物種沒(méi)有參考基因組或轉(zhuǎn)錄組時(shí)，需要將讀段自行組裝成可用的參考序列，然后將所有測(cè)序讀段通過(guò)從頭組裝生成重疊群和單一序列，再進(jìn)行后續(xù)數(shù)據(jù)分析。目前，有多種組裝程序可供分析使用[30,31]，如Newbler[32]、MIRA[33]、CAP3[34]、Seq-Man[35]、TGICL[36]、stackPACK[37]、Velvet[38]、AbySS[39]、ALLPATHS2[40]、Oases[41]、SOAPdenovo[42]、Multiple-k method[43]、Scaffolding using translation mapping (STM)、Trinity[44]和PCAP[45]等。

2.2.4 表達(dá)量化及差異分析

經(jīng)過(guò)讀段定位后，可根據(jù)讀段在轉(zhuǎn)錄本上的分布情況預(yù)測(cè)基因豐度。一般來(lái)說(shuō)，通過(guò)軟件分析獲得的注釋文件中會(huì)含有轉(zhuǎn)錄本分布信息，能夠通過(guò)分析讀段的匹配情況來(lái)識(shí)別新的轉(zhuǎn)錄本。在高等動(dòng)物的生命活動(dòng)調(diào)控過(guò)程中，存在可變剪切(alterna-tively splice, AS)等一系列復(fù)雜的調(diào)控系統(tǒng)，致使全部轉(zhuǎn)錄物不能直接對(duì)應(yīng)基因。此外，還存在不同轉(zhuǎn)錄本外顯子共享、讀段不能跨越多個(gè)剪接位點(diǎn)等問(wèn)題。因此，需要先進(jìn)行可變剪接識(shí)別[46]?，F(xiàn)階段，常用轉(zhuǎn)錄本識(shí)別軟件包括Stringtie[47]、cufflinks[48]、CIDANE[49]、GRIT[50]、TransComb[51]、iReckon[52]、SLIDE[53]、Montebello[54]、Augustus[55]、IsoLassocia-tion[56]、Scripture[57]、Traph[58]和MITIE等。人()和斑馬魚()等基因信息注釋完整的物種可直接進(jìn)行基因預(yù)測(cè)定量分析，但對(duì)于一些基因注釋信息不完整的物種需要先進(jìn)行轉(zhuǎn)錄本的預(yù)測(cè)。目前已完成的大量轉(zhuǎn)錄組測(cè)序數(shù)據(jù)快速完善了遺傳數(shù)據(jù)庫(kù)，反復(fù)檢驗(yàn)校對(duì)了基因注釋信息，推動(dòng)了分子生物學(xué)的研究進(jìn)展。

測(cè)序過(guò)程中測(cè)序深度、基因片段大小、運(yùn)用算法、實(shí)驗(yàn)批次等因素極易造成誤差，所以在定量時(shí)應(yīng)使用標(biāo)準(zhǔn)化的方法消除差異，最常用的樣本內(nèi)標(biāo)準(zhǔn)化方法包括RPKM[59]、FPKM、TMP[60]和KPKM等。根據(jù)比對(duì)到基因上的reads或者fragments數(shù)目，進(jìn)行基因長(zhǎng)度和測(cè)序總量歸一化后，即可統(tǒng)計(jì)表達(dá)量。常用的轉(zhuǎn)錄本定量軟件包括HTSeq[61]、feat-ureCounts[62]、StringTie、Cufflinks、RSEM[63]、Sailfish[64]、kallisto[65]、IsoLasso和NURD[66]等。為確定在不同條件或平臺(tái)的表達(dá)差異情況，經(jīng)過(guò)量化和標(biāo)準(zhǔn)化表達(dá)水平的基因仍需利用測(cè)序覆蓋度估測(cè)表達(dá)豐度分布[67,68]。

RNA-seq數(shù)據(jù)的差異表達(dá)(differential expression, DE)分析往往以單變量的方式進(jìn)行，很難擬合出一個(gè)普遍適用的統(tǒng)計(jì)模型，因此根據(jù)生物統(tǒng)計(jì)學(xué)原理將DE分析軟件進(jìn)行下列分析對(duì)比(表1)。為使表達(dá)分布符合統(tǒng)計(jì)方法的假設(shè)，RNA-seq數(shù)據(jù)會(huì)被某種方式歸一化。按照歸一化在DE分析中的前后順序可將常用軟件分為兩類：第一類是先進(jìn)行標(biāo)準(zhǔn)化處理再通過(guò)統(tǒng)計(jì)學(xué)原理計(jì)算表達(dá)差異，包括TMM[69]、DESeq[70]、PoissonSeq[71]和UpperQuartile[72]等；第二類是將歸一化作為DE分析的步驟，在處理時(shí)直接執(zhí)行歸一化，如FPKM、RPKM和TPM，但可信度較差。

表1 差異表達(dá)分析軟件

2.3 基因功能分析

生物信息學(xué)主要用數(shù)學(xué)及統(tǒng)計(jì)學(xué)的方法對(duì)生物信息進(jìn)行存儲(chǔ)、分配、檢索及分析，是一門集合生命科學(xué)、自然科學(xué)與信息科學(xué)的新學(xué)科。目前，生物信息學(xué)在生命科學(xué)研究中應(yīng)用十分廣泛，在全轉(zhuǎn)錄組研究中也占據(jù)重要地位。在全轉(zhuǎn)錄組學(xué)分析中主要以差異基因的功能富集、候選基因的功能預(yù)測(cè)、調(diào)控網(wǎng)絡(luò)的構(gòu)建等分析為主。當(dāng)前以GO、COG及KEGG數(shù)據(jù)庫(kù)為基礎(chǔ)的分析工具包、軟件和網(wǎng)站眾多，如基于R語(yǔ)言的GSA、PADOG、SAFE、Globaltest、Sigpathway、GAGE、GSVA、PLAGE、ZSCORE、SSGSEA、MRGSE、ANCOVA、CAMERA、SPIA、TopoGSA、ToPASeq、NetGSA、NEA和GOGANPA等工具包；DAVID、GOstat、GenMAPP、GOMiner、Onto-Express、EnrichNet和NOA等網(wǎng)站；MetaCore、Cytoscape和GSEA等軟件。以上分析軟件各有千秋，功能也不盡相同，在數(shù)據(jù)處理與分析時(shí)，應(yīng)根據(jù)實(shí)驗(yàn)?zāi)窟x擇合適的方法。如需更系統(tǒng)地反映分子調(diào)控機(jī)制，還可結(jié)合其他基因組學(xué)手段如全基因組測(cè)序、RNA甲基化、表觀修飾等數(shù)據(jù)進(jìn)行聯(lián)合分析。

3 全轉(zhuǎn)錄組學(xué)技術(shù)在畜牧業(yè)中的應(yīng)用

3.1 全轉(zhuǎn)錄組測(cè)序在豬中的應(yīng)用

豬()的飼養(yǎng)對(duì)畜牧業(yè)發(fā)展乃至民生都至關(guān)重要。不同品種豬的肌纖維類型對(duì)豬肉品質(zhì)有很大影響，但潛在的分子機(jī)制仍不清楚。近年來(lái)，有關(guān)豬非編碼RNA的功能研究開展得很多。為探究miRNA在豬骨骼肌中的作用，Mcdaneld等[75]分別選取了增殖中的衛(wèi)星細(xì)胞(4~代)、胚胎(60、90、105日齡)、出生胎兒和成年豬組織進(jìn)行了全轉(zhuǎn)錄組研究，發(fā)現(xiàn)了12個(gè)新型miRNA與肌肉生長(zhǎng)發(fā)育密切相關(guān)，其中肌肉特異性miR-206在衛(wèi)星細(xì)胞中近乎缺失，但在其他發(fā)育階段表達(dá)量卻很高；miR-1在成年豬中的表達(dá)豐度最高；miR-133在胚胎期和初生時(shí)的豐度極低甚至檢測(cè)不到，但miR-368、miR-376和miR-423-5p在初生豬中極高；miR-432表達(dá)量在胚胎發(fā)育早期最高，隨后逐漸降低。該研究是對(duì)豬胚胎發(fā)育期骨骼肌miRNA較為全面的轉(zhuǎn)錄組分析，為深入探究豬骨骼肌miRNA調(diào)控機(jī)制提供了理論基礎(chǔ)。沈一飛[76]利用RNA-seq和Small RNA-seq技術(shù)對(duì)瘦肉型約克夏豬和脂肪型金華豬進(jìn)行了甲狀腺組織的mRNA、lncRNA和miRNA共表達(dá)鑒定與比較，結(jié)果發(fā)現(xiàn)差異表達(dá)mRNA 492個(gè)，差異表達(dá)lncRNA 48個(gè)和miRNA 18個(gè)。通過(guò)功能分析和聚類，共有256個(gè)RNA(其中18個(gè)miRNA、1個(gè)lncRNA)參與到同一個(gè)調(diào)控網(wǎng)絡(luò)中。在整個(gè)調(diào)控網(wǎng)絡(luò)中，發(fā)現(xiàn)ssc- miRNA-221-5p、ssc-miRNA-708-5p、ssc-miRNA-532- 3p和novel_12等發(fā)揮重要調(diào)控作用。Li等[77]對(duì)從胎兒期到成年期期間獲得的組織混合物制備的10個(gè)小RNA測(cè)序文庫(kù)中的榮昌豬(Rongchang pigs) miRNA進(jìn)行全面檢測(cè)，通過(guò)哺乳動(dòng)物miRNA、前體發(fā)夾(pre-miRNA)、高覆蓋率豬基因組裝配(2009年4月)和表達(dá)序列標(biāo)簽(EST)的分析，將豬miRNAome的所有組成部分?jǐn)U展到867個(gè)pre- miRNAs (623個(gè)基因組坐標(biāo))，編碼1004個(gè)miRNA，其中777個(gè)是獨(dú)特的。對(duì)47個(gè)組織特異性樣品中選定的30個(gè)miRNA進(jìn)行qRT-PCR實(shí)驗(yàn)，發(fā)現(xiàn)測(cè)序數(shù)據(jù)和試驗(yàn)結(jié)果一致。Sun等[78]對(duì)長(zhǎng)白豬(Landrace)和蘭塘豬(Lantang pigs)背最長(zhǎng)肌進(jìn)行了全轉(zhuǎn)錄組測(cè)序研究，從22 469個(gè)編碼轉(zhuǎn)錄物中篩選出547個(gè)差異表達(dá)mRNA，通過(guò)生物信息學(xué)分析挖掘出與肌肉生長(zhǎng)發(fā)育相關(guān)的 17個(gè)基因。此外還發(fā)現(xiàn)差異表達(dá)的5566個(gè)lncRNA和4360個(gè)circRNA。其中，3376個(gè)lncRNA和1401個(gè)circRNA在Lantang文庫(kù)中上調(diào)表達(dá)，而1590個(gè)lncRNA和2959個(gè)circRNA下調(diào)表達(dá)。通過(guò)結(jié)合匹配的miRNA譜分析測(cè)序數(shù)據(jù)，鑒定出26種參與ceRNA網(wǎng)絡(luò)的海綿載體，包括19種lncRNA、40種circRNA和9種mRNA。全轉(zhuǎn)錄組研究提供了一種全新的分析方法，對(duì)解析豬肌肉生長(zhǎng)發(fā)育規(guī)律和疾病發(fā)生機(jī)制具有重要意義。

3.2 全轉(zhuǎn)錄組測(cè)序在禽類中的應(yīng)用

我國(guó)家禽遺傳資源豐富，地方品種各具特色，但由于總體生產(chǎn)水平較低，許多優(yōu)良性狀利用效率不高。隨著分子生物學(xué)和生物信息學(xué)的發(fā)展與融合，通過(guò)各種技術(shù)和手段，已鑒定出一批與生長(zhǎng)、繁殖等重要經(jīng)濟(jì)性狀相關(guān)的分子標(biāo)記和候選基因。Yu等[79,80]選取了6和10日齡的雞胚左(L)、右(R)卵巢作為樣品，通過(guò)全轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)卵巢的退化進(jìn)行了研究。在6R樣品中產(chǎn)生31 066 414個(gè)序列讀數(shù)，在6L樣品中產(chǎn)生31 900 200個(gè)序列讀數(shù)，在10R樣品中產(chǎn)生31 400 070個(gè)序列讀數(shù)，在10L樣品中讀取35 504746個(gè)序列。這些數(shù)據(jù)中，大約73.33%的序列可以定位到雞胚卵巢的參考基因組序列上。雞胚卵巢的發(fā)育受許多基因和信號(hào)通路的調(diào)控，通過(guò)對(duì)差異表達(dá)基因比對(duì)和通路功能聚類發(fā)現(xiàn)了22個(gè)與卵巢發(fā)育和退化相關(guān)的基因。其中轉(zhuǎn)錄水平排名前20的卵巢基因可能與碳代謝、吞噬體及鈣信號(hào)傳導(dǎo)密切相關(guān)。Glazov等[81]在已發(fā)現(xiàn)的miRNA基礎(chǔ)上進(jìn)行了深度挖掘，將3個(gè)小RNA文庫(kù)同時(shí)比較并嚴(yán)格區(qū)分真正miRNA前體與結(jié)構(gòu)相似的RNA，共鑒定到361個(gè)新的miRNA、88個(gè)新的miRNA候選物、18個(gè)mirtrons (包括6個(gè)新的非典型mirtron候選物)和21個(gè)mirtron候選物。為識(shí)別潛在的禽類特異性miRNA，同時(shí)與人、狗()、斑馬魚、爪蟾()和蜥蜴()進(jìn)行保守進(jìn)化分析，結(jié)果表明只有6種miRNA在非禽類脊椎動(dòng)物中具有保守性，剩余miRNA可能具有鳥類和/或雞系特異性。Li等[80]分別從孵育10天、12天、14天和18天的白來(lái)航雞(White Leghorns)蛋胚中收集骨骼肌(胸大肌)，經(jīng)過(guò)RNA-seq技術(shù)測(cè)序并與已發(fā)布的數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，篩選到281個(gè)新型基因間lncRNA，對(duì)這些lncRNA進(jìn)行保守性分析，并利用UCSC數(shù)據(jù)庫(kù)評(píng)分，結(jié)果表明以上lncRNA的保守性均高于隨機(jī)非編碼序列，但遠(yuǎn)低于蛋白質(zhì)編碼基因。該研究是首例有關(guān)白來(lái)航雞骨骼肌lncRNA的分析，鑒定出的新型lncRNA極大豐富了雞ncRNA數(shù)據(jù)庫(kù)。

近年來(lái)，鵝()以其適應(yīng)性強(qiáng)、生長(zhǎng)快、營(yíng)養(yǎng)成分豐富和投入要求低而備受關(guān)注。此前，Kang等[82]通過(guò)抑制性消減雜交(supp-ression subtractive hybridization, SSH)方法鑒定了一些與鵝從產(chǎn)卵階段到產(chǎn)蛋階段繁殖相關(guān)的差異表達(dá)基因，Guo等[83]也使用相同的方法在產(chǎn)蛋階段和育雛階段發(fā)現(xiàn)了若干差異表達(dá)的基因。Xu等[84]采用短讀序列技術(shù)(Illumina)對(duì)10只380日齡的雌性浙東白鵝(Goose)的卵巢樣本進(jìn)行了從頭轉(zhuǎn)錄組裝配，使用Illumina RNA-seq和DGE深度測(cè)序并繪制出鵝卵巢組織的轉(zhuǎn)錄組圖譜，得到67 315 996個(gè)100 bp的短讀數(shù)，組裝成130 514個(gè)獨(dú)特序列。基于已知蛋白質(zhì)的BLAST結(jié)果，分析鑒定到52 642個(gè)目標(biāo)序列。該研究分析了鵝產(chǎn)蛋、育雛期間的轉(zhuǎn)錄變化情況，鑒定出大量與卵泡發(fā)育和生殖相關(guān)的基因。

病理性肥胖是鴨養(yǎng)殖業(yè)所面臨的重要問(wèn)題之一，其分子機(jī)制仍然未知。Chen等[85]對(duì)家養(yǎng)鴨與野鴨兩個(gè)品種腹部脂肪進(jìn)行轉(zhuǎn)錄組測(cè)序分析，預(yù)測(cè)了23 699個(gè)未注釋基因，確定了753個(gè)差異表達(dá)基因。在北京鴨()中，一些與脂質(zhì)代謝的相關(guān)基因(、和等)和致癌基因(、和等)上調(diào)表達(dá)，而與腫瘤抑制和免疫相關(guān)的基因(、、和等被下調(diào)，這些數(shù)據(jù)表明鴨的腫瘤發(fā)生可能與病理性肥胖密切相關(guān)。此外，發(fā)現(xiàn)280 576個(gè)單核苷酸變異在兩個(gè)品種之間存在差異，包括8641個(gè)異構(gòu)變異，富含參與脂質(zhì)和免疫相關(guān)途徑的基因，表明與鴨的代謝功能和免疫相關(guān)功能密切相關(guān)。

3.3 全轉(zhuǎn)錄組測(cè)序在反芻動(dòng)物中的應(yīng)用

近年來(lái)國(guó)內(nèi)牛羊肉市場(chǎng)需求不斷增加，這對(duì)牛羊的育種和養(yǎng)殖工作提出了更高的要求，只有充分了解牛羊生長(zhǎng)發(fā)育與繁殖等性狀的發(fā)生機(jī)制，才能提高生產(chǎn)效率，全轉(zhuǎn)錄組學(xué)則為其提供了全新的研究手段。Di等[86]在灘羊(Tan)和小尾寒羊(STH)的卵巢中鑒定出483個(gè)miRNA (包括97個(gè)已知的、369個(gè)保守的和17個(gè)預(yù)測(cè)的新miRNA)?；贙EGG分析，一些差異表達(dá)miRNA的靶基因參與生殖激素相關(guān)途徑(如類固醇生物合成、雄激素和雌激素代謝和GnRH信號(hào)傳導(dǎo)途徑)以及卵泡和黃體發(fā)育相關(guān)途徑，這對(duì)綿羊的繁育具有重要意義。Chang等[87]使用新一代測(cè)序技術(shù)(Solexa高通量測(cè)序技術(shù))研究了綿羊黃體期卵巢組織，鑒定出267種新的miRNA，并利用qRT-PCR和Northern印跡證實(shí)了在綿羊卵巢和睪丸中表達(dá)的一種新型miRNA (ovis_aries_ovary m0033_3p)。根據(jù)序列和結(jié)構(gòu)的一致性，推測(cè)ovis_ aries_ovarym0033_3p具有類似于hsa-miR-214-3p的功能，能夠參與細(xì)胞存活、胚胎發(fā)育、繁育生殖和卵巢癌抗性的精細(xì)調(diào)節(jié)。張世芳等[88]采用Solexa技術(shù)對(duì)5頭特克賽爾羊(Texel)進(jìn)行miRNA深度測(cè)序，獲得了16 532 850條原始序列讀數(shù)。通過(guò)與哺乳動(dòng)物成熟miRNA數(shù)據(jù)庫(kù)、miRNA前體序列、綿羊基因組數(shù)據(jù)庫(kù)的比對(duì)分析，更新miRNA前體序列庫(kù)至1529條，編碼的miRNA成熟體序列增至1999條。Miao等[89]對(duì)道賽特綿羊(Dorset)和小尾寒羊(Han)卵巢組織測(cè)序，鑒定出可能參與繁殖力調(diào)節(jié)的候選基因，這些候選基因參與各種細(xì)胞活動(dòng)，如代謝級(jí)聯(lián)、催化功能和信號(hào)轉(zhuǎn)導(dǎo)。此外，通過(guò)miRNA譜分析鑒定了每組綿羊特有的特異性miRNA，發(fā)現(xiàn)若干與生殖力調(diào)控相關(guān)的miRNA。

Billerey等[90]檢測(cè)了9頭利木贊牛犢(Limousin)的胸肌樣本，每個(gè)文庫(kù)約獲得14~45百萬(wàn)個(gè)配對(duì)末端讀數(shù)，發(fā)現(xiàn)418種lincRNA (large intergenic non-coding RNAs)，與已知的10 775種蛋白編碼基因存在顯著差異。Sun等[91]利用Ribo-Zero RNA-seq技術(shù)深度剖析了胚胎、犢牛和成年牛骨骼肌的轉(zhuǎn)錄組譜，發(fā)現(xiàn)犢牛和成年牛之間的表達(dá)水平高度相關(guān)。在胚胎期有數(shù)百個(gè)基因顯著表達(dá)，但在出生后至少減少了10倍，表明這些基因在肌肉發(fā)育中具有潛在作用。此外，該研究首次分析了牛骨骼肌中全部轉(zhuǎn)錄異構(gòu)體，發(fā)掘出36 694個(gè)新型異構(gòu)體，檢測(cè)到185 036個(gè)SNP和12 428個(gè)短插入缺失(InDel)位點(diǎn)。研究發(fā)現(xiàn)可變剪接事件、SNP和InDel的數(shù)量在胚胎中比在犢牛和成年牛中更多，這表明基因表達(dá)在胚胎中最活躍。Cánovas等[92]通過(guò)轉(zhuǎn)錄組學(xué)測(cè)序?qū)?頭荷斯坦奶牛(Holstein cow)的乳樣品進(jìn)行了SNP篩選，共檢測(cè)到19 175個(gè)差異表達(dá)基因，100 734個(gè)SNP位點(diǎn)，其中33 045個(gè)與荷斯坦奶牛SNP位點(diǎn)重合，這些SNP位于泌乳期間表達(dá)基因的編碼區(qū)中，可用于荷斯坦奶牛乳用性狀的基因變異分析和關(guān)聯(lián)研究。

在梅花鹿()的育種研究中，Yao等[93]對(duì)鹿茸進(jìn)行了轉(zhuǎn)錄組測(cè)序，組裝出89 001個(gè)獨(dú)特序列(平均大小450 bp)，發(fā)現(xiàn)了一些高度表達(dá)的基因參與鹿茸快速生長(zhǎng)的調(diào)節(jié)，包括轉(zhuǎn)錄因子、信號(hào)分子和細(xì)胞外基質(zhì)蛋白。這些數(shù)據(jù)是當(dāng)前梅花鹿最全面的序列資源，為鹿的分子遺傳學(xué)和功能基因組學(xué)的研究提供了理論基礎(chǔ)。

3.4 全轉(zhuǎn)錄組測(cè)序在馬屬動(dòng)物中的應(yīng)用

全轉(zhuǎn)錄組研究在馬屬動(dòng)物中開展得較少，目前主要是通過(guò)轉(zhuǎn)錄組學(xué)測(cè)序技術(shù)進(jìn)行基因挖掘、注釋和功能預(yù)測(cè)。Xie等[94]從頭組裝了驢()白細(xì)胞的轉(zhuǎn)錄組，鑒定出264 714個(gè)不同序列，預(yù)測(cè)了38 949個(gè)蛋白質(zhì)片段。通過(guò)比較驢、馬()和野馬()的蛋白質(zhì)序列，將驢蛋白片段與哺乳動(dòng)物表型相關(guān)聯(lián)。通過(guò)比較驢和馬的外耳性狀相關(guān)蛋白，鑒定出3種與耳形大小相關(guān)的蛋白質(zhì)HIC1、PRKRA和KMT2A。Scott等[95]通過(guò)轉(zhuǎn)錄組測(cè)序?qū)︸R的lncRNA進(jìn)行注釋，發(fā)現(xiàn)了20 800新型轉(zhuǎn)錄本，證明了lncRNA獨(dú)有的特征，包括低表達(dá)、低外顯子多樣性和低水平的序列保守性。該研究結(jié)果所提供的候選基因可作為日后lncRNA注釋的基準(zhǔn)。

近年來(lái)發(fā)表了大量關(guān)于馬運(yùn)動(dòng)機(jī)能、骨骼發(fā)育的文獻(xiàn)報(bào)道，一部分是通過(guò)馬組織的RNA-seq數(shù)據(jù)改進(jìn)蛋白質(zhì)編碼基因的結(jié)構(gòu)注釋，另一部分是對(duì)RNA序列的分析。如Park等[96]對(duì)6匹純種馬運(yùn)動(dòng)前后的血液和肌肉進(jìn)行全轉(zhuǎn)錄測(cè)序，通過(guò)與前人的研究對(duì)比，發(fā)現(xiàn)超過(guò)19 417個(gè)新型單基因簇，鑒定出189 973個(gè)單核苷酸位點(diǎn)變異(single nucleotide variants, SNV)。使用差異表達(dá)分析，確定了多個(gè)運(yùn)動(dòng)調(diào)節(jié)基因：血液中有62個(gè)上調(diào)基因和80個(gè)下調(diào)基因，肌肉中有878個(gè)上調(diào)基因和285個(gè)下調(diào)基因。結(jié)果表明，在差異表達(dá)的基因中有91個(gè)轉(zhuǎn)錄因子編碼基因，其中包括56個(gè)功能未知的轉(zhuǎn)錄因子候選物可能與早期調(diào)節(jié)運(yùn)動(dòng)機(jī)制相關(guān)；此外，還發(fā)現(xiàn)了新型RNA表達(dá)模式：同一基因的不同可變剪接形式在運(yùn)動(dòng)前后表現(xiàn)出反向表達(dá)模式。該研究首次提供了馬轉(zhuǎn)錄組數(shù)據(jù)和較為全面的分析結(jié)果，包括運(yùn)動(dòng)前后表達(dá)的基因，以及與運(yùn)動(dòng)相關(guān)的候選基因：6個(gè)運(yùn)動(dòng)相關(guān)基因和91個(gè)早期調(diào)節(jié)轉(zhuǎn)錄因子，3個(gè)高SNV密度的基因，以及4個(gè)交替表達(dá)的剪接體。

3.5 全轉(zhuǎn)錄組測(cè)序在其他特種經(jīng)濟(jì)動(dòng)物中的應(yīng)用

特種經(jīng)濟(jì)動(dòng)物養(yǎng)殖已成為調(diào)整農(nóng)村產(chǎn)業(yè)結(jié)構(gòu)、發(fā)展特色經(jīng)濟(jì)的新亮點(diǎn)，為了發(fā)揮皮毛的最大經(jīng)濟(jì)價(jià)值，科研人員開始研究關(guān)于被毛顏色的調(diào)控機(jī)制。牛曉艷等[97]首先對(duì)不同毛色的獺兔(Rex rabbit)進(jìn)行測(cè)序和差異基因分析，找到12 408個(gè)差異表達(dá)基因，然后通過(guò)KEGG分析將得到的差異基因聚類富集到相關(guān)代謝通路上，結(jié)果發(fā)現(xiàn)8個(gè)與黑素細(xì)胞分化相關(guān)的差異基因。宋興超[98]利用RNA-seq技術(shù)對(duì)水貂()被毛色素沉積機(jī)理進(jìn)行了研究，鑒定出不同時(shí)期被毛黑素含量的變化，并根據(jù)水貂皮膚組織中成熟黑素細(xì)胞的分布特點(diǎn)，開展SNPs檢測(cè)，將 不同基因突變體與水貂毛色表型進(jìn)行關(guān)聯(lián)分析。mRNA定量表達(dá)驗(yàn)證結(jié)果表明，和等基因參與了黑素細(xì)胞發(fā)育、黑素小體前體形成、黑素小體轉(zhuǎn)運(yùn)和真黑和褐黑色素合成等生物學(xué)過(guò)程。

4 結(jié)語(yǔ)與展望

全轉(zhuǎn)錄組學(xué)以其精準(zhǔn)、系統(tǒng)、直觀的技術(shù)優(yōu)勢(shì)為畜禽重要經(jīng)濟(jì)性狀功能基因的挖掘、鑒定與驗(yàn)證提供了新的技術(shù)平臺(tái)和手段，并已廣泛地運(yùn)用在臨床醫(yī)學(xué)、藥學(xué)、生物學(xué)、水產(chǎn)學(xué)和農(nóng)林學(xué)等多個(gè)領(lǐng)域，為人類疾病研究、新藥研發(fā)和動(dòng)植物育種等開辟了新的思路。但是，全轉(zhuǎn)錄組學(xué)在畜牧領(lǐng)域的研究較其他領(lǐng)域而言起步較晚，研究還不夠深入，尤其在羊上，轉(zhuǎn)錄組的研究目前還主要集中在小RNA測(cè)序和基因注釋上。本研究團(tuán)隊(duì)將開展綿羊全轉(zhuǎn)錄組研究，并對(duì)其生長(zhǎng)性狀、肉用性狀等重要經(jīng)濟(jì)性狀進(jìn)行解析和應(yīng)用。

[1] Lockhart DJ, Winzeler EA. Genomics, gene expression and DNA arrays.,2000, 405(6788): 827–836.

[2] Lindberg J, Lundeberg J. The plasticity of the mammalian transcriptome.,2010, 95(1): 1–6.

[3] Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics.,2009, 10(1): 57–63.

[4] Wong ML, Medrano JF. Real-time PCR for mRNA quantitation.,2005, 39(1): 75–85.

[5] Liu C, LI X, Chen LL. Methods for genome-wide characterization of long noncoding RNAs.,2016, (6): 745–752.

[6] Yang L, Duff MO, Graveley BR, Carmichael GG, Chen LL. Genomewide characterization of non-polyadenylated RNAs.,2011, 12(2): R16.

[7] Carthew RW, Sontheimer EJ. Origins and mechanisms of miRNAs and siRNAs.,2009, 136(4): 642–655.

[8] Chen LL. Linking long noncoding RNA localization and function.,2016, 41(9): 761–772.

[9] Yin QF, Yang L, Zhang Y, Xiang JF, Wu YW, Carmichael G, Chen LL. Long noncoding RNAs with snoRNA ends.,2012, 48(2): 219–230.

[10] Zhang Y, Zhang X, Chen T, Xiang J, Yin Q, Xing Y, Zhu S, Yang L, Chen L. Circular intronic long noncoding RNAs.,2013, 51(6): 792–806.

[11] Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L. Complementary sequence-mediated exon circularization.,2014, 159(1): 134–147.

[12] Qi YX, Liu YB, Rong WH. RNA-Seq and its applications: a new technology for transcriptomics.,2011, 33(11): 1191–1201.祁云霞, 劉永斌, 榮威恒. 轉(zhuǎn)錄組研究新技術(shù): RNA- Seq及其應(yīng)用. 遺傳, 2011, 33(11): 1191–1201.

[13] Hong QY, Bi XJ, Wang DN, Li ZZ, Yu H, Xia NS, Li SW. Research progress on RNA-Seq technology.,2017, 37(6): 443–448.洪奇陽(yáng), 畢行建, 王大寧, 李子真, 俞海, 夏寧邵, 李少偉. 轉(zhuǎn)錄組測(cè)序技術(shù)研究進(jìn)展. 中國(guó)生化藥物雜志, 2017, 37(6): 443–448.

[14] Patel RK, Jain M. NGS QC toolkit: a toolkit for quality control of next generation sequencing data.,2012, 7(2): e30619.

[15] Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data.,2016, 32(2): 292–294.

[16] Yang X, Liu D, Liu F, Wu J, Zou J, Xiao X, Zhao F, Zhu B. HTQC: a fast quality control toolkit for Illumina sequencing data.,2013, 14: 33.

[17] Zhou Q, Su X, Wang A, Xu J, Ning K. QC-Chain: fast and holistic quality control method for next-generation sequencing data.,2013, 8(4): e60234.

[18] Ward J, Cole C, Febrer M, Barton GJ. Almost significant: simplifying quality control of high-throughput sequencing data.,2016, 32(24): 3850–3851.

[19] Zhang M, Sun H, Fei Z, Zhan F, Gong X, Gao S. Fastq_clean: an optimized pipeline to clean the Illumina sequencing data with quality control. In: IEEE International Conference on Bioinformatics and Biomedicine. 2015, 44–48.

[20] Lo CC, Chain PSG. Rapid evaluation and quality control of next generation sequencing data with FaQCs.,2014, 15(1): 366.

[21] Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq.,2011, 8(6): 469–477.

[22] Yang IS, Kim S. Analysis of whole transcriptome sequencing data: workflow and software.,2015, 13(4): 119–125.

[23] Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores.,2008, 18(11): 1851–1858.

[24] Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform., 2009, 25(4): 1754–1760.

[25] Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.,2009, 10(3): R25.

[26] Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq.,2009, 25(9): 1105–1111.

[27] Wang K, Singh D, Zeng Z, Coleman SJ, Huang Y, Savich GL, He X, Mieczkowski P, Grimm SA, Perou CM, MacLeod JN, Chiang DY, Prins JF, Liu J. MapSplice: accurate mapping of RNA-seq reads for splice junction discovery.,2010, 38(18): e178.

[28] Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner.,2013, 29(1): 15–21.

[29] Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads.,2010, 26(7): 873–881.

[30] Kumar S, Blaxter ML. Comparing de novo assemblers for 454 transcriptome data.,2010, 11: 571.

[31] Garg R, Patel RK, Tyagi AK, Jain M. De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification., 2011, 18(1): 53–63.

[32] Margulies M, Egholm M, Altman WE, Attiya S, Bader JS, Bemben LA, Berka J, Braverman MS, Chen YJ, Chen Z. Genome sequencing in microfabricated high-density picolitre reactors.,2005, 437(7057): 376–380.

[33] Burlibasa C, Vasiliu D, Vasiliu M. Genome sequence assembly using trace signals and additional sequence information. In: German Conference on Bioinformatics. 1999, 45–56.

[34] Huang X, Madan A. CAP3: a DNA sequence assembly program.,1999, 9(9): 868–877.

[35] Swindell SR, Plasterer TN. SEQMAN. Contig assembly.,1997, 70: 75–89.

[36] Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, Tsai J, Quackenbush J.gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets.,2003, 19(5): 651–652.

[37] Miller RT, Christoffels AG, Gopalakrishnan C, Burke J, Ptitsyn AA. A comprehensive approach to clustering of expressed human gene sequence: the sequence tag alignment and consensus knowledge base.,1999, 9(11): 1143–1155.

[38] Zerbino DR, Birney E. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs.,2008, 18(5): 821–829.

[39] Jackman SD, Birol ?. Assembling genomes using short- read sequencing technology.,2010, 11(1): 202–202.

[40] MacCallum I, Przybylski D, Gnerre S, Burton J, Shlyakhter I, Gnirke A, Malek J, McKernan K, Ranade S, Shea TP, Williams L, Young S, Nusbaum C, Jaffe DB. ALLPATHS 2: small genomes assembled accurately and with high continuity from short paired reads.,2009, 10(10): R103–R103.

[41] Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels.,2012, 28(8): 1086–1092.

[42] Li R, Li Y, Kristiansen K, Wang J, Wang J. SOAP: short oligonucleotide alignment program.,2008, 24(5): 713–714.

[43] Surget-Groba Y, Montoya-Burgos JI. Optimization of de novo transcriptome assembly from next-generation sequencing data.,2010, 20(10): 1432–1440.

[44] Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome.,2011, 29(7): 644–652.

[45] Huang X, Wang J, Aluru S, Yang SP, Hillier L. PCAP: a whole-genome assembly program.,2003, 13(9): 2164–2170.

[46] Huh JW, Kim YH, Park SJ, Kim DS, Lee SR, Kim KM, Jeong KJ, Kim JS, Song BS, Sim BW, Kim SU, Kim SH, Chang KT. Large-scale transcriptome sequencing and gene analyses in the crab-eating macaque () for biomedical research.,2012, 13: 163.

[47] Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads.,2015, 33(3): 290–295.

[48] Roberts A, Pimentel H, Trapnell C, Pachter L. Identification of novel transcripts in annotated genomes using RNA-Seq.,2011, 27(17): 2325–2329.

[49] Canzar S, Andreotti S, Weese D, Reinert K, Klau GW. CIDANE: comprehensive isoform discovery and abundance estimation.,2016, 17: 16.

[50] Boley N, Stoiber MH, Booth BW, Wan KH, Hoskins RA, Bickel PJ, Celniker SE, Brown JB. Genome-guided transcript assembly by integrative analysis of RNA sequence data.,2014, 32(4): 341–346.

[51] Liu JT, Yu T, Tao J, Li GJ. TransComb: genome-guided transcriptome assembly via combing junctions in splicing graphs.,2016, 17: 213.

[52] Mezlini AM, Smith EJ, Fiume M, Buske O, Savich GL, Shah S, Aparicio S, Chiang DY, Goldenberg A, Brudno M. iReckon: simultaneous isoform discovery and abundance estimation from RNA-seq data.,2013, 23(3): 519–529.

[53] Li JJ, Jiang CR, Brown JB, Huang H, Bickel PJ. Sparse linear modeling of next-generation mRNA sequencing (RNA-Seq) data for isoform discovery and abundance estimation.,2011, 108(50): 19867–19872.

[54] Hiller D, Wong WH. Simultaneous isoform discovery and quantification from RNA-seq.,2013, 5(1): 100–118.

[55] Stanke M, Keller O, Gunduz I, Hayes A, Waack S, Morgenstern B. AUGUSTUS: ab initio prediction of alternative transcripts.,2006, 34(web server issue): 435–439.

[56] Li W, Feng J, Jiang T. IsoLasso: A LASSO regression approach to RNA-seq based transcriptome assembly.,2011, 18(11): 1693–1707.

[57] Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Xian A, Fan L, Koziol MJ, Gnirke A, Nusbaum C. Ab initio reconstruction of cell type–specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs.,2010, 28(5): 503–510.

[58] Tomescu, Alexandru I, Kuosmanen, Anna, Makinen, Veli, Rizzi R. A novel min-cost flow method for estimating transcript expression with RNA-Seq.,2013, 14(Suppl.5): S15.

[59] Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq.,2008, 5(7): 621–628.

[60] Pachter L. Models for transcript quantification from RNA-Seq. 2013.

[61] Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data.,2015, 31(2): 166–169.

[62] Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.,2014, 30(7): 923–930.

[63] Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.,,2011, 12: 323.

[64] Rob P, Mount SM, Kingsford C. Sailfish enables alignment-free isoform quantification from RNA-seq reads using lightweight algorithms.,2013, 32(5): 462–464.

[65] Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification.,2016, 34(5): 525–527.

[66] Bullard JH, Purdom E , Hansen KD, Dudoit S. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments.,2010, 11: 94.

[67] Ma X, Zhang X. NURD: an implementation of a new method to estimate isoform expression from non-uniform RNA-seq data.,2013, 14: 220.

[68] Oshlack A, Wakefield MJ. Transcript length bias in RNA-seq data confounds systems biology.2009, 4: 14.

[69] Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data.,2010, 11(3): 1–9.

[70] Anders S, Huber W. Differential expression analysis for sequence count data.,2010, 11: R106.

[71] Li J, Witten DM, Johnstone IM, Tibshirani R. Normalization, testing, and false discovery rate estimation for RNA- sequencing data.,2012, 13(3): 523–538.

[72] Rosenbaum PR, Rubin DB. Reducing bias in observational studies using subclassification on the propensity score.,1984, 79(387): 516–524.

[73] Law CW, Chen Y, Wei S, Smyth GK. voom: precision weights unlock linear model analysis tools for RNA-seq read counts.,2014, 15(2): R29.

[74] Li J, Tibshirani R. Finding consistent patterns: A nonparametric approach for identifying differential expression in RNA-Seq data., 2011, 22(5): 519–536.

[75] Mcdaneld TG, Smith TP, Doumit ME, Miles JR, Coutinho LL, Sonstegard TS, Matukumalli LK, Dan JN, Wiedmann RT. MicroRNA transcriptome profiles during swine skeletal muscle development.,2009, 10: 77.

[76] Shen YF. Sequencing and characterization of mRNA, lncRNA, and miRNA in thytoid gland of Yorkshire and Jinhua Pigs[D]. Zhejiang University, 2016.沈一飛. 約克夏豬和金華豬甲狀腺組織mRNA、lncRNA和miRNA測(cè)序及功能分析[學(xué)位論文]. 浙江大學(xué), 2016.

[77] Li M, Xia Y, Gu Y, Zhang K, Lang Q, Chen L, Guan J, Luo Z, Chen H, Li Y, Li Q, Li X, Jiang AA, Shuai S, Wang J, Zhu Q, Zhou X, Gao X, Li X. MicroRNAome of porcine pre- and postnatal development.,2010, 5(7): e11541.

[78] Sun J, Xie M, Huang Z, Li H, Chen T, Sun R, Wang J, Xi Q, Wu T, Zhang Y. Integrated analysis of non-coding RNA and mRNA expression profiles of 2 pig breeds differing in muscle traits.,2017, 95(3): 1092–1103.

[79] Jianning YU, Yan L, Chen Z, Hui LI, Ying S, Zhu H, Shi Z. Investigating right ovary degeneration in chick embryos by transcriptome sequencing., 2017, 63(3): 295–303.

[80] Li T, Wang S, Wu R, Zhou X, Zhu D, Zhang Y. Identification of long non-protein coding RNAs in chicken skeletal muscle using next generation sequencing.,2012, 99(5): 292–298.

[81] Glazov EA, Cottee PA, Barris WC, Moore RJ, Dalrymple BP, Tizard ML. A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach., 18(6): 957–964.

[82] Kang B, Guo JR, Yang HM, Zhou RJ, Liu JX, Li SZ, Dong CY. Differential expression profiling of ovarian genes in prelaying and laying geese.,2009, 88(9): 1975–1983.

[83] Guo J, Tang QP, Zhang SJ, Ma YH, Lu HL. Identification of broodiness-related geese genes by suppression subtractive hybridization.,2011, 42(10): 1477–1784.郭軍, 湯青萍, 章雙杰, 馬月輝, 陸火林, 蘇建東, 鄒劍敏, 陳寬維, 李慧芳. 利用抑制消減雜交技術(shù)篩選鵝就巢行為相關(guān)基因. 畜牧獸醫(yī)學(xué)報(bào),2011, 42(10): 1477– 1484.

[84] Xu Q, Zhao WM, Chen Y, Tong YY, Rong GH, Huang ZY, Zhang Y, Chang GB, Wu XS, Chen GH. Transcriptome profiling of the goose () ovaries identify laying and broodiness phenotypes.,2013, 8(2): e55496.

[85] Chen L, Luo J, Li JX, Li JJ, Wang DQ, Tian Y, Lu LZ. Transcriptome analysis of adiposity in domestic ducks by transcriptomic comparison with their wild counterparts.,2015, 46(3): 299–307.

[86] Di R, He J, Song S, Tian D, Liu Q, Liang X, Ma Q, Sun M, Wang J, Zhao W, Cao G, Wang J, Yang Z, Ge Y, Chu M. Characterization and comparative profiling of ovarian microRNAs during ovine anestrus and the breeding season.,2014, 15: 899.

[87] Chang W, Wang J, Tao D, Zhang Y, Jianzhong HE, Shi C. Identification of a novel miRNA from the ovine ovary by a combinatorial approach of bioinformatics and experiments.,2015, 77(12): 1617–1624.

[88] Zhang SF, Wei CH, Lu J, Zhang XN, Zhou XL, Zhang SZ, Wang GK, Cao JX, Zhao FP, Zhang L, Du LX. Identification of the microRNAome in texel sheep by deep sequencing., 2013, 40(9): 19–22.張世芳, 魏彩虹, 陸健, 張小寧, 周鑫磊, 張淑珍, 王光凱, 曹家雪, 趙福平, 張莉, 杜立新. 深度測(cè)序鑒定綿羊microRNA轉(zhuǎn)錄組. 中國(guó)畜牧獸醫(yī), 2013, 40(9): 19– 22.

[89] Miao X, Qin QL. Genome-wide transcriptome analysis of mRNAs and microRNAs in Dorset and Small Tail Han sheep to explore the regulation of fecundity.,2015, 402: 32–42.

[90] Billerey C, Boussaha M, Esquerré D, Rebours E, Djari A, Meersseman C, Klopp C, Gautheret D, Rocha D. Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing.,2014, 15: 499.

[91] Sun X, Li M, Sun Y, Cai H, Li R, Wei X, Lan X, Huang Y, Lei C, Chen H. The developmental transcriptome landscape of bovine skeletal muscle defined by Ribo-Zero ribonucleic acid sequencing.2015, 93(12): 5648–5658.

[92] Cánovas A, Rincon G, IslasTrejo A, Wickramasinghe S, Medrano JF. SNP discovery in the bovine milk transcriptome using RNA-Seq technology.,2010, 21(11–12): 592–598.

[93] Yao B, Yu Z, Mei Z, Liu M, Liu H, Li J. De novo characterization of the antler tip of Chinese Sika deer transcriptome and analysis of gene expression related to rapid growth.,2012, 364(1–2): 93–100.

[94] Xie FY, Feng YL, Wang HH, Ma YF, Yang Y, Wang YC, Shen W, Pan QJ, Yin S, Sun YJ, Ma JY.assembly of the donkey white blood cell transcriptome and a comparative analysis of phenotype-associated genes between donkeys and horses.,2015, 10(7): e0133258.

[95] Scott EY, Mansour T, Bellone RR, Brown CT, Mienaltowski MJ, Penedo MC, Ross PJ, Valberg SJ, Murray JD, Finno CJ. Identification of long non-coding RNA in the horse transcriptome.,2017, 18(1): 511.

[96] Park KD, Park J, Ko J, Kim BC, Kim HS, Ahn K, Do KT, Choi H, Kim HM, Song S, Lee S, Jho S, Kong HS, Yang YM, Jhun BH, Kim C, Kim TH, Hwang S, Bhak J, Lee HK, Cho BW. Whole transcriptome analyses of six thoroughbred horses before and after exercise using RNA-Seq.,2012, 13: 473.

[97] Niu XY, Ren KL, Cao L, Li YP, Zhen JJ, Feng GL, Huang SF. Transcriptome analysis on coat color related genes in rex rabbits.,2016, 36(2): 1–6.牛曉艷, 任克良, 曹亮, 李燕平, 鄭建婷, 馮國(guó)亮, 黃淑芳. 利用轉(zhuǎn)錄組測(cè)序方法研究獺兔毛色相關(guān)基因. 中國(guó)草食動(dòng)物科學(xué), 2016, 36(2): 1–6.

[98] Song XC. Study on the mechanisms of hair pigmentation and skin transcriptome annotation in mink (neovison vison) based onhigh throughput RNA-sequencing[D]. Chinese Acadamy of Agricultural Sciences, 2016. 宋興超. 水貂被毛色素沉積機(jī)理及基于高通量RNA- seq皮膚轉(zhuǎn)錄組注釋研究[學(xué)位論文]. 中國(guó)農(nóng)業(yè)科學(xué)院, 2016.

Application of whole transcriptomics in animal husbandry

Tianpei Shi, Li Zhang

RNA is a polymeric molecule which is involved in various biological processes including the coding, decoding, regulation, and expression of genes. Whole transcriptome sequencing is the dominant method for studying RNA functions which assists researchers to deepen the exploration and analysis of gene structure and function and to reveal intrinsic links between gene expression and life phenomena. To date, extensive research has been done in animal husbandry models including swine and poultry by using high-throughput RNA sequencing technology. These studies included transcript structure, gene expression level and non-coding region function. In this review, we briefly introduce the general processes and strategies of RNA-seq and small RNA-seq technologies, and summarize the various achievements of the application of whole transcriptomics in animal husbandry.

gene expression; whole transcriptome; RNA-seq; high-throughput RNA sequencing

2018-07-31;

2018-10-22

國(guó)家自然科學(xué)基金項(xiàng)目(編號(hào)：U1503285) 和中國(guó)農(nóng)業(yè)科學(xué)院基本科研業(yè)務(wù)費(fèi)重大項(xiàng)目(編號(hào)：Y2017XM02)資助[Supported by the National Natural Science Foundation of China(No. U1503285) and Chinese Academy of Agricultural Sciences(No.Y2017XM02 )]

石田培，在讀碩士研究生，專業(yè)方向：動(dòng)物遺傳育種與繁殖。E-mail: 1337684764@qq.com

張莉，博士，研究員，博士生導(dǎo)師，研究方向：動(dòng)物遺傳育種與繁殖。E-mail: zhangli07@caas.cn

10.16288/j.yczz.18-218

2018/11/15 9:28:51

URI: http://kns.cnki.net/kcms/detail/11.1913.R.20181115.0928.002.html

(責(zé)任編委: 趙要風(fēng))